计算相同的结果 MATLAB ODE45 和 Python
Computing the Same Results MATLAB ODE45 and Python
我试图了解如何在 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 相同的结果
enter image description here
这是我正在编写的一系列脚本中的一部分，我不希望将我的代码放在 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
不幸的是，根据您上次 post 进行的最新更新，我又回到了起点。 ODE 和 solve_ivp 提供了相同的答案，但这仍然不是解决方案。
- 在
solve_ivp调用中你定义的函数有误，正确的方法是
fun=lambda t, y: function(t, y, mu)
- 函数
r.integrate 中的 ode 求解器似乎最多执行一个内部步骤。要到达最终时间 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], [0, 0, 0]])
Df = np.block([[OxO, Ind],
[GMatrix, K]])
DfA = np.matmul(Df, A)
x,y,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 时，在索引和括号位置上存在许多小错误。 Re-organized 并且有更多的结构化变量，函数看起来像
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 结果一致。
我试图了解如何在 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 相同的结果
enter image description here
这是我正在编写的一系列脚本中的一部分，我不希望将我的代码放在 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
不幸的是，根据您上次 post 进行的最新更新，我又回到了起点。 ODE 和 solve_ivp 提供了相同的答案，但这仍然不是解决方案。
- 在
solve_ivp调用中你定义的函数有误，正确的方法是
fun=lambda t, y: function(t, y, mu)
- 函数
r.integrate中的ode求解器似乎最多执行一个内部步骤。要到达最终时间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], [0, 0, 0]])
Df = np.block([[OxO, Ind],
[GMatrix, K]])
DfA = np.matmul(Df, A)
x,y,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 时，在索引和括号位置上存在许多小错误。 Re-organized 并且有更多的结构化变量，函数看起来像
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 结果一致。