How to solve instability in ODE problem?

New to Julia
, ,
1

I am trying to implement a model that is implemented in Matlab in Julia using the DifferentialEquations package. In Julia, the solver stops after a couple of iterations with a warning that an instability has been detected. Here is my (simplified) code:

function h(dz, z, p, t)
    
    H, HC, xn = p

    # Get value of xn at time t
    idx = max(Int(ceil(t/0.05)), 1)
    
    z[:,2:4] = exp.(z[:,2:4])

    fv = z[:,3].^(1 ./ H[:,4])

    ff = (1 .- (1 .- H[:,5])) .^ (1 ./ z[:,2]) ./ H[:,5]
    
    dz[:,1] = xn[idx] .- H[:,1].*z[:,1] - H[:,2].*(z[:,2] .- 1)
    dz[:,2] = z[:,1] ./ z[:,2] 
    dz[:,3] = (z[:,2] - fv + HC*z[:,5]) ./ (H[:,3].*z[:,3]) 
    dz[:,4] = (ff.*z[:,2] - fv.*z[:,4]./z[:,3] + HC*z[:,6])./(H[:,3].*z[:,4])
    dz[:,5] = (-z[:,5] + z[:,3] .- 1)./(H[:,7])
    dz[:,6] = (-z[:,6] + z[:,4] .- 1)./(H[:,7])
end

z0 = zeros(2, 6) # initial condition
tspan = (0.0, 600) # time span
H = [0.65  0.41  0.98  0.32  0.34  0.4584  0.5;
        0.65  0.41  0.98  0.32  0.34  0.4584  0.5]
HC = [0.0  0.5;
          0.0  0.0]
p = [H, HC, xn] # parameters
ts = range(0.05, step=0.05, stop=600) # sample times

ode = ODEProblem(h, z0, tspan, p)
sol = solve(ode, Tsit5(), saveat=ts, tstops=ts)

The input parameter xn is the solution of another ODE problem; for simplicity, assume a vector of 0s and 1s:

xn = zeros(12000,2)
xn[6000:6004] .= 1

The problem occurs with the exponentiation of z[:,2:4] inside the function. (If I remove this line from the function and exponentiate the corresponding elements of the initial states z0, the ODE can be solved. However, that’s not the model I’m trying to implement.) In the first 2 iterations the solutions z are exactly 0.0 and 1.0, then they start fluctuating slightly, until after a few more iterations (about 8) there is a deviation from 1 that is large enough that the values jump on the next iteration and then become infinitely large.

Does anyone have a suggestion what I could do to deal with this instability?

2

is this a stiff ODE? If so, you need a stiff solver. Try replacing Tsit5() with AutoTsit5(Rosenbrock23()) which will dynamically switch to a stiff solver if it thinks it needs to. For more details, see ODE Solvers · DifferentialEquations.jl

3

I’m pretty sure your ODE definition should not include something which (selectively?) exponentiates the state values. I think, definitionally, your ODE should be of the form x' = f(x,p,t). In fact, I am not surprised that, since every time h is evaluated your state is being (selectively) exponentiated, your state quickly diverges.

In this case, it may help to see the original Matlab code you are re-implementing, or even link to the mathematical model you are using, since I can’t tell from what is included why this exponentiation step is necessary. Hope that is helpful!

1 Like
4

Yeah, don’t do that. The only reason why your MATLAB code even works is because this does not do what you expect in MATLAB and it actually ignores this. We don’t ignore your changes, so it will cause bad things. See: