问题描述
我一直在尝试使用 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,y[6],y[7],x2_double_dot(y,y2_double_dot(y,y[10],y[11],x3_double_dot(y,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,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
解决方法
我将设置和导数计算浓缩为
masses = [1,1,1]
r1,v1 = [ 0,0],[-2**0.5,-2**0.5]
r2,v2 = [-1,[0.5**0.5,0.5**0.5]
r3,v3 = [ 1,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,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,v3])
solution = integrate.RK45(odesys,t0,u0,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],'-o',lw=1,ms=3,label="body 1")
plt.plot(u[4],u[5],'-x',label="body 2")
plt.plot(u[8],u[9],'-s',label="body 3")
代替，例如，u[4] 中的数据视图转换后的 u_values 也可以使用使用原始结果数据结构的 u_values[:][4]。
随着数据的变化，质心在很大程度上保持不变，第一个物体比其他两个小，重力常数增加
masses = [0.01,[ 0.5**0.5,[-0.5**0.5,-0.5**0.5]
G = 4
由此产生的动态是
