# Instability detected in solving ODE with DifferentialEquations.jl

**URL:** <https://discourse.julialang.org/t/instability-detected-in-solving-ode-with-differentialequations-jl/99792>\
**Category:** Numerics\
**Tags:** differentialequation\
**Created:** [June 2, 2023, 10:44pm UTC](https://discourse.julialang.org/t/instability-detected-in-solving-ode-with-differentialequations-jl/99792 "2023-06-02T22:44:49Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Rezair](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rezair/32/48155_2.png) [@Rezair](https://discourse.julialang.org/u/Rezair)\
**Post date:** [June 2, 2023, 10:44pm UTC](https://discourse.julialang.org/t/instability-detected-in-solving-ode-with-differentialequations-jl/99792/1 "2023-06-02T22:44:49Z")

</div>

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:

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 3, 2023, 12:40am UTC](https://discourse.julialang.org/t/instability-detected-in-solving-ode-with-differentialequations-jl/99792/2 "2023-06-03T00:40:50Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Rezair](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rezair/32/48155_2.png) [@Rezair](https://discourse.julialang.org/u/Rezair)\
**Post date:** [June 3, 2023, 11:15am UTC](https://discourse.julialang.org/t/instability-detected-in-solving-ode-with-differentialequations-jl/99792/3 "2023-06-03T11:15:20Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 3, 2023, 3:36pm UTC](https://discourse.julialang.org/t/instability-detected-in-solving-ode-with-differentialequations-jl/99792/4 "2023-06-03T15:36:43Z")

</div>

> [@Rezair](#):
>
> ```julia
> ind = floor(Int,t/p.dT)+1
> x, y, z = p.steadystate[1,ind], p.steadystate[2,ind], p.steadystate[3,ind]
> 
> ```

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:

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

```

and then:

```julia
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.

---

<div class="post-metadata">

**Author:** ![Rezair](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rezair/32/48155_2.png) [@Rezair](https://discourse.julialang.org/u/Rezair)\
**Post date:** [June 3, 2023, 5:53pm UTC](https://discourse.julialang.org/t/instability-detected-in-solving-ode-with-differentialequations-jl/99792/5 "2023-06-03T17:53:30Z")

</div>

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

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 3, 2023, 5:57pm UTC](https://discourse.julialang.org/t/instability-detected-in-solving-ode-with-differentialequations-jl/99792/6 "2023-06-03T17:57:24Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Rezair](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rezair/32/48155_2.png) [@Rezair](https://discourse.julialang.org/u/Rezair)\
**Post date:** [June 3, 2023, 6:03pm UTC](https://discourse.julialang.org/t/instability-detected-in-solving-ode-with-differentialequations-jl/99792/7 "2023-06-03T18:03:53Z")

</div>

Sure, here is the MATLAB code:

```julia
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);

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 9, 2023, 10:23pm UTC](https://discourse.julialang.org/t/instability-detected-in-solving-ode-with-differentialequations-jl/99792/8 "2023-06-09T22:23:22Z")

</div>

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](https://frankschae.github.io/post/shadowing/) where you normally need to use some kinds of renormalization to account for the instability of the tangent space.