【发布时间】：2021-01-17 19:28:16
【问题描述】：
我试图了解如何在 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, y, 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], [0, 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, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0, 0, 0, 0, 0, 0, 1.0,
0.9919755553772727, 0.0, -0.0018716577540106951,
0.0, -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,
rtol=1e-13,atol=1e-22)
new_time = states.t
new_temp = states.y[:,-1]
index2 = [41, 40, 39, 38, 37, 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，dense_out = true enter image description here
但是这种方法也不能提供正确的结果。 以下是该方法的结果： enter image description here
更新：当我使用 MATLAB 使用的第一个时间步在纸上手动计算 RK45 时，我得到了相同的答案。此外，当我强制时间序列使用第一个时间间隔时，我会得到与 solve_ivp->RK45 相同的答案，但结果是密集的。然而，即使使用 MATLAB 中的相同完整时间序列，我得到的结果也与 MATLAB 不同。
@Lutz Lehmann 在对各种不同的方法进行了一些研究和测试之后，您认为 r.integrate 只集成一次是正确的。为了在每个点上进行积分，需要一个循环。此外，我能够让 ode 和 solve_ivp 得到相同的答案（尽管它是错误的答案）。使用 solve_ivp 时，我必须执行以下操作，这在使用 ode 时给出了相同的答案。
r = integrate.solve_ivp(fun=lambda t, y: function(t, y, mu),
t_span=time, y0=y, method='RK45', dense_output=True,
rtol=1e-13, atol=1e-22)
i = 0
while r.t[i] < tf:
r = integrate.solve_ivp(fun=lambda t, y: function(t, y, mu),
t_span=time, y0=y, method='RK45', dense_output=True,
rtol=1e-13, atol=1e-22)
print(r.t[i])
i += 1
new_time = r.t
new_temp = r.y[:, -1]
index2 = [41, 40, 39, 38, 37, 36]
new_temp = np.delete(new_temp, index2)
print(new_temp)
r = ode(function)
r.set_integrator('vode', method='bdf', rtol=1e-13, atol=1e-22, 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, 40, 39, 38, 37, 36]
new_temp = np.delete(new_temp, index2)
test = new_temp.reshape(6,6)
print(test)
我应该注意到，与使用 ode 相比，使用 solve_ivp 的方法要慢得多。产生相同结果的速度差异可能意味着 ode 是首选方法（不确定）。
这就是我得到的解决方案。 enter image description here
不幸的是，这意味着根据您上次发布的最新更新，我回到了我开始的地方。 ODE 和 solve_ivp 提供了相同的答案，但这仍然不是解决方案。
【问题讨论】：
-
您可以轻松调试 ODE 函数——传递相同的输入并查看它们是否从 Matlab 到 Python 产生相同的输出。如果可行，请检查 Ode 求解器——检查您是否使用相同的方法，您的语法是否正确。不过很难回答您的问题，因为您的帖子格式不太好（格式不一致，空格很多，难以阅读）。
-
Python 也使用自适应时间步进。就像您可以在 Matlab 中的指定时间强制评估一样。快速使用情况的主要区别在于
solve_ivp->RK45没有Refine选项，如ode45（默认值为4，附加点的插值）。因此，必须使用密集输出来填补空白以绘制合理的图表。 -
@David 在调试函数时，初始时间步长确实会产生相同的结果，但是错误的是 ode 求解器。我不知道如何更正解决方案以产生相同的结果。对不起，当我复制并通过代码表单蜘蛛时，它似乎出现了格式，引入了很多空格。如果您删除空格并更正格式，它将在任何给定的 python ide 中运行。
-
@LutzLehmann 我确实尝试了具有密集输出的solve_ivp -> rk45，不幸的是，与所提供的方法和matlab ode45相比，它产生的结果误差更大。我担心自适应步长大小不一样，并且由于 ode 的刚度，它会产生不正确的结果。如果您愿意，我可以发布使用 rk45 和 dense_output = True 的 solve_ivp 的结果。
-
@LutzLehmann 我无法在 cmets 中发布您建议的方法，因此我将其添加到原始帖子中。现在可以在问题的底部看到方法和输出。
标签： python matlab scipy ode orbital-mechanics