# Adjust solution obtained from DifferentialEquations.jl

**URL:** <https://discourse.julialang.org/t/adjust-solution-obtained-from-differentialequations-jl/31452>\
**Category:** General Usage\
**Created:** [November 24, 2019, 3:34pm UTC](https://discourse.julialang.org/t/adjust-solution-obtained-from-differentialequations-jl/31452 "2019-11-24T15:34:34Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [November 24, 2019, 3:34pm UTC](https://discourse.julialang.org/t/adjust-solution-obtained-from-differentialequations-jl/31452/1 "2019-11-24T15:34:34Z")

</div>

Hello,  
I am using DifferentialEquations.jl to solve these ordinary differential equations:  
 ![20191108_064445_001](https://global.discourse-cdn.com/julialang/original/3X/b/6/b6b3761f9858b684133a337f6c4095991cc523d5.jpeg)  
where dN/dt is bacteria susceptible, dI/dt is bacteria infected and dV/dt is phages. I wrote the following code:

```julia
function dynamo!(du, u, p, t)
    μ, κ, φ, ω, η, β = p
    #=
    du[1] = susceptible
    du[2] = infected
    du[3] = phages
    =#
    du[1] = ((μ * u[1]) * (1 - ((u[1]+u[2])/κ))) - (φ * u[1] * u[3]) - (ω * u[1])
    du[2] = (φ * u[1] * u[3]) - (η * u[2]) - (ω * u[2])
    du[3] = (β * η * u[2]) - (φ * u[1] * u[3]) - (ω * u[3])
end

# set parameters
mu = 0.16 # maximum growth rate
kappa = 2.2*10^7 # maximum population density
phi = 10.0^-9 # adsorption rate
beta = 50.0 # burst size
omega = 0.05 # outflow
eta = 1.0 # lyse rate
tmax = 4000.0 # time span 0-tmax
s0 = 50000.0 # initial susceptible population
i0 = 0.0 # initial infected population
v0 = 80.0 # initial phage population

u0 = [s0, i0, v0]
p = [mu, kappa, phi, omega, eta, beta]
tspan = (0.0, tmax)

# set solver
prob = ODEProblem(dynamo!, u0, tspan, p)
soln = solve(prob)

```

and it works, but the plot I get is this:  
 ![p25](https://global.discourse-cdn.com/julialang/original/3X/6/4/644f0bb1fec590255ba64d51050c7dfa77bcaf5b.png)  
instead of this:

 ![20191124_162236](https://global.discourse-cdn.com/julialang/original/3X/a/f/af617e0f87c79ff1871f0d19df5ad45b1e382788.jpeg)  
I guess the difference is in the algorithm chosen by the solver (unless there are flaws in the conversion of the equations).  
Can I tune the solver in order to get the expected results? and if I change the algorithm, won’t I introduce a bias in the results, for they will be what they should be but what they are supposed to be? In this case, I have a certain solution available, but what for real life data?  
Thank you

---

<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 24, 2019, 3:45pm UTC](https://discourse.julialang.org/t/adjust-solution-obtained-from-differentialequations-jl/31452/2 "2019-11-24T15:45:26Z")

</div>

It is a little confusing that you use different symbols in your code compared to the model picture you include.

Anyway, observe that your plot (I don’t know which package you used) uses linear scaling of the ordinate axis (“y axis”), while the plot you compare with uses logarithmic scaling of the ordinate axis.

If you use the `Plots` package, you can change to logarithmic scaling by adding keyword `yscale = :log10` in the plot function (if my memory serves me well).

---

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [November 24, 2019, 6:46pm UTC](https://discourse.julialang.org/t/adjust-solution-obtained-from-differentialequations-jl/31452/3 "2019-11-24T18:46:10Z")

</div>

I though adjustment of the scale was automatic, sorry. But you were right: it was a problem of scale. However, with `yscale = :log10` I got  
 ![dynSim](https://global.discourse-cdn.com/julialang/original/3X/3/e/3e363e287503b4b991e4b4cd892b6226e5a64270.png)  
With `ax = gca(); ax.set_yscale("log")` I got:  
 ![2](https://global.discourse-cdn.com/julialang/original/3X/9/4/9413f49327aabc1585bae735626d187d89f89766.png)  
Which is consistent with the reference figure. I don’t know why the colours of the lines are different from those of the legend though. Thank you

---

<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 24, 2019, 8:28pm UTC](https://discourse.julialang.org/t/adjust-solution-obtained-from-differentialequations-jl/31452/4 "2019-11-24T20:28:31Z")

</div>

I’ve modified your code slightly… First, I use the following packages:

```julia
using Plots, LaTeXStrings; pyplot()
using DifferentialEquations

```

Next, I leave your function as it is, but change 3 of the following values:

```julia
# set parameters
mu = 0.16 # maximum growth rate
kappa = 2.2e7 # maximum population density
phi = 1e-9 # adsorption rate
beta = 50.0 # burst size
omega = 0.05 # outflow
eta = 1.0 # lyse rate
tmax = 4000.0 # time span 0-tmax
s0 = 50000.0 # initial susceptible population
i0 = 1e-6 # initial infected population
v0 = 80.0 # initial phage population

```

I’ve changed the following:

1. `kappa = 2.2*10^7` → `kappa = 2.2e7`. The latter is the proper way to write 2.2\cdot 10^7. Your method is potentially dangerous for the following reason: `10^7` is treated as an _integer_ in Julia. When you multiply _floating point_ value `2.2` by integer `10^7` you _do_ get a floating point number. So, in this case: no problem. But consider that `10^18` becomes `1000000000000000000` =10^{18}, while `10^19` becomes `-8446744073709551616`\neq 10^{19}… what happens here is that 10^{19} is larger than the largest integer that Julia normally handles, and what you see is integer overflow.
2. In general, `xen` where `x` is an _integer_ or \* floating point\* number, `n` is a positive or negative integer, and `e` is … `e` – creates a floating point number x\cdot 10^n. `xen` is the proper Julia way to write x\cdot 10^n, while `x*10^n` is potentially dangerous.
3. I have thus also changed `phi = 10.0^-7` → `phi = 1e-7`.
4. I have also changed `i0 = 0` → `i0 = 1e-6`. The reason is that with `i0 = 0`, you get an error message with logarithmic scaling because the logarithm of 0 obviously is -\infty.

With the `Plots` package and `LaTeXStrings`, plotting can be done as follows:

```julia
plot(soln,label=["Susceptible bacteria" "Infected bacteria" "Bacteriophage"])
plot!(yscale=:log10,ylim=(1e2,1e9))
plot!(xlabel="time [h]", ylabel="concentration [1/mL]",title=L"\eta = 1.0")

```

leading to:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/7/9/7980e9e9887c961a0909de8cba118b0385f57573.png)  
Personally, I’d have scaled the time axis to make the plot easier to “read”. I don’t know how to do that with the data structure of `soln` (Maybe Chris R has a trick up his sleeve). One way to do it is to take out the solutions into a _matrix_:

```julia
Soln = reduce(hcat,soln.u)';

```

and plot the solution as follows:

```julia
plot(soln.t/24/7,Soln,lw=2,label=["Susceptible bacteria" "Infected bacteria" "Bacteriophage"])
plot!(yscale=:log10,ylim=(1e2,1e9))
plot!(xlabel="time [week]", ylabel="concentration [1/mL]",title=L"\eta = 1.0")

```

leading to:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/6/26556ba968ac3f3dab8ab83f74794e44a4cd34c5.jpeg)

The key advantage of the first way to plot the solution is that the solution `soln` includes an interpolation function which is used to generate a smooth plot, while in the second approach (`Soln`) `plot` uses linear interpolation between the data points. (Perhaps not so easy to see here…)

---

<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 25, 2019, 10:39pm UTC](https://discourse.julialang.org/t/adjust-solution-obtained-from-differentialequations-jl/31452/5 "2019-11-25T22:39:30Z")

</div>

The trick up the sleeve… how to scale the time axis (or other things) while keeping the DifferentialEquations interpolator for the solution…

So – keep the `soln`, i.e., solution as described above, and forget about the `Soln` matrix unpacking.

First, I pick out which variables in the solution I want to plot:

```julia
plot(soln,vars=[(0,1),(0,2),(0,3)],label=["Susceptible bacteria" "Infected bacteria" "Bacteriophage"])
plot!(yscale=:log10,ylim=(1e2,1e9))
plot!(xlabel="time [h]", ylabel="concentration [1/mL]",title=L"\eta = 1.0")

```

Here, this gives the same plot as above (`plot(soln, label=...)`). But I have specified that the first plot involves the mapping s\_0 \rightarrow s\_1 where s is my shorthand notation for `soln` and s\_0 = time (independent variable) while s\_1 is dependent variable 1; indicated by tuple `(0,1)`. Next, `(0,2)` is s\_0 \rightarrow s\_2, while `(0,3)` is s\_0 \rightarrow s\_3. Thus, `vars=[(0,1), (0,2), (0,3)]` implies to plot each variable (1,2,3) as a function of time (0).

But I can also include a mapping of the plots. Suppose I create a function for scaling hour to week. Thus:

```julia
t_h2w(t,u) = (t/24/7,u)

```

takes a pair of time and solution (`(t,u)`) and returns the same pair with the first argument (time) scaled from hour to week.

I then inform the `plot` command about this transformation as follows:

```julia
plot(soln,vars=[(t_h2w,0,1),(t_h2w,0,2),(t_h2w,0,3)],label=["Susceptible bacteria" "Infected bacteria" "Bacteriophage"])
plot!(yscale=:log10,ylim=(1e2,1e9))
plot!(xlabel="time [w]", ylabel="concentration [1/mL]",title=L"\eta = 1.0")

```

With this strategy, I get a weird result…

 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/7/2737ae31dbbe8a9197098e298e8e5eadae0d8cff.png)  
So… what happened is that plots changed the time axis to weeks (or rather: scaled t \rightarrow t/24/7, but did not change the time axis limit. To fix this, I simply need to do as follows:

```julia
 plot(soln,vars=[(t_h2w,0,1),(t_h2w,0,2),(t_h2w,0,3)],label=["Susceptible bacteria" "Infected bacteria" "Bacteriophage"])
plot!(yscale=:log10,ylim=(1e2,1e9),xlim=(0,tmax/24/7))
plot!(xlabel="time [w]", ylabel="concentration [1/mL]",title=L"\eta = 1.0")

```

and presto!

 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/a/0a6eb32075e2fb250e6210c305ccc7abe5e6d4fd.png)  
The advantage of this approach over my previous suggestion of extracting data into matrix `Soln` is that we have used the interpolating functions of `soln`, thus getting a smoother solution.