python中四阶Runge-Kutta法求解3个耦合非线性微分方程

2024-04-18 08:54:45 发布

您现在位置：Python中文网/ 问答频道 /正文
6397 3

我试图绘制一个带电粒子在Reissner–Nordström黑洞（带电黑洞）周围的轨道

我有三个二阶微分方程和三个一阶微分方程。由于问题的性质，每个导数都以适当的时间而不是时间t表示。运动方程如下

2 first order differential equation second order differential equations

3 second order differential equations

1 first order differential equation(there should be a negative multiplied by everything under the square root.

我用四阶龙格库塔方法来积分轨道。我的困惑，以及我最可能犯错误的地方，来自这样一个事实，通常当你有一个二阶耦合微分方程，你把它简化成两个一阶微分方程。然而，在我的问题中，我得到了3个一阶微分方程及其相应的二阶微分方程。我假设，既然给了我这些一阶方程，我就不必减少二阶。这些方程是非线性的，这一事实使事情更加复杂

我相信我可以用龙格库塔来解决这些问题，但是我不确定我的运动方程的实现。当我运行代码时，我得到一个错误，一个负数在F2的平方根之下，但是情况不应该是这样，因为F2应该正好等于零（毫无疑问，这是由F1引起的精度问题）。然而，即使我取所有的绝对值，在F1，F2，F3的平方根下。。。我的角动量L和能量E是不守恒的。我主要希望有人对我在龙格-库塔循环中使用微分方程的方式进行评论，并告诉我应该如何简化二阶微分方程

import matplotlib.pyplot as plt
import numpy as np
import math as math
#=============================================================================
h=1
M         = 1                  #Mass of RN blackhole
r         = 3*M                #initital radius of particle from black hole
Q         = 0                  #charge of particle
r_s       = 2*M                #Shwar radius
S         = 0                  # initial condition for RK4
V         = .5                 # Initial total velocity of particle
B         = np.pi/2            #angle of initial velocity
V_p       = V*np.cos(B)        #parallel velocity
V_t       = V*np.sin(B)        #transverse velocity
t         = 0
Theta     = 0
E         = np.sqrt(Q**2-2*r*M+r**2)/(r*np.sqrt(1-V**2))
L         = V_t*r/(np.sqrt(1-V**2))
r_dot     = V_p*np.sqrt(r**2-2*M+Q**2)/(r*np.sqrt(1-V**2))
Theta_dot = V_t/(r*np.sqrt(1-V**2))
t_dot     = E*r**2/(r**2-2*M*r+Q**2)

#=============================================================================
while(r>2*M and r<10*M):   #Runge kutta while loop
    A1 = 2*(Q**2-M*r) * r_dot*t_dot / (r**2-2*M*r+Q**2)                                           #defines T double dot fro first RK4 step
    B1 = -2*Theta_dot*r_dot / r                                                                   #defines theta double dot for first RK4 step
    C1 = (r-2*M*r+Q**2)*(Q**2-M*r)*t_dot**2 / r**5 + (M*r-Q**2)*r_dot**2 / (r**2-2*M*r+Q**2)      #defines r double dot for first RK4 step
    D1 = E*r**2/(r**2-2*M*r+Q**2)                                                                 #defines T dot for first RK4 step
    E1 = L/r**2                                                                                   #defines theta dot for first RK4 step
    F1 = math.sqrt(-(1-r_s/r+Q**2/r**2) * (1-(1-r_s/r+Q**2/r**2)*D1**2 + r**2*E1**2))              #defines r dot for first RK4 step
    
    t_dot_1     = t_dot     + (h/2) * A1
    Theta_dot_1 = Theta_dot + (h/2) * B1
    r_dot_1     = r_dot     + (h/2) * C1
    t_1         = t         + (h/2) * D1
    Theta_1     = Theta     + (h/2) * E1
    r_1         = r         + (h/2) * F1
    S_1           = S         + (h/2)
    
    A2 = 2*(Q**2-M*r_1) * r_dot_1*t_dot_1 / (r_1**2-2*M*r_1+Q**2)                                                   
    B2 = -2*Theta_dot_1*r_dot_1 / r_1                                                                               
    C2 = (r_1-2*M*r_1+Q**2)*(Q**2-M*r_1)*t_dot_1**2 / r_1**5 + (M*r_1-Q**2)*r_dot_1**2 / (r_1**2-2*M*r_1+Q**2)      
    D2 = E*r_1**2/(r_1**2-2*M*r_1+Q**2)                                                                                   
    E2 = L/r_1**2                                                                                                     
    F2 = np.sqrt(-(1-r_s/r_1+Q**2/r_1**2) * (1-(1-r_s/r_1+Q**2/r_1**2)*D2**2 + r_1**2*E2**2))                                 
    
    t_dot_2     = t_dot     + (h/2) * A2
    Theta_dot_2 = Theta_dot + (h/2) * B2
    r_dot_2     = r_dot     + (h/2) * C2
    t_2         = t         + (h/2) * D2
    Theta_2     = Theta     + (h/2) * E2
    r_2         = r         + (h/2) * F2
    S_2           = S         + (h/2)
    
    
    A3 = 2*(Q**2-M*r_2) * r_dot_2*t_dot_2 / (r_2**2-2*M*r_2+Q**2)                                                   
    B3 = -2*Theta_dot_2*r_dot_2 / r_2                                                                               
    C3 = (r_2-2*M*r_2+Q**2)*(Q**2-M*r_2)*t_dot_2**2 / r_2**5 + (M*r_2-Q**2)*r_dot_2**2 / (r_2**2-2*M*r_2+Q**2)      
    D3 = E*r_2**2/(r_2**2-2*M*r_2+Q**2)                                                                                   
    E3 = L/r_2**2                                                                                                     
    F3 = np.sqrt(-(1-r_s/r_2+Q**2/r_2**2) * (1-(1-r_s/r_2+Q**2/r_2**2)*D3**2 + r_2**2*E3**2))                                 
    
    t_dot_3     = t_dot     + (h/2) * A3
    Theta_dot_3 = Theta_dot + (h/2) * B3
    r_dot_3     = r_dot     + (h/2) * C3
    t_3         = t         + (h/2) * D3
    Theta_3     = Theta     + (h/2) * E3
    r_3         = r         + (h/2) * F3 
    S_3           = S       + (h/2)
    
    
    
    A4 = 2*(Q**2-M*r_3) * r_dot_3*t_dot_3 / (r_3**2-2*M*r_3+Q**2)                                                   
    B4 = -2*Theta_dot_3*r_dot_3 / r_3                                                                              
    C4 = (r_3-2*M*r_3+Q**2)*(Q**2-M*r_3)*t_dot_3**2 / r_3**5 + (M*r_3-Q**2)*r_dot_3**2 / (r_3**2-2*M*r_3+Q**2)      
    D4 = E*r_3**2/(r_3**2-2*M*r_3+Q**2)                                                                                   
    E4 = L/r_3**2                                                                                                     
    F4 = np.sqrt(-(1-r_s/r_3+Q**2/r_3**2) * (1-(1-r_s/r_3+Q**2/r_3**2)*D3**2 + r_3**2*E3**2))                                #defines r dot for first RK4 step
    
    t_dot     = t_dot     + (h/6.0) * (A1+(2.*A2)+(2.0*A3) + A4)
    Theta_dot = Theta_dot + (h/6.0) * (B1+(2.*B2)+(2.0*B3) + B4)
    r_dot     =  r_dot    + (h/6.0) * (C1+(2.*C2)+(2.0*C3) + C4)
    t         =  t        + (h/6.0) * (D1+(2.*D2)+(2.0*D3) + D4)
    Theta     = Theta     + (h/6.0) * (E1+(2.*E2)+(2.0*E3) + E4)
    r         =  r        + (h/6.0) * (F1+(2.*F2)+(2.0*F3) + F4)
    S          = S+h
    
    print(L,r**2*Theta_dot)
    
    
    plt.axes(projection = 'polar')
    plt.polar(Theta, r, 'g.')

Tags： offorstepnpsqrtdotf2f1
1条回答
网友
1楼 · 发布于 2024-04-18 08:54:45

以您提供的三个二阶微分方程为例。这些是经过适当时间参数化的测地线方程。但是，原始度量是旋转不变的（即SO（3）不变），因此它有一组简单的守恒定律，加上度量的守恒（即适当时间的守恒）。这意味着t和theta的二阶微分方程可以一次积分，从而得到t和theta的两个一阶微分方程组和r的一个二阶微分方程组：

dt/ds = c_0 * r**2 / (r**2 - 2*M*r + Q**2)

dtheta/ds = c_1 / r**2

d**2r/ds**2 = ( (r**2-2*M*r + Q**2)*(Q**2 - M*r)/r**5) * (dt/ds)**2 
                   + ( (M*r - Q**2) /(r**2 - 2*M*r + Q**2) ) * (dr/ds)**2

这里你可以用不同的方法，其中一种方法是通过将上面的前两个方程代入轨道上的度量等于1的方程，来推导r的一阶运动微分方程的方程。但是你也可以直接到这里，把dt/ds方程的右边插入r的第三个方程，将系统表示为

dt/ds = c_0 * r**2 / (r**2 - 2*M*r + Q**2)

dtheta/ds = c_1 / r**2

d**2r/ds**2 = ( c_0**2*(Q**2 - M*r)/(r*(r**2-2*M*r + Q**2)))
                + ( (M*r - Q**2) /(r**2 - 2*M*r + Q**2) ) * (dr/ds)**2

为了避免使用平方根和复杂度（平方根也是昂贵的计算，而有理函数是简单且快速的代数计算），定义四个一阶微分方程的等效系统

dt/ds = c_0 * r**2 / (r**2 - 2*M*r + Q**2)

dtheta/ds = c_1 / r**2

dr/ds = u

du/ds = ( c_0**2*(Q**2 - M*r)/(r*(r**2-2*M*r + Q**2)))
                + ( (M*r - Q**2) /(r**2 - 2*M*r + Q**2) ) * u**2

借助t, theta, r的初始条件及其导数dt/dt, dtheta/dt, dr/dt，您可以计算第一个和第二个方程中使用的常数c_0和c_1，然后计算u = dr/dt的初始条件

相关问题 更多 >

    编程相关推荐