# Instability detected in solving ODE with DifferentialEquations.jl

I am trying to solve an equation in Julia but I get Unstable retcode, I have solved the same eqaution in MATLAB with no issues. here is a simplified version of the code:

``````using DifferentialEquations

struct RosslerMSF
alpha::Float64
beta::Float64
gamma::Float64
dT::Float64
end

Jac(f, alpha, gamma) = [0 -1 -1; 1 alpha 0; f[3] 0 f[1]-gamma]

function (p::RosslerMSF)(du,u,par,t)
ind = floor(Int,t/p.dT)+1
J = Jac([x, y, z], p.alpha, p.gamma)

du[1:3] = J*[u[1], u[2], u[3]]
du[4:12] = J*[u[4] u[7] u[10];u[5] u[8] u[11];u[6] u[9] u[12]]
end

function Rossler(du,u,p,t)
x, y, z = u[1], u[2], u[3]

du[1] = -y-z
du[2] = x + 0.2*y
du[3] = 0.2 + (x - 9)*z
end

init = [0.1 -0.5 2]
sol = solve(ODEProblem(Rossler,init,(0,2200)),Tsit5(),dt=0.001,saveat = 200:0.001:2200)

funx = RosslerMSF(0.2,0.2,9,[sol[:,1,:]; sol[:,2,:]; sol[:,3,:]],0.001)
initx = [init..., 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
solx = solve(ODEProblem(funx, initx, (0, 2000)), Tsit5(), dt=0.001)
``````

What value of `u` did you use to verify the equations? Also, did you try MATLABDiffEq.jl to double check if the solution is the same?

For the first function `Rossler` both MATLAB and Julia work and their results agree, but for second equation `RosslerMSF` Julia fails to solve the equation while MATLAB solves it correctly.
I donâ€™t know much about MATLABDiffEq.jl, I tried to use it to solve the equations and for the `Rossler` function the results agreed in MATLAB and Julia, but for second function `RosslerMSF`, which causes the unstability problem I got so many errors and couldnâ€™t solve it using MATLABDiffEq.jl.

Thatâ€™s not going to give a well-defined solution at high accuracy because itâ€™s discontinuous. Youâ€™d want to set `tstops` to each discontinuity. You havenâ€™t set the tolerances so youâ€™d only expect around 2 digits to be the same by default, and then with the discontinuity all bets are kind of off. Try setting the tolerances of `abstol=1e-9,reltol=1e-9` on both and see what you get.

But also, I donâ€™t understand why youâ€™re doing it like that. If you instead just passed `sol` as a parameter and used `sol(t)` it would be continuous and more accurate. Basically:

``````funx = RosslerMSF(0.2,0.2,9,sol,0.001)
``````

and then:

``````function (p::RosslerMSF)(du,u,par,t)
J = Jac([x, y, z], p.alpha, p.gamma)

du[1:3] = J*[u[1], u[2], u[3]]
du[4:12] = J*[u[4] u[7] u[10];u[5] u[8] u[11];u[6] u[9] u[12]]
end
``````

would be simpler and be the same? Youâ€™d want to remove `saveat` in the first solve so youâ€™d get fully continuous derivatives rather than a linear interpolation.

1 Like

I applied the changes you mentioned, but I still get the same instability warning.
Here is the code after making those changes:

``````using DifferentialEquations

struct RosslerMSF
alpha::Float64
beta::Float64
gamma::Float64
end

Jac(f, alpha, gamma) = [0 -1 -1; 1 alpha 0; f[3] 0 f[1]-gamma]

function (p::RosslerMSF)(du,u,par,t)
J = Jac([x, y, z], p.alpha, p.gamma)

du[1:3] = J*[u[1], u[2], u[3]]
du[4:12] = J*[u[4] u[7] u[10];u[5] u[8] u[11];u[6] u[9] u[12]]
end

function Rossler(du,u,p,t)
x, y, z = u[1], u[2], u[3]

du[1] = -y-z
du[2] = x + 0.2*y
du[3] = 0.2 + (x - 9)*z
end

init = [0.1 -0.5 2]
sol = solve(ODEProblem(Rossler,init,(0,2200)),Tsit5(),abstol=1e-9,reltol=1e-9,dt=0.001)

funx = RosslerMSF(0.2,0.2,9,sol)
initx = [init..., 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
solx = solve(ODEProblem(funx, initx, (0, 2000)), Tsit5(), abstol=1e-9, reltol=1e-9, dt=0.001)

``````

What solver are you using in MATLAB? Can you include the MATLAB code?

Sure, here is the MATLAB code:

``````init = [0.1,-0.5,2];
[~,V]=ode45(@(~,X) [        -X(2)-X(3) ;
X(1)+0.2*X(2) ;
0.2+(X(1)-9)*X(3)], (0:0.001:2200),init);

X = V(200001:end,1);
Y = V(200001:end,2);
Z = V(200001:end,3);

columnize = @(M) M(:);

Jac = @(F)       [     0,   -1,      -1  ;
1,  0.2,       0  ;
F(3),    0,  F(1)-9 ];

msf_func = @(t,F,ind) [Jac([X(ind),Y(ind),Z(ind)])*[F(1);F(2);F(3)];
columnize(Jac([X(ind),Y(ind),Z(ind)])*[ F(4), F(7), F(10);
F(5), F(8), F(11);
F(6), F(9), F(12)]) ];
initi = [init,1,0,0,0,1,0,0,0,1];
[~,LE] = ode45(@(t,F) msf_func(t,F,floor(t/0.001+1)),0:0.001:2000,initi);
``````

After some analysis in a Slack thread, we see that this system is chaotic and so the first solve is not going to translate well between systems. This seems to effect the second solve because the stability of the second solve is then dependent on where in the attractor the solution is, which is effectively random between systems because of the chaos. But in general, it doesnâ€™t seem like this second system is guaranteed to be stable since itâ€™s a forward sensitivity calculation, which is known to be unstable on chaotic systems. For more information on calculating such sensitivities see Shadowing Methods for Forward and Adjoint Sensitivity Analysis of Chaotic Systems | FS where you normally need to use some kinds of renormalization to account for the instability of the tangent space.