问题描述
我试图了解如何在Python中获得与在MATLAB中相同的结果。附件是我尝试过的源代码，对于两种不同的方法，结果是不正确的。代码的底部是MATLAB的预期解决方案。任何解决此问题的帮助将不胜感激。
from scipy.integrate import ode
from scipy import integrate
import numpy as np
def function2(x,mu):
x,y,z = x
r1 = np.sqrt((x + mu) ** 2 + (y ** 2) + (z ** 2))
r2 = np.sqrt((x - (1 - mu)) ** 2 + (y ** 2) + (z ** 2))
u1_x = 1 - (1 - mu) * (1 / (r1 ** 3) - 3 * ((x + mu) ** 2) / (r1 ** 5)) - \
mu * (1 / (r2 ** 3) - 3 * ((x - (1 - mu)) ** 2) / (r2 ** 5))
u2_y = 1 - (1 - mu) * (1 / (r1 ** 3)) - 3 * y ** 2 / (r1 ** 5) - \
mu * (1 / r2 ** 3 - 3 * y ** 2 / r2 ** 5)
u3_z = (-1) * (1 - mu) * (1 / r1 ** 3) - 3 * z ** 2 / r1 ** 5 - mu * \
(1 / r2 ** 3 - 3 * z ** 2 / r2 ** 5)
u1_y = 3 * (1 - mu) * y * (x + mu) / r1 ** 5 + \
3 * mu * y * (z - (1 - mu)) / r2 ** 5
u1_z = 3 * (1 - mu) * z * (x + mu) / r1 ** 5 + \
3 * mu * z * (x - (1 - mu)) / r2 ** 5
u2_z = 3 * (1 - mu) * y * z / r1 ** 5 + 3 * mu * y * z / r2 ** 5
u3_y = u2_z
u2_x = u1_y
u3_x = u1_z
gmatrix = np.array([[u1_x,u1_y,u1_z],[u2_x,u2_y,u2_z],[u3_x,u3_y,u3_z]])
return gmatrix
def function(t,mu):
x = y[36:39]
GMatrix = function2(x,mu)
OxO = np.zeros([3,3])
Ind = np.identity(3)
K = np.array([[0,2,0],[-2,[0,0]])
Df = np.bmat([[OxO[0],Ind[0]],[OxO[1],Ind[1]],[OxO[2],Ind[2]],[GMatrix[0],K[0]],[GMatrix[1],K[1]],[GMatrix[2],K[2]]])
Df = np.reshape(Df,(6,6))
A_temp = np.squeeze(np.array(y))
A_temp = A_temp.flatten()
B_temp = [0]*42
for i in range(len(A_temp)):
B_temp[i] = A_temp[i]
B_temp = B_temp[:-6]
B_temp = np.array(B_temp)
A = B_temp.reshape(6,6)
DfA = np.matmul(Df,A)
a = [0] * 36
b = np.squeeze(np.array(DfA))
b = b.flatten()
for i in range(len(b)):
a[i] = b[i]
r1 = np.sqrt((mu+y[36])**2 + (y[37]**2) + (y[38]**2))
r2 = np.sqrt((1-mu-y[36])**2 + (y[37]**2) + (y[38]**2))
m1 = 1 - mu
m2 = mu
c = [y[39],y[40],y[41],y[36] + 2 * y[40] + m1 * (-mu - y[36]) / (r1**3) + m2 * (1-mu-y[36]) / (r2**3),y[37] - 2 * y[39] - m1 * (y[37]) / (r1**3) - m2 * y[37] / (r2**3),-m1 * y[38] / (r1**3) - m2 * y[38] / (r2**3)]
ydot = a + c
return ydot
将集成ODE的driver：
if __name__ == '__main__':
t0 = 0
tf = 1.450000000000000
mu = 3.054248395728148e-06
x_n = [1.0,0.0,1.0,0.9919755553772727,-0.0018716577540106951,-0.0117506137115032,0.0]
#meth = 'adams'
meth = 'bdf'
r = ode(function).set_integrator('vode',method=meth,rtol=1e-13,atol=1e-22,with_jacobian=False)
r.set_initial_value(x_n,t0).set_f_params(mu)
r.integrate(tf)
temp = r.y
index2 = [41,40,39,38,37,36]
temp = np.delete(temp,index2)
temp = temp.reshape(6,6)
time = [t0,tf]
states = integrate.solve_ivp(fun=lambda t,y:
function(t,x_n,mu),t_span=time,y0=x_n,method='LSODA',dense_output=True,atol=1e-22)
new_time = states.t
new_temp = states.y[:,-1]
index2 = [41,36]
new_temp = np.delete(new_temp,index2)
new_temp = new_temp.reshape(6,6)
print(new_temp)
print(temp)
所需解决方案// MATLAB ode45和ode113相同的结果
这是我正在编写的一系列脚本的一部分，并且不希望我的代码不在MATLAB中。我知道MATLAB的答案是正确的，因为最终解决方案提供了我尝试创建的所需轨道。我还要注意，似乎MATLAB正在使用自适应步长，而不是像Python np.linspace(start,end,step)
建议的方法是ivp_solver rk45，其中density_out = true enter image description here
但是，此方法也无法提供正确的结果。 这是该方法的结果： enter image description here
更新：当我在MATLAB所使用的第一步下手动计算纸上的RK45时，我得到了相同的答案。另外，当我强迫时间序列使用第一个时间间隔时，我得到与solve_ivp-> RK45相同的答案，且密集。但是，即使使用来自MATLAB的相同的全时序列，我也会得到与MATLAB不同的结果。
@Lutz Lehmann在对各种不同的方法进行了研究和测试之后，您会发现正确的答案是r.integrate仅集成一次。为了在每个点集成，需要一个循环。另外，我能够将ode和resolve_ivp设为相同的答案（尽管这是错误的答案）。当使用resolve_ivp时，我必须执行以下操作，这在使用ode时给出了相同的答案。
r = integrate.solve_ivp(fun=lambda t,y: function(t,y0=y,method='RK45',atol=1e-22)
i = 0
while r.t[i] < tf:
r = integrate.solve_ivp(fun=lambda t,y: function(t,atol=1e-22)
print(r.t[i])
i += 1
new_time = r.t
new_temp = r.y[:,index2)
print(new_temp)
r = ode(function)
r.set_integrator('vode',method='bdf',with_jacobian=False)
r.set_initial_value(y,t0)
r.set_f_params(mu)
r.integrate(tf)
t = []
Y = [y]
while r.t < tf:
r.integrate(tf,step=True)
Y = np.vstack((Y,[r.y]))
t.append([r.t])
new_temp = Y[-1,:]
index2 = [41,index2)
test = new_temp.reshape(6,6)
print(test)
我应该注意，与使用ode相比，使用solve_ivp的方法要慢得多。产生相同结果的速度差异可能意味着ode是首选方法（不确定）。
这是我得到的解决方案。 enter image description here
不幸的是，这意味着根据您上次发布的最新信息，我回到了开始的位置。 ODE和solve_ivp提供了相同的答案，但这仍然不是解决方案。
解决方法
- 在
solve_ivp调用中，您定义了错误的函数，正确的方法是
fun=lambda t,y: function(t,y,mu)
- 函数
ode中的r.integrate求解器似乎最多执行一个内部步骤。要达到最终时间tf，您必须循环调用：
while r.t < tf: r.integrate(tf)
然后，两者的结果在前10位左右是足够一致的。有关与Matlab结果差异的来源，请参见最后一节。
PS：您可以大大减少第二个功能，实际上并不需要所有的拼合和复制。我改写为
def function(t,u,mu):
A = np.reshape(u[:36],[6,6])
x = u[36:39]
v = u[39:]
GMatrix = function2(x,mu)
OxO = np.zeros([3,3])
Ind = np.identity(3)
K = np.array([[0,2,0],[-2,[0,0]])
Df = np.block([[OxO,Ind],[GMatrix,K]])
DfA = np.matmul(Df,A)
x,z = x
vx,vy,vz = v
r1 = np.sqrt((mu+x)**2 + (y**2) + (z**2))
r2 = np.sqrt((1-mu-x)**2 + (y**2) + (z**2))
m1 = 1 - mu
m2 = mu
a = [x + 2 * vy - m1 * (mu + x) / (r1**3) + m2 * (1-mu-x) / (r2**3),y - 2 * vx - m1 * y / (r1**3) - m2 * y / (r2**3),-m1 * z / (r1**3) - m2 * z / (r2**3)]
ydot = np.concatenate([DfA.flatten(),v,a])
return ydot
在function2中势的Hessean的计算中，索引和括号位置存在许多小错误。重新组织并且具有更多结构变量，该函数看起来像
def function2(x,mu):
ce1,ce2 = -mu,1-mu
m1,m2 = 1-mu,mu
r1 = ((x[0]-ce1)**2+x[1]**2+x[2]**2)**0.5
r2 = ((x[0]-ce2)**2+x[1]**2+x[2]**2)**0.5
u1_x = 1 - m1 * (1 / r1**3 - 3 * (x[0] - ce1)**2 / r1**5) \
- m2 * (1 / r2**3 - 3 * (x[0] - ce2)**2 / r2**5)
u1_y = 3 * m1 * x[1] * (x[0] - ce1) / r1**5 \
+ 3 * m2 * x[1] * (x[0] - ce2) / r2**5
u1_z = 3 * m1 * x[2] * (x[0] - ce1) / r1**5 \
+ 3 * m2 * x[2] * (x[0] - ce2) / r2**5
u2_x = u1_y
u2_y = 1 - m1 * (1 / r1**3 - 3 * x[1]**2 / r1**5) \
- m2 * (1 / r2**3 - 3 * x[1]**2 / r2**5)
u2_z = 3 * m1 * x[1] * x[2] / r1**5 + 3 * m2 * x[1] * x[2] / r2**5
u3_x = u1_z
u3_y = u2_z
u3_z = - m1 * (1 / r1**3 - 3 * x[2]**2 / r1**5) \
- m2 * (1 / r2**3 - 3 * x[2]**2 / r2**5)
GMatrix = np.array([[u1_x,u1_y,u1_z],[u2_x,u2_y,u2_z],[u3_x,u3_y,u3_z]])
return GMatrix
这给出了结果
[ 23.55077315 -0.39182483 3.67078856 4.77188265 4.55322364 0.54862501]
[ -8.519012 -0.71609406 -1.5344374 -2.12546806 -1.4395364 -0.23934585]
[ -0.37941138 0.11874396 -1.39417439 -0.10224713 -0.13959515 0.19607223]
[ 43.56974185 -1.72339855 6.85563491 8.62759039 8.39374083 0.9221739 ]
[-28.39640391 -0.10433561 -4.47605353 -5.56490582 -6.06643015 -0.69034888]
据我所知与Matlab结果一致。
设置时间 控制面板
错误1：Request method ‘DELETE‘ not supported 错误还原：...
