Unhelpful interpolation for solution from differential equation

Specific Domains Modelling & Simulations
, ,
1

I am trying to use the solution from my differential equation to find local maxima. I am following the example found here, but am running into an issue. I have found that the interpolation of my solution is very poor. Here is a plot of the first 25 points of the numerical solution (scatter points) and the first 25 time units of the interpolated solution (line):

t_interp = 0:0.05:25
u_interp = map(i -> soln(i, idxs=3), t_interp)
t = soln.t[1:25]
u = map(i -> i[3], soln.u[1:25])

plot_1
plot_11200×900 36.7 KB

This interpolation is not helpful, and even a linear interpolation would be more accurate than the one shown here. Is there anything I can do to manipulate this interpretation to be more helpful?

Model:

function system_model(
    u::SVector{6, Float64},
    p_set::SVector{14, Float64},
    t::Float64 = 0.0,
)::SVector{6, Float64}

    Nᵣ, Nᵤ, Iᵣ, Iᵤ, P, Yᵣ = u
    Sᵣ, Sᵤ = Nᵣ - Iᵣ, Nᵤ - Iᵤ
    
    gᵣ, αᵤᵣ,
        hᵣ, μᵣ, χᵣ, uᵣ,
    gᵤ, αᵣᵤ,
        hᵤ, μᵤ, χᵤ, uᵤ,
        βᵤ, ζᵤ = p_set
    
    dNᵣdt = gᵣ * Nᵣ - (Nᵣ + αᵣᵤ * Nᵤ) * Nᵣ - μᵣ * Iᵣ
    dIᵣdt = Sᵣ * P - (Nᵣ + αᵣᵤ * Nᵤ) * Iᵣ - (hᵣ + μᵣ) * Iᵣ

    return SA[
        # dNᵣ/dt
        dNᵣdt,
        # dNᵤ/dt
        isinf(ζᵤ) ? 0. : (gᵤ - ζᵤ * t) * Nᵤ - (Nᵤ + αᵤᵣ * Nᵣ) * Nᵤ - μᵤ * Iᵤ,
        # dIᵣ/dt
        dIᵣdt,
        # dIᵤ/dt
        isinf(ζᵤ) ? 0. : βᵤ * Sᵤ * P - (Nᵤ + αᵤᵣ * Nᵣ) * Iᵤ - (hᵤ + μᵤ + ζᵤ * t) * Iᵤ,
        # dP/dt
        χᵣ * Iᵣ + χᵤ * Iᵤ - (uᵣ * Nᵣ + uᵤ * Nᵤ) * P - P,
        # dYᵣ/dt ("driven" state variable)
        (dIᵣdt - Yᵣ * dNᵣdt) / Nᵣ
    ]
end

Solver:

using StaticArrays: SA
using OrdinaryDiffEq: solve, ODEProblem, Tsit5, DiscreteCallback
using SciMLBase: ODESolution


_t_span = (0., 500.)
_secondary_host_tol = 1e-10

# makes secondary host go to zero when its population is low enough
_condition(u, _, _) = 0. < u[2] < _secondary_host_tol
_affect!(∫) = begin ∫.u = SA[∫.u[1], 0.0, ∫.u[3], 0.0, ∫.u[5], ∫.u[6]] end
_cb_extinction = DiscreteCallback(_condition, _affect!)

function _solve_system_ivp(
    u₀::AbstractVector{Float64},
    p_set::AbstractVector{Float64}
)::ODESolution
    prob = ODEProblem(system_model, SA[u₀...], _t_span, SA[p_set...])
    return solve(prob, Tsit5(), callback=_cb_extinction)
end
2

Whoops, by setting saveat=0.05 I see that this is a result of solver tolerances, and this interpolation is correct for what the solver found. Reducing absolute and relative tolerance resulted in much smoother output.

3

Could you please edit the code with your fix? Otherwise,
it doesn’t work much as a “Solution” to beginners or lurkers
still learning basics of DE by watching. Thanks! :slight_smile:

1 Like
4

change

return solve(prob, Tsit5(), callback=_cb_extinction)

to

return solve(prob, Tsit5(), callback=_cb_extinction, reltol=1e-6, abstol=1e-6)

or even lower/stricter tolerances such as 1e-9 or 1e-12.

1 Like