I am trying to understand how to get the same results in Python as I get in MATLAB. Attached is the source code of what I have tried, results being incorrect for the two different methods. At the bottom of the code is the expected solution as the result of MATLAB. Any help with this problem would be greatly appreciated.
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
driver that will integrate the ODE(s):
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)
desired solution // MATLAB ode45 & ode113 same result
This is part of a greater series of scripts that I am writing and would prefer not have my code in MATLAB. I know the MATLAB answer is correct because the the end solutions provides the desired orbit I am trying to create. I should also note that it would appear that MATLAB is using adaptive step sizes and not a predefined time series one would create like in Python np.linspace(start,end,step)
A suggested method was the ivp_solver rk45 with dense_out = true enter image description here
however this method also does not provide the correct results. here are the results to that method: enter image description here
Update: When I manually calculate RK45 on paper with the first time step used by MATLAB I get the same answer. Also, when I force the time series to use the first time interval I get the same answer with the solve_ivp->RK45 with dense out. However even when using the same full time series from MATLAB I get results different from MATLAB.
@Lutz Lehmann After doing some research and testing of a variety of different methods you are correct in that r.integrate only integrate once. In order to integrate at each point a loop is required. Additionally, I was able to get ode and solve_ivp to same answer (although it is the wrong answer). When using solve_ivp I had to do the following which gave me the same answer when using 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)
I should note that the method using solve_ivp is much slower compared to using ode. The difference in speed yielding the same result probably means that ode is the preferred method (not sure).
This was the solution I got. enter image description here
Unfortunately what means is that based on this latest update conducted from your last post I am back to where I started. ODE and solve_ivp provide the same answer however this is still not the solution.
solve_ivp call you defined the function wrong, the correct way isfun=lambda t, y: function(t, y, mu)
ode solver in the function r.integrate seems to perform at most one internal step. To reach the final time tf, you have to loop this call:while r.t < tf: r.integrate(tf)
The results of both then are sufficiently coincident, in the leading 10 or so digits. For the source to the differences to the Matlab results see the last section.
PS: You can drastically reduce the second function, all the flattening and copying is not really necessary. I rewrote it as
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
In the computation of the Hessean of the potential in function2 there were many small errors in indices and parentheses placement. Re-organized and with more structuring variables the function can look like
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
This gives the results
[ 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]
which as far as I can see coincide with the Matlab results.