Instability detected in solving ODE with DifferentialEquations.jl

Specific Domains Numerics
1

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
    steadystate::Matrix{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
    x, y, z = p.steadystate[1,ind], p.steadystate[2,ind], p.steadystate[3,ind]
    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)
2

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?

3

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.

4

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)
    x, y, z = p.steadystate(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
5

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
    steadystate::ODESolution
end

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

function (p::RosslerMSF)(du,u,par,t)
    x, y, z = p.steadystate(t+200)
    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)
6

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

7

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);
8

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.