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

Specific Domains Numerics
1

Hello,
I am simulating a microbial consortium (2 bacteria, 2 phages) but I get a MethodError during the solving step. Since I have duplicated items (two growth rates, two infection rates etc) I used vectors of 2 elements for each parameter. I think there is some problem in extracting numbers from these vectors. The implementation is:

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]            # initial infected population
r0      = [0, 0]            # initial resistant population
u0 = [x0, i0, r0, v0]  # S, I, R, V
parms = (μ=mu, κ=kappa, ω=omega, δ=delta, η=eta, β=beta, λ=lambda)
prob = ODEProblem(FUNK!, u0, (0.0, 1e4), parms)
soln = solve(prob, Rosenbrock23())

julia> soln = solve(prob, Rosenbrock23())
ERROR: MethodError: no method matching real(::String)
Closest candidates are:
  real(::ChainRulesCore.AbstractThunk) at ~/.julia/packages/ChainRulesCore/GUvJT/src/tangent_types/thunks.jl:35
  real(::Complex) at ~/src/julia/share/julia/base/complex.jl:72
  real(::LinearAlgebra.Symmetric{var"#s859", S} where {var"#s859"<:Real, S<:(AbstractMatrix{<:var"#s859"})}) at ~/src/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/symmetric.jl:353
  ...

How can I solve the problem?
Thank you

PS:
what is the difference between [1, 2] and (1, 2)?

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
3

Thank you! I am getting problems with u0 now. Only a subset of soln gives results, the other is simply a vector of 0s. Even if I use u0 = [1e6, 0, 0, 1e4, 1e7, 0, 0, 1e3] I get the same problem…

4

You have 8 variables but wrote ODEs for only 6.

1 Like
5

Yes, I had two values more. It worked fine:

CASE_0
CASE_01920×1440 148 KB