I have a system with 5 linear ODEs. The code runs and the results are just what I expect. In a second system, I expand the first system by duplicating the five equations. There are now 10 equations where there is no connection between the first and the second five equations.
I expected that the the second system gives me identical results to the first system, but duplicated. Instead, I get the following error:
┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://diffeq.sciml.ai/stable/solvers/ode_solve/#Stiff-Problems).
└ @ SciMLBase C:\Users\arvka\.julia\packages\SciMLBase\OfgnG\src\integrator_interface.jl:491
retcode: MaxIters
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 999977-element Vector{Float64}:
I tried larger maxiters (1e7) but it didn’t help.
Here is the first system:
using Pkg
Pkg.activate(".")
using LinearAlgebra
using DifferentialEquations
u0 = [1800.0 1800.0 1800.0 1800.0 1800.0]
p = Any[[5.846581865622962, 0.998611, 0.971057, 0.683772, 0.00629], [0.001, 0.001, 0.001, 0.001, 0.000322], 64000]
function model1(du, u, p, t)
nsites = 1
r, d, K = p
counter = 0
for site in 1:nsites
# Stages 1
counter += 1
stage = 1
prev_stage = 5
du[counter] = (r[stage] * u[site, prev_stage]) -
(r[stage+1] * u[site, stage]) -
(d[stage] * u[site, stage])
# Stage 2
counter += 1
stage = 2
prev_stage = 1
du[counter] = (r[stage] * u[site, prev_stage]) -
(r[stage+1] * u[site, stage]) -
(d[stage] * u[site, stage])
#stage 3
counter += 1
stage = 3
prev_stage = 2
du[counter] = (r[stage] * u[site, prev_stage] * ((K - u[site, stage]) / K)) -
(r[stage+1] * u[site, stage]) -
(d[stage] * u[site, stage])
# stage 4
counter += 1
stage = 4
prev_stage = 3
du[counter] = (r[stage] * u[site, prev_stage] * ((K - u[site, stage]) / K)) -
(r[stage+1] * u[site, stage]) -
(d[stage] * u[site, stage])
# stage 5
counter += 1
stage = 5
prev_stage = 4
du[counter] = (r[stage] * u[site, prev_stage] * ((K - u[site, stage]) / K)) -
(0.0 * u[site, stage]) -
(d[stage] * u[site, stage])
end
end
tspan = (1.0, 200.0)
prob = ODEProblem(model1, u0, tspan, p)
sol = solve(prob)
And here is the second system. Note that the only difference in this case is nsites=2, instead of 1 and u0 is a matrix with two rows.
u02 = [1800.0 1800.0 1800.0 1800.0 1800.0; 1800.0 1800.0 1800.0 1800.0 1800.0]
p2 = [[5.846581865622962, 0.998611, 0.971057, 0.683772, 0.00629], [0.001, 0.001, 0.001, 0.001, 0.000322], 64000.0]
function model2(du, u, p, t)
nsites = 2
r, d, K = p
counter = 0
for site in 1:nsites
# Stages 1
counter += 1
stage = 1
prev_stage = 5
du[counter] = (r[stage] * u[site, prev_stage]) -
(r[stage+1] * u[site, stage]) -
(d[stage] * u[site, stage])
# Stage 2
counter += 1
stage = 2
prev_stage = 1
du[counter] = (r[stage] * u[site, prev_stage]) -
(r[stage+1] * u[site, stage]) -
(d[stage] * u[site, stage])
#stage 3
counter += 1
stage = 3
prev_stage = 2
du[counter] = (r[stage] * u[site, prev_stage] * ((K - u[site, stage]) / K)) -
(r[stage+1] * u[site, stage]) -
(d[stage] * u[site, stage])
# stage 4
counter += 1
stage = 4
prev_stage = 3
du[counter] = (r[stage] * u[site, prev_stage] * ((K - u[site, stage]) / K)) -
(r[stage+1] * u[site, stage]) -
(d[stage] * u[site, stage])
# stage 5
counter += 1
stage = 5
prev_stage = 4
du[counter] = (r[stage] * u[site, prev_stage] * ((K - u[site, stage]) / K)) -
(0.0 * u[site, stage]) -
(d[stage] * u[site, stage])
end
end
tspan = (1.0, 200.0)
prob2 = ODEProblem(model2, u02, tspan, p2)
sol = solve(prob2)