Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Large error when using Vern7 and Rodas5 #1060

Open
andreasnoack opened this issue Mar 11, 2020 · 1 comment
Open

Large error when using Vern7 and Rodas5 #1060

andreasnoack opened this issue Mar 11, 2020 · 1 comment

Comments

@andreasnoack
Copy link
Contributor

andreasnoack commented Mar 11, 2020

The problem below fails with a domain error because u[2] becomes negative. Since u[2] is determined by a linear subproblem I don't think it's expected. It's not just a bit negative but ~-155 for Vern7 and -1e37 for Rodas5. The solvers Tsit5, DP8, VCABM, and Rosenbrock23 all work fine.

using OrdinaryDiffEq

function hcv!(du, u, p, t)
    du[1] = -exp(p.logKa) + 180*isinfusion(t)
    du[2] = exp(p.logKa)*u[1] - exp(p.logKe)*u[2]
    du[3] = p.s - u[3]*(p.e*u[5] + p.d)
    du[4] = p.e*u[5]*u[3] - exp(p.logδ)*u[4]
    du[5] = p.p/((u[2]/exp(p.logVd)/exp(p.logEC50) + 1e-100)^exp(p.logn) + 1)*u[4] - exp(p.logc)*u[5]
end

rfx = (
    ηKa = 2.5507831516698394,
    ηKe = 0.42510262896450046,
    ηVd = 0.6897380500448431,
    ηn = -1.063040053371937,
    ηδ = -0.19115630736819755,
    ηc = -1.3073359078470523,
    ηEC50 = -0.7018572581731202
    )

param = (
    p        = 100.0,
    d        = 0.001,
    e        = 1e-7,
    s        = 20000.0,
    logKa   = log(0.80)  + rfx.ηKa,
    logKe   = log(0.15)  + rfx.ηKe,
    logVd   = log(100.0) + rfx.ηVd,
    logn    = log(2.0)   + rfx.ηn,
    logδ    = log(0.20)  + rfx.ηδ,
    logc    = log(7.0)   + rfx.ηc,
    logEC50 = log(0.12)  + rfx.ηEC50,
)

u0 = [
    0.0,
    0.0,
    exp(param.logc + param.logδ)/(param.p*param.e),
    (param.s*param.e*param.p - param.d*exp(param.logc + param.logδ))/(param.p*exp(param.logδ)*param.e),
    (param.s*param.e*param.p - param.d*exp(param.logc + param.logδ))/(exp(param.logc + param.logδ)*param.e)
    ]

tspan = (0.0, 40.0)

dosetimes = [8*i for i in 0:3]

function isinfusion(t)
    any(dosetime -> dosetime <= t <= dosetime + 1, dosetimes)
end

prob = ODEProblem(hcv!, u0, tspan, param)

_tstops = sort([dosetimes; dosetimes .+ 1])

sol = solve(prob, Vern7(), tstops=_tstops) # fails
sol = solve(prob, Rodas5(), tstops=_tstops) # fails
sol = solve(prob, Tsit5(), tstops=_tstops) # works
sol = solve(prob, DP8(), tstops=_tstops) # works
sol = solve(prob, VCABM(), tstops=_tstops) # works
@ChrisRackauckas
Copy link
Member

@YingboMa can you take a dive into this?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants