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