Я пытаюсь построить орбиты заряженной частицы вокруг черной дыры Рейсснера – Нордстрема (заряженной черной дыры).
У меня есть три дифференциальных уравнения 2-го порядка, а также 3 дифференциальных уравнения первого порядка. Из-за природы проблемы каждая производная выражается в терминах собственного времени, а не времени t. Уравнения движения следующие.
2 дифференциальных уравнения первого порядка, дифференциальные уравнения второго порядка
3 дифференциальных уравнения второго порядка
Я использую метод Рунге-Кутты 4-го порядка для интеграции орбит. Мое замешательство и то, в чем я, скорее всего, ошибаюсь, связано с тем, что обычно, когда у вас есть связанное дифференциальное уравнение второго порядка, вы сводите его к двум дифференциальным уравнениям первого порядка. Однако в моей задаче мне были даны 3 дифференциальных уравнения первого порядка вместе с соответствующими им дифференциальными уравнениями второго порядка. Я предположил, что, поскольку мне были даны эти уравнения первого порядка, мне вообще не пришлось бы уменьшать второй порядок. Тот факт, что эти уравнения нелинейны, еще больше усложняет ситуацию.
Я уверен, что могу использовать Runge kutta для решения таких задач, однако я не уверен в своей реализации уравнений движения. Когда я запускаю код, я получаю сообщение об ошибке, что отрицательное число меньше квадратного корня из F2, однако этого не должно быть, потому что F2 должен быть точно равен нулю (несомненно, проблема точности, возникающая из-за F1). Однако даже когда я беру абсолютное значение всего под квадратным корнем из F1, F2, F3 ... мой угловой момент L и энергия E не сохраняются. В основном я хотел бы, чтобы кто-то прокомментировал, как я использую свое дифференциальное уравнение внутри цикла Рунге-кутта, и проинформировал меня, как я должен уменьшить дифференциальные уравнения 2-го порядка.
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.')
r^2-2Mr+Q^2. Находятся ли полярные координатыrиthetaотносительно фиксированного центра или чего-то совершенно другого? Кажется, что первая система уравнений является частью вывода, показывающим, как скорость распределяется по полярным компонентам. Уравнения движения входят только во втором наборе. Третье, кажется, относится к некоему первому интегралу, обычно к энергии? Моделирование должно быть завершено только со вторым набором. - person Lutz Lehmann schedule 17.03.2021