Julia DifferentialEquations: increase maxiter

Hello,
I have set an ODEProblem object, but the procedure is stopped (although I still get t a solution) because a higher maxiters is needed.
How can I check what is the default value and what is the syntax to increase the value?
Thank you

https://diffeq.sciml.ai/stable/basics/common_solver_opts/#Miscellaneous

EDIT: In my modest experience, running into maxiters may be a sign that something is wrong with your model. If it is stopping at some discontinuity you may need to do something structural like using callbacks to address a change in the model behavior at that point.

Yes. Or you’re using a non-stiff ODE solver on a stiff problem. Usually if you hit maxiters you’re doing something wrong. Though there are cases where you just truly need more iterations.

That makes sense. I have the following:

function doubleInfect!(du, u, p, t)
    μ, ν, κ, φ, ω, β, ρ, τ, η = p
    #=
    du[1] = susceptible
    du[2] = infected
    du[3] = resistant
    du[4] = phages
    =#
    TOT = u[1]+u[2]+u[3]
    du[1] = (μ * u[1]) * (1 - TOT/κ) - (φ * u[1] * u[3]) - (ω * u[1])
    du[2] = (φ * u[1] * u[3]) - (η * u[2]) - (ω * u[2])
    du[3] = ((ν * u[4]) * (1 - TOT/κ)) - (ω * u[4])
    du[4] = (β * η * u[2]) - (ρ * φ * u[1] * u[3]) - (ω * u[3])
end
mu    = 0.47        # maximum growth rate susceptible (B. longum)
nu    = 0.72        # maximum growth rate resistant (E. coli)
kappa = 2.2*10^7    # maximum population density
phi   = 10.0^-9     # adsorption rate
psi   = 10.0^-9     # adsorption rate alternative
omega = 0.05        # outflow
eta   = 1.0         # lyse rate
Eta   = 1.0         # lyse rate alternative
beta  = 50.0        # burst size
Beta  = 50.0        # burst size alternative
rho   = 0.6         # reinfection rate
tau   = 3.62        # latency time
tmax  = 4000.0      # time span 0-tmax
s0    = 50000.0     # initial susceptible population
i0    = 1.0e-9      # initial infected population
v0    = 80.0        # initial phage population
0    = 500.0       # initial resistant population
u0 = [s0, i0, r0, v0]
parms = [mu, nu, kappa, phi, omega, beta, rho, tau, eta]
tspan = (0.0, tmax)
prob = ODEProblem(doubleInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()))

I changed the algorithm from Tsit5() to AutoVern7(Rodas5()) and AutoTsit5(Rosenbrock23()), but I always get:

julia> prob = ODEProblem(doubleInfect!, u0, tspan, parms)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 4000.0)
u0: [50000.0, 1.0e-9, 500.0, 80.0]

julia> soln = solve(prob, AutoVern7(Rodas5()))
┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/3iigH/src/integrator_interface.jl:329
retcode: MaxIters
Interpolation: automatic order switching interpolation

What is the error here? How can I solve it? I get a solution, but I don’t like having errors dangling around…
Thank you

When I run it, it seems to blow up at 492.9 seconds (time in the simulation).

Are you sure you wrote down the equations you wanted? Your model doesn’t make too much sense.

du[3] = ((ν * u[4]) * (1 - TOT/κ)) - (ω * u[4])

simplifies with your values to du[3] = 0.68*u[4], so if the concentration of phages is always positive, then then the number of resistant diverges to infinity. As it diverges, TOT goes to infinity so (1 - TOT/κ) is going to go negative, so you’re going to then get very large negative terms. So your ODE won’t stay positive in u[1] (also the φ * u[1] * u[3] term). Once those two are negative, u[2] goes negative and…

I assume this isn’t what you wanted… but it’s what you told the computer to do. My guess is that one of those u[4] values is really u[3], so your equation should’ve been:

du[3] = ((ν * u[4]) * (1 - TOT/κ)) - (ω * u[3])

You are right, there we some problems with the indices. The model comes from the book Weitz’s Quantitative viral ecology, and the ODEs are:
1
which I converted in:


and then I added a bacterium that is not affected by phages:
3
hence:

function doubleInfect!(du, u, p, t)
    μ, ν, κ, φ, ω, β, τ, η = p
    #=
    du[1] = susceptible
    du[2] = infected
    du[3] = resistant
    du[4] = phages
    =#
    TOT = u[1]+u[2]+u[4]
    du[1] = (μ * u[1]) * (1 - TOT/κ) - (φ * u[1] * u[3]) - (ω * u[1])
    du[2] = (φ * u[1] * u[3]) - (η * u[2]) - (ω * u[2])
    du[3] = (β * η * u[2]) - (φ * u[1] * u[3]) - (ω * u[3])
    du[4] = (ν * u[4]) * (1 - TOT/κ) - 0                 - (ω * u[4])
    
end
# parameters
mu    = 0.47        # maximum growth rate susceptible (B. longum)
nu    = 0.72        # maximum growth rate resistant (E. coli)
kappa = 2.2*10^7    # maximum population density
phi   = 10.0^-9     # adsorption rate
psi   = 10.0^-9     # adsorption rate alternative
omega = 0.05        # outflow
eta   = 1.0         # lyse rate
Eta   = 1.0         # lyse rate alternative
beta  = 50.0        # burst size
Beta  = 50.0        # burst size alternative
rho   = 0.6         # reinfection rate
tau   = 3.62        # latency time
tmax  = 4000.0      # time span 0-tmax
s0    = 50000.0     # initial susceptible population
i0    = 0.0         # initial infected population
v0    = 80.0        # initial phage population
r0    = 0.0      # initial resistant population
u0 = [s0, i0, r0, v0]
parms = [mu, nu, kappa, phi, omega, beta, tau, eta]
tspan = (0.0, tmax)
# delay infection
condition(u, t, integrator) = t==1000           # time of inoculum
affect!(integrator) = integrator.u[3] += 500    # amount of inoculum
cb = DiscreteCallback(condition,affect!)
# implement model
prob = ODEProblem(doubleInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[1000])

This time:

julia> prob = ODEProblem(doubleInfect!, u0, tspan, parms)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 4000.0)
u0: [50000.0, 0.0, 0.0, 80.0]

julia> soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[1000])
retcode: Success
Interpolation: automatic order switching interpolation
t: 105-element Array{Float64,1}:

So the problem was indeed in the equations. Thank you