DifferentialEquations: mysterious problem

New to Julia
1

I’m simulating a model with 8 states using DifferentialEquations and get the expected result:

image
image895×589 39.9 KB

[Red line is simulation; scatter plot is experimental values.]

Next, I double the model by simulating two identical models (16 states) which are individually independent. When I re-do the simulation, there is a deviation:

image
image904×599 46.4 KB

These results are from Jupyter notebook. I get the same results if I run the code in VC. I have tried to reduce reltol and abstol, without any effect.

Here is the ODE for the first plot:

function mod1!(dx,x,p,t)
    S,E,I,C,A,R,CC,xm = x
    k0_i,k_e,k_ca,eta,k_r,tm,Tm,N,u = p
    k_i = k0_i*xm
    k_c = eta*k_ca
    k_a = (1-eta)*k_ca
    #
    dS = -k_i*(I+A)*S/N
    dE = k_i*(I+A)*S/N - k_e*E
    dI = k_e*E - k_ca*I
    dC = k_c*I - k_r*C
    dA = k_a*I - k_r*A
    dR = k_r*(C+A)
    dCC = k_c*I
    dxm = (u(t+tm) - xm)/Tm
    dx .= [dS,dE,dI,dC,dA,dR,dCC,dxm]
end

In this model, u is an ConstantInterpolation function. I call the function as follows:

x0 = [S0,E0,I0,C0,A0,R0,CC0,xm0]
#
p = [k0_i,k_e,k_ca,eta,k_r,tm,Tm,N,u]
#
prob = ODEProblem(mod1!,x0,tspan,p)
sol = solve(prob)

The model for the second plot is:

function mod2!(dx,x,p,t)
    S,E,I,C,A,R,CC,xm,S2,E2,I2,C2,A2,R2,CC2,xm2 = x
    k0_i,k0_i2,k_e,k_ca,eta,k_r,tm,Tm,tm2,Tm2,N,u,N2,u2 = p
    k_i = k0_i*xm
    k_i2 = k0_i2*xm2
    k_c = eta*k_ca
    k_a = (1-eta)*k_ca
    #
    dS = -k_i*(I+A)*S/N
    dE = k_i*(I+A)*S/N - k_e*E
    dI = k_e*E - k_ca*I
    dC = k_c*I - k_r*C
    dA = k_a*I - k_r*A
    dR = k_r*(C+A)
    dCC = k_c*I
    dxm = (u(t+tm) - xm)/Tm
    #
    dS2 = -k_i2*(I2+A2)*S2/N2
    dE2 = k_i2*(I2+A2)*S2/N2 - k_e*E2
    dI2 = k_e*E2 - k_ca*I2
    dC2 = k_c*I2 - k_r*C2
    dA2 = k_a*I2 - k_r*A2
    dR2 = k_r*(C2+A2)
    dCC2 = k_c*I2
    dxm2 = (u2(t+tm2) - xm2)/Tm2
    #
    dx .= [dS,dE,dI,dC,dA,dR,dCC,dxm,dS2,dE2,dI2,dC2,dA2,dR2,dCC2,dxm2]
end

I call the function as follows:

x0 = [S0,E0,I0,C0,A0,R0,CC0,xm0,S0,E0,I0,C0,A0,R0,CC0,xm0]
#
p = [k0_i,1.5*k0_i,k_e,k_ca,eta,k_r,tm,Tm,tm,Tm,N,u,N,u]
#
tin = 0.
tfin = 245
tspan = (tin,tfin)
prob = ODEProblem(mod2!,x0,tspan,p)
sol = solve(prob)

Observe that in mod2!, the two models (S… xm, and S2… xm2) are (i) identical, and (ii) independent of each other. The only difference when simulating them is that for the second part of mod2!, I change parameter k0_i by multiplying it with 1.5.

If I set k0_i2 = k0_i, then I get the correct result for the simulation, but since the models are independent, I should get the same result no matter what k0_i2 is??

Question: What am I doing wrong??
[I can provide a full set of code, if necessary, but it is rather lengthy.]

2

Doubling up would increase the error estimate because of the norm and thus would decrease the dt and cause the solution to be more accurate. If you use a lower solver tolerance the solution should be indistinguishable.

1 Like
3

Ah – thanks for suggestion. I need to insert reltol = 1e-6 for mod1!! That causes that model to deviate from my carefully tuned input function u.

To get better fit between the results of mod1 and mod2, I also need to use some similar reltol for mod2!…
[I thought the problem was with `mod2!`, so I only tested changing `reltol` for that model, but the problem was with both models…]