Day_04: Solving differential equations

Foreward note:

After ignoring the very obvious advice to install Julia from the official website or repo, I finally uninstalled the Arch package (via pacman -S julia) and the Anaconda package (via conda install -c conda-forge julia) and installed Julia v1.7.2 through

• git clone https://github.com/JuliaLang/julia.git

• cd julia

• make

• ./julia I thought my Arch system would be up-to-date enough for the Julia packages or maybe Jupyter would work better if I used the Anaconda package. All wrong, just follow the well-documented path and install the correct packages. Learn from my mistakes.

Today, I wanted to solve some differential equations. The DifferentalEquations documentation has some great initial examples. I chose the 1D exponential model to start,

$\frac{du}{dt} = \lambda u$

where $\lambda$ is the exponential growth/decay parameter. For positive values, the more you have, the more you get. For negative values, the more you have, the more you lose.

using Plots, LaTeXStrings
using DifferentialEquations

[ Info: Precompiling Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80] (cache misses: wrong dep version loaded (2))

Failed to precompile Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80] to "/home/ryan/.julia/compiled/v1.11/Plots/jl_bTPiSc".

Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IO, internal_stdout::IO, keep_loaded_modules::Bool; flags::Cmd, cacheflags::Base.CacheFlags, reasons::Dict{String, Int64}, isext::Bool)
    @ Base ./loading.jl:3085
  [3] (::Base.var"#1082#1083"{Base.PkgId})()
    @ Base ./loading.jl:2492
  [4] mkpidlock(f::Base.var"#1082#1083"{Base.PkgId}, at::String, pid::Int32; kwopts::@Kwargs{stale_age::Int64, wait::Bool})
    @ FileWatching.Pidfile ~/projects/julia/usr/share/julia/stdlib/v1.11/FileWatching/src/pidfile.jl:95
  [5] #mkpidlock#6
    @ ~/projects/julia/usr/share/julia/stdlib/v1.11/FileWatching/src/pidfile.jl:90 [inlined]
  [6] trymkpidlock(::Function, ::Vararg{Any}; kwargs::@Kwargs{stale_age::Int64})
    @ FileWatching.Pidfile ~/projects/julia/usr/share/julia/stdlib/v1.11/FileWatching/src/pidfile.jl:116
  [7] #invokelatest#2
    @ ./essentials.jl:1057 [inlined]
  [8] invokelatest
    @ ./essentials.jl:1052 [inlined]
  [9] maybe_cachefile_lock(f::Base.var"#1082#1083"{Base.PkgId}, pkg::Base.PkgId, srcpath::String; stale_age::Int64)
    @ Base ./loading.jl:3609
 [10] maybe_cachefile_lock
    @ ./loading.jl:3606 [inlined]
 [11] _require(pkg::Base.PkgId, env::String)
    @ Base ./loading.jl:2488
 [12] __require_prelocked(uuidkey::Base.PkgId, env::String)
    @ Base ./loading.jl:2315
 [13] #invoke_in_world#3
    @ ./essentials.jl:1089 [inlined]
 [14] invoke_in_world
    @ ./essentials.jl:1086 [inlined]
 [15] _require_prelocked(uuidkey::Base.PkgId, env::String)
    @ Base ./loading.jl:2302
 [16] macro expansion
    @ ./loading.jl:2241 [inlined]
 [17] macro expansion
    @ ./lock.jl:273 [inlined]
 [18] __require(into::Module, mod::Symbol)
    @ Base ./loading.jl:2198
 [19] #invoke_in_world#3
    @ ./essentials.jl:1089 [inlined]
 [20] invoke_in_world
    @ ./essentials.jl:1086 [inlined]
 [21] require(into::Module, mod::Symbol)
    @ Base ./loading.jl:2191

Setting up the differential equation

In this step, I built the differential equation as a function of u, p, and t. Here the

• u is the state of the differential equation, starting at an initial condition, u0 = 1/2.

• params are the parameters in the function, in this case its my value of $\lambda$

• t is the current time

f(u,params,t) = params*u;
u0 = 1/2;
tspan = (0.0,1.0);
p = 1.01;
prob = ODEProblem(f,u0,tspan,p)

ODEProblem with uType Float64 and tType Float64. In-place: false
timespan: (0.0, 1.0)
u0: 0.5

Solving the differential equation

Now, prob is an ODEProblem that contains

• the function, f(u, params, t)

• the initial value (at t = 0), u0

• the timespan, (0, 1.0)

• and the parameters, params = p

Solving this oridinary differential equation is now done with solve as such,

sol = solve(prob)

retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 5-element Vector{Float64}:
 0.0
 0.09964258706516003
 0.3457024247583422
 0.6776921908052249
 1.0
u: 5-element Vector{Float64}:
 0.5
 0.552938681151017
 0.7089376245893466
 0.9913594502399236
 1.3728004409033037

Checking the output

The output from solve is really nice. You can

• check the time steps, t: 5-element Vector {Float64}

• check the solutions, u: 5-element Vector {Float64}

• interpolate between values, sol(0.5)

• plus more

For example, here I plot the time vs current value. Then, add the interpolated point at t = 0.5

Plots.theme(:cooper)
halftime = 0.5
plot(sol.t, sol.u, markershape = :utriangle, label = "solutions")
scatter!([halftime], [sol(halftime)], markershape = :square, label = "interpolated")
plot!(xlabel = "years", ylabel = "current value")

Changing parameters in model

Next, I wanted to change the value of $\lambda$. Below I set up a for-loop to change $\lambda = 1.1,~ 1.01,~ 0,~ -1,~ -2$. I have to create a new problem in each case before running the solver and plotting the results. I also create smoother lines by interpolating with 11 time steps.

Note: when I plot inside a for-loop, I needed to create an initial empty plot and assign it to a variable, plts. Then, I added plots inside the loop using plot!. Finally, displaying the result with a call to plts

t = range(0, 1, 11)

plts = plot(bg=:white, xlabel = "time (years)", ylabel = "current value")
for p in [1.1, 1.01, 0, -1, -2]
    prob = ODEProblem(f,u0,tspan,p)
    sol = solve(prob)
    t = range(0, 1, 11)
    plot!(t, sol(t), label = string(L"$\lambda$","= $p"))
end
plts

Animating the effect of change

Its amazing what the effect of exponential growth can have on your health, knowledge, bank account, etc. If each day you’re increasing your knowledge by 1%, then by the end of the year you don’t have 100+365% knowledge, you have 3,778% of the knowledge you had a year ago. When you continuously integrate, it grow faster.

Here I build an animation from exponential decay, $\lambda = -2$ to exponential growth, $\lambda = 2$ every 0.1 steps.

First, I set up the plot background, labels, and ylimits. Then, use the @gif to loop through each value of $\lambda$ for parameter, p.

plts = plot(bg=:white, 
    xlabel = "time (years)", 
    ylabel = "current value", 
    ylims = (0, 4))

@gif for p in -2:0.1:2
    prob = ODEProblem(f,u0,tspan,p)
    sol = solve(prob)
    t = 0:0.1:1
    plot!(t, sol(t), label = "")#string(L"$\lambda$","= $p"))
end

┌ Info: Saved animation to 
│   fn = /home/ryan/Documents/Career_docs/cooperrc-gh-pages/Julia-learning/tmp.gif
└ @ Plots /home/ryan/.julia/packages/Plots/D9pfj/src/animation.jl:114

Wrapping up

Today, I had planned to look at a more advanced set of differential equations, but I ran into Pkg version issues because I did not follow the standard practices that were suggested. I love the “free” interpolation in the DifferentialEquations solver. I have made so many numerical solutions that I solve with unnecessary precision and number of timesteps because I wanted to compare specific values at given times. These first-order growth/decay models converge quickly, but sometimes I need to know a value at $t=0.6$, so I resolve the differential equation with smaller steps. It makes much more sense to interpolate.