# DifferentialEquations: mysterious problem

**URL:** <https://discourse.julialang.org/t/differentialequations-mysterious-problem/49836>\
**Category:** New to Julia\
**Created:** [November 9, 2020, 12:43pm UTC](https://discourse.julialang.org/t/differentialequations-mysterious-problem/49836 "2020-11-09T12:43:22Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [November 9, 2020, 12:43pm UTC](https://discourse.julialang.org/t/differentialequations-mysterious-problem/49836/1 "2020-11-09T12:43:22Z")

</div>

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

 ![image](https://global.discourse-cdn.com/julialang/original/3X/6/2/62e8d06f731ccdc1eccbf0c43d0f1acd869e5a1a.png)  
[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](https://global.discourse-cdn.com/julialang/original/3X/a/0/a0a9256369fd0e2266632c02a5e723d18abb0646.png)

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:

```julia
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:

```julia
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:

```julia
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:

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 9, 2020, 1:09pm UTC](https://discourse.julialang.org/t/differentialequations-mysterious-problem/49836/2 "2020-11-09T13:09:23Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [November 9, 2020, 1:24pm UTC](https://discourse.julialang.org/t/differentialequations-mysterious-problem/49836/3 "2020-11-09T13:24:45Z")

</div>

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