How to solve `no method matching real(::String)` when solving ODEs?

Specific Domains Numerics
2

Arrays of arrays do not have vector semantics. Instead, use ComponentArrays.jl for this kind of nested model. Here’s an example:

using OrdinaryDiffEq, ComponentArrays

function FUNK!(du, u, p, t) 
    μ = p.μ   # bacterial growth rate
    κ = p.κ   # bacterial carrying capacity
    ω = p.ω   # system wash-out rate
    δ = p.δ   # phagial infection rate
    η = p.η   # phagial lysis rate (inverse latency)
    β = p.β   # phagial burst size
    λ = p.λ   # phagial decay rate
    #=
    du[1] = species A sensible
    du[2] = species A infected
    du[3] = phage a
    du[4] = species B sensible
    du[5] = species B infected
    du[6] = phage b
    =#
    N = u[1] + u[2] + u[4] + u[5]   # total bacteria
    ρ = 1 - N/κ                     # logistic term
    ϡ = (δ[1]*u[1]*u[3])            # upsampi: infected bacteria A
    ς = (δ[2]*u[4]*u[6])            # varsigma: infected bacteria B

    du[1] = (μ[1]*u[1]*ρ)       - ϡ                     - (ω*u[1])
    du[2] =                       ϡ     - (η[1]*u[2])   - (ω*u[2])
    du[3] = (η[1]*β[1]*u[2])    - ϡ     - (λ[1]*u[3])   - (ω*u[3])

    du[4] = (μ[2]*u[4]*ρ)       - ς                     - (ω*u[4])
    du[5] =                       ς     - (η[2]*u[5])   - (ω*u[5])
    du[6] = (η[2]*β[2]*u[4])    - ς     - (λ[2]*u[5])   - (ω*u[6])
end
kappa   = 2.2e7             # carrying capacity (from Weitz)
omega   = 0.05              # outflow (from Weitz)
mu      = [0.16, 0.31]      # growth rate susceptible (pathogen)
beta    = [50.0, 75.0]      # burst size           (from Weitz)
delta   = [1e-9, 1e-9]      # adsorption rate (from Weitz)
eta     = [1, 1]            # reciprocal of latency   (from Weitz)
lambda  = [0, 0]            # decay rate in hours     (from Weitz)
x0      = [1e6, 1e5]        # initial bacterial density
v0      = [1e5, 1e4]        # initial phage population
i0      = [0.0, 0]            # initial infected population
r0      = [0.0, 0]            # initial resistant population
u0 = ComponentArray(x=x0, i=i0, r=r0, v=v0)  # S, I, R, V
parms = ComponentArray(μ=mu, κ=kappa, ω=omega, δ=delta, η=eta, β=beta, λ=lambda)
prob = ODEProblem(FUNK!, u0, (0.0, 1e4), parms)
soln = solve(prob, Rosenbrock23())

There are no ODEs on the v` terms though, so I assume something is missing from the model, but this solves.

1 Like