I am using DifferentialEquations.jl to solve these ordinary differential equations:
where dN/dt is bacteria susceptible, dI/dt is bacteria infected and dV/dt is phages. I wrote the following code:
function dynamo!(du, u, p, t) μ, κ, φ, ω, η, β = p #= du = susceptible du = infected du = phages =# du = ((μ * u) * (1 - ((u+u)/κ))) - (φ * u * u) - (ω * u) du = (φ * u * u) - (η * u) - (ω * u) du = (β * η * u) - (φ * u * u) - (ω * u) end # set parameters mu = 0.16 # maximum growth rate kappa = 2.2*10^7 # maximum population density phi = 10.0^-9 # adsorption rate beta = 50.0 # burst size omega = 0.05 # outflow eta = 1.0 # lyse rate 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 u0 = [s0, i0, v0] p = [mu, kappa, phi, omega, eta, beta] tspan = (0.0, tmax) # set solver prob = ODEProblem(dynamo!, u0, tspan, p) soln = solve(prob)
and it works, but the plot I get is this:
instead of this:
I guess the difference is in the algorithm chosen by the solver (unless there are flaws in the conversion of the equations).
Can I tune the solver in order to get the expected results? and if I change the algorithm, won’t I introduce a bias in the results, for they will be what they should be but what they are supposed to be? In this case, I have a certain solution available, but what for real life data?