I am trying to transfer a Matlab model to Julia, but am running into instability issues. I have found this discussion:
[Handling Instability When Solving ODE Problems]
and tried some of the suggested actions, but am struggling to diagnose my specific problem. I think that I just lack understanding – I am new to Julia and to working with differential equations – and hope that someone can point me in the right direction.
I have the following (simplified slightly):
using DifferentialEquations, LinearAlgebra, Random, Plots, BenchmarkTools
function f_h!(dz, z, p, t)
xn, PfH, PfHC = p
# Get xn at time t
dt = 0.05
idx = max(Int(ceil(t/dt)), 1)
# Exponentiate state variables
z[:,2:4] = exp.(z[:,2:4])
# Compute fv and ff
fv = z[:,3].^(1 ./ H[:,3])
ff = (1 .- (1 .- H[:,5])) .^ (1 ./ z[:,2]) ./ H[:,5]
# Compute state equations
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
PfH = [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]
PfHC = [0.0 -0.693147; 0.0 0.0]
z0 = 0.001*ones(2, 6) # initial condition
tspan = (0.0, 600.0) # time span
p = [xn, PfH, PfHC] # parameters
ts = range(0.05, step=0.05, stop=600.0) # sample times
ode_h = ODEProblem(f_h!, z0, tspan, p)
sol = solve(ode_h, Tsit5(), saveat=ts)
plot(sol)
which fails on the first time step with:
┌ Warning: Instability detected. Aborting
└ @ SciMLBase /home/birtet/.julia/packages/SciMLBase/n3U0M/src/integrator_interface.jl:351
The input argument xn is the solution of another ODE Problem (sampled at 0.05 s time steps), which does not suffer from instability issues:
using DifferentialEquations, LinearAlgebra, Random, Plots, BenchmarkTools
function f_n!(x, p, t)
A, B, C, u = p
# Get input at time t
dt = 0.05
idx = max(Int(ceil(t/dt)), 1)
ut = u[idx,:]
# Modulations by inputs
J = A
for i in 1:size(B, 3)
J += B[:,:,i]*ut[i]
end
J*x + C*ut
end
A = [-1.1 0.0; 0.0 -1.1]
B = zeros(2, 2, 2)
C = [0.3 0.0; 0.0 0.3]
x0 = 0.001*ones(2, 1) # initial condition
tspan = (0.0, 600.0) # time span
p = [A, B, C, u] # parameters
ts = range(0.05, step=0.05, stop=600.0) # sample times
ode_n = ODEProblem(f_n!, x0, tspan, p)
sol = solve(ode_n, Tsit5(), saveat=ts)
xn = sol.u
plot(sol);
(I know that in this particular example the B-matrix seems unnecessary, but its elements aren’t always set to 0.) In f_n! the parameter u represents exogenous/experimental inputs. For reasons of brevity I have not included the code to generate u, but it is a 12000×2 Matrix{Int64} of 0s and 1s.
I have confirmed that the derivative functions in f_h! are exactly the same in MATLAB and Julia. (As a side note, in the MATLAB version f_n! and f_h! are combined into a single ODE problem, but I don’t think that should make a difference - or at least I got the same instability warning when I tried this in Julia.)
When I tried solving ode_h with CVODE_BF(), I got the following error:
[CVODE ERROR] CVode
At t = 0 and h = 4.95839e-12, the corrector convergence test failed repeatedly or with |h| = hmin.
When using radau(), there is no error or warning message, but the output is filled with NaNs. Based on the answers to the other post, I think this means that there is an issue with one (or more) of my derivative functions, but I don’t understand what it is. The denominators are all non-zero (at least initially). I’d appreciate any insights into what could be causing the problem.
(I also wanted to read up on defining a Jacobian, but this link no longer works: