What is wrong with this system of ODEs?

Specific Domains
,
1

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)
2

I think the problem is how you access du or how you defined u02:

See this example

julia> A = [1 2 3; 4 5 6]
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

julia> A[2]
4 # instead of 2

This leads to problems when you acces du like this:

du[counter]

you could use this instead.

du[site,stage]

Here is the full code with nsites=2 that works correctly I think.

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[site,stage] = (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[site,stage] = (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[site,stage] = (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[site,stage] = (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[site,stage] = (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)
3

Oh I didn’t know that du has the same shape as u. This solved the problem. Thank you!