使用 rk4 绘制三个物体的轨道

Question

我一直在尝试使用 RK4 方法绘制三个粒子的轨迹。我无法在这段时间内生成一系列结果，因为它会显示以下错误消息：

  File "C:\Users\Local\Runge-Kutta 4 Code.py", line 65, in <module>
    solution.step()

  File "F:\Anaconda3\lib\site-packages\scipy\integrate\_ivp\base.py", line 170, in step
    raise RuntimeError("Attempt to step on a failed or finished "

RuntimeError: Attempt to step on a failed or finished solver.

我怀疑我的初始“y_0”值有问题，但我可能是错的。任何帮助都将不胜感激。我的代码如下：

import numpy as np
import matplotlib.pyplot as plt
import scipy as scipy

from scipy import integrate
from numpy import asarray
from numpy import savetxt

# Physical constants
mass_vector = np.array([1, 1, 1])

r_vec_1 = np.array([0, 0])
v_vec_1 = np.array([-np.sqrt(2), -np.sqrt(2)])

r_vec_2 = np.array([-1, 0])
v_vec_2 = np.array([np.sqrt(2) / 2, np.sqrt(2) / 2])

r_vec_3 = np.array([1, 0])
v_vec_3 = np.array([np.sqrt(2) / 2, np.sqrt(2) / 2])


#Initial x acceleration ODE's
def x1_double_dot(y, mass_vector):
    #G can be omitted for scale purposes (should not be compared with realistic data)
    return ((mass_vector[1]*(y[4]-y[0])/((y[0]-y[4])**2 + (y[1]-y[5])**2)**(3/2)) +
            (mass_vector[2]*(y[8]-y[0])/((y[0]-y[8])**2 + (y[1]-y[9])**2)**(3/2)))
def x2_double_dot(y, mass_vector):
    return ((mass_vector[0]*(y[0]-y[4])/((y[4]-y[0])**2 + (y[5]-y[1])**2)**(3/2)) +
            (mass_vector[2]*(y[8]-y[4])/((y[4]-y[8])**2 + (y[5]-y[9])**2)**(3/2)))
def x3_double_dot(y, mass_vector):
    return ((mass_vector[0]*(y[0]-y[8])/((y[8]-y[0])**2 + (y[9]-y[1])**2)**(3/2)) +
            (mass_vector[1]*(y[4]-y[8])/((y[8]-y[4])**2 + (y[9]-y[5])**2)**(3/2)))
#Initial y acceleration ODE's
def y1_double_dot(y, mass_vector):
    return ((mass_vector[1]*(y[5]-y[1])/((y[0]-y[4])**2 + (y[1]-y[5]**2)**(3/2))) +
            (mass_vector[2]*(y[9]-y[1])/((y[0]-y[8])**2 + (y[1]-y[9])**2)**(3/2)))
def y2_double_dot(y, mass_vector):
    return ((mass_vector[0]*(y[1]-y[5])/((y[4]-y[0])**2 + (y[5]-y[1]**2)**(3/2))) +
            (mass_vector[2]*(y[9]-y[5])/((y[4]-y[8])**2 + (y[5]-y[9])**2)**(3/2)))
def y3_double_dot(y, mass_vector):
    return ((mass_vector[0]*(y[1]-y[9])/((y[8]-y[0])**2 + (y[9]-y[1]**2)**(3/2))) +
             (mass_vector[1]*(y[5]-y[9])/((y[8]-y[4])**2 + (y[9]-y[5])**2)**(3/2)))

#This is my X(t) at time zero
y_0 = np.concatenate((r_vec_1, v_vec_1, r_vec_2, v_vec_2, r_vec_3, v_vec_3))
y = y_0

#This is my F(X) at time zero
def fun(t,y):
    return np.array([y[2], y[3], x1_double_dot(y, mass_vector), y1_double_dot(y, mass_vector),
              y[6], y[7], x2_double_dot(y, mass_vector), y2_double_dot(y, mass_vector),
              y[10], y[11], x3_double_dot(y, mass_vector), y3_double_dot(y, mass_vector)])

# collect data
t_values = []
y_values = []

#Time start, step, and finish point
t0,tf,t_step = 0, 2, 0.1
nsteps = int((tf - t0)/t_step)

solution = integrate.RK45(fun, t0, y_0, tf, first_step=t_step)

#The loop for running the Runge-Kutta method over some time period.
for step in range(nsteps):
    solution.step()
    y_values.append(solution.y[0])
    # break loop after modeling is finished
    if solution.status == 'finished':
        break

Answer 1

我将设置和导数计算压缩为

masses = [1,1,1]
r1,v1 = [ 0,0], [-2**0.5,-2**0.5]
r2,v2 = [-1,0], [0.5**0.5, 0.5**0.5]
r3,v3 = [ 1,0], [0.5**0.5, 0.5**0.5]
G = 1

def odesys(t,u):
    def force(a): return G*a/sum(a**2)**1.5
    r1,v1,r2,v2,r3,v3 = u.reshape([-1,2])
    m1,m2,m3 = masses
    f12, f13, f23 = force(r1-r2), force(r1-r3), force(r2-r3)
    a1,a2,a3 = -m2*f12-m3*f13, m1*f12-m3*f23, m1*f13+m2*f23
    return np.concatenate([v1,a1,v2,a2,v3,a3])

然后执行我基本上复制了你的代码，只是添加了一些对图形有用的选项

from scipy import integrate

#Time start, step, and finish point
t0,tf,t_step = 0, 2, 0.1
nsteps = int((tf - t0)/t_step)
u0 = np.concatenate([r1,v1,r2,v2,r3,v3])

solution = integrate.RK45(odesys, t0, u0, tf, first_step=0.2*t_step, max_step=t_step)

# collect data
t_values = [t0]
u_values = [u0]

#The loop for running the Runge-Kutta method over some time period.
for step in range(nsteps):
    solution.step()
    t_values.append(solution.t)
    u_values.append(solution.y)
    # break loop after modeling is finished
    if solution.status == 'finished':
        break

没有报告错误，解决方案绘制为

我的scipy版本是1.4.1。

剧情通过

获得

u = np.asarray(u_values).T
# x1,y1 = u[0], u[1],  x2,y2 = u[4],u[5]
plt.plot(u[0],u[1],'-o',lw=1, ms=3, label="body 1")
plt.plot(u[4],u[5],'-x',lw=1, ms=3, label="body 2")
plt.plot(u[8],u[9],'-s',lw=1, ms=3, label="body 3")

而不是，例如 u[4] 在转换 u_values 中的数据视图后，也可以使用 u_values[:][4] 使用原始结果数据结构。

随着数据的改变，质心基本保持不变，第一个天体比另外两个小，重力常数增加

masses = [0.01, 1, 1]
r1,v1 = [ 0,0], [-2**0.5,-2**0.5]
r2,v2 = [-1,0], [ 0.5**0.5, 0.5**0.5]
r3,v3 = [ 1,0], [-0.5**0.5,-0.5**0.5]
G = 4

生成的动态是

使用 rk4 绘制三个物体的轨道

Plot orbit of three bodies using rk4

python

scipy

runge-kutta

orbital-mechanics