# DifferentialEquations: mysterious problem

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

[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:

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.]

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

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…]