0

I want to integrate the same differential equation with different initial conditions and I want to compare each points that I have obtained. Unfortunately if I do:

sol = solve_ivp(my_model, [0, leng], [1,1,1], method="LSODA", args=(rho, sigma, beta), dense_output=False)

sol_2 = solve_ivp(my_model, [0, leng], [1+0.001,1+0.001, 1-0.001], method="LSODA", args=(rho, sigma, beta), dense_output=False)

The number of points obtained in these two integration processes is different.

I don't want to use the interpolation procedure since I would like to work with real data.

Can I set the number of points of the solution(s)?

Siderius
  • 174
  • 2
  • 14
  • May I ask for "my_model" and leng? – Rafael Valero Feb 06 '21 at 21:43
  • Is your doubt on the interpolated values related to the behavior observed in https://stackoverflow.com/questions/21859870/absolute-error-of-ode45-and-runge-kutta-methods-compared-with-analytical-sol? – Lutz Lehmann Feb 06 '21 at 21:46
  • 2
    You could put both instances of the system into the same ODE function, so that you integrate a 6 dimensional state. Then you get the step sizes minimized across the instances, but the values at the same time. So in total with slightly less error than in the individual integration. – Lutz Lehmann Feb 06 '21 at 21:54
  • @RafaelValero : The test model is, by parameter names and dimension, quite probably the Lorenz attractor. With the given differences, the solutions should visibly separate rather quickly, latest at leng=20, possibly already with leng=5. So to observe non-chaotic, regular error behavior, one would have to set `leng=2`. – Lutz Lehmann Feb 06 '21 at 21:58
  • @RafaelValero yes the separation is ok but I want two trajectories with the same amount of points. It is not strange to obtain a difference in the number of points roughly equal to 10^4. – Siderius Feb 06 '21 at 22:04

1 Answers1

1

One variant to compare instances differing in initial conditions is to collate them all into one big system where they are all solved simultaneously all with the same step sequence

def multi_Lorenz(t,u):
    return np.concatenate([Lorenz(t,uu) for uu in u.reshape([-1,3])])

u0 = np.concatenate([[1+k*1e-5,1+k*1e-5,1-k*1e-5] for k in range(11)])
res = solve_ivp(multi_Lorenz,[1,25],u0, method="LSODA")

plt.figure(figsize=(14,6))
plt.plot(res.t, res.y[::3].T)
plt.grid(); plt.show()

plot of the solutions

The visible split of the solution family is at t=22, however in the differences to the k=0 graph as reference, the increase to the tenfold of the initial difference happens almost immediately in the spike around t=1 and remains in that range up to t=12

N = max(k for k,t in enumerate(res.t) if t<12)
plt.figure(figsize=(14,6))
plt.plot(res.t[:N], (res.y[3::3,:N]-res.y[0,:N]).T)
plt.grid(); plt.show()

enter image description here

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51