{
"cells": [
{
"cell_type": "markdown",
"id": "29093ce0",
"metadata": {},
"source": [
"# Day_04: Solving differential equations\n",
"\n",
"> ## Foreward note:\n",
"> After ignoring the _very_ obvious advice to install Julia from the official website or repo, \n",
"> 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 \n",
"> - `git clone https://github.com/JuliaLang/julia.git`\n",
"> - `cd julia`\n",
"> - `make`\n",
"> - `./julia`\n",
"> 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_. \n",
"\n",
"Today, I wanted to solve some differential equations. The [`DifferentalEquations` documentation](https://diffeq.sciml.ai/stable/tutorials/ode_example/#Example-2:-Solving-Systems-of-Equations) has some great initial examples. I chose the 1D exponential model to start, \n",
"\n",
"$\\frac{du}{dt} = \\lambda u$\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "ce546225",
"metadata": {},
"outputs": [],
"source": [
"using Plots, LaTeXStrings\n",
"using DifferentialEquations"
]
},
{
"cell_type": "markdown",
"id": "a1f110a2",
"metadata": {},
"source": [
"## Setting up the differential equation\n",
"\n",
"In this step, I built the differential equation as a function of `u`, `p`, and `t`. Here the\n",
"- `u` is the state of the differential equation, starting at an initial condition, `u0 = 1/2`. \n",
"- `params` are the parameters in the function, in this case its my value of $\\lambda$\n",
"- `t` is the current time"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "461cc339",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\u001b[36mODEProblem\u001b[0m with uType \u001b[36mFloat64\u001b[0m and tType \u001b[36mFloat64\u001b[0m. In-place: \u001b[36mfalse\u001b[0m\n",
"timespan: (0.0, 1.0)\n",
"u0: 0.5"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(u,params,t) = params*u;\n",
"u0 = 1/2;\n",
"tspan = (0.0,1.0);\n",
"p = 1.01;\n",
"prob = ODEProblem(f,u0,tspan,p)"
]
},
{
"cell_type": "markdown",
"id": "de135797",
"metadata": {},
"source": [
"## Solving the differential equation\n",
"\n",
"Now, `prob` is an `ODEProblem` that contains\n",
"- the function, `f(u, params, t)`\n",
"- the initial value (at t = 0), `u0`\n",
"- the timespan, `(0, 1.0)`\n",
"- and the parameters, `params = p`\n",
"\n",
"Solving this oridinary differential equation is now done with `solve` as such,"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "1e758e05",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"retcode: Success\n",
"Interpolation: specialized 4th order \"free\" interpolation, specialized 2nd order \"free\" stiffness-aware interpolation\n",
"t: 5-element Vector{Float64}:\n",
" 0.0\n",
" 0.09964258706516003\n",
" 0.3457024247583422\n",
" 0.6776921908052249\n",
" 1.0\n",
"u: 5-element Vector{Float64}:\n",
" 0.5\n",
" 0.552938681151017\n",
" 0.7089376245893466\n",
" 0.9913594502399236\n",
" 1.3728004409033037"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sol = solve(prob)"
]
},
{
"cell_type": "markdown",
"id": "d1ae611c",
"metadata": {},
"source": [
"## Checking the output\n",
"\n",
"The output from `solve` is _really_ nice. You can \n",
"- check the time steps, `t: 5-element Vector {Float64}`\n",
"- check the solutions, `u: 5-element Vector {Float64}`\n",
"- interpolate between values, `sol(0.5)`\n",
"- [plus more](https://diffeq.sciml.ai/stable/basics/solution/#solution)\n",
"\n",
"For example, here I plot the time vs current value. Then, add the interpolated point at t = 0.5"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "285b277f",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 4,
"metadata": {
"filenames": {
"image/svg+xml": "/home/ryan/Documents/Career_docs/cooperrc-gh-pages/_build/jupyter_execute/Julia-learning/day_04_7_0.svg"
}
},
"output_type": "execute_result"
}
],
"source": [
"Plots.theme(:cooper)\n",
"halftime = 0.5\n",
"plot(sol.t, sol.u, markershape = :utriangle, label = \"solutions\")\n",
"scatter!([halftime], [sol(halftime)], markershape = :square, label = \"interpolated\")\n",
"plot!(xlabel = \"years\", ylabel = \"current value\")"
]
},
{
"cell_type": "markdown",
"id": "56af4d66",
"metadata": {},
"source": [
"## Changing parameters in model\n",
"\n",
"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. \n",
"\n",
"> 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`"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "cb886430",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 5,
"metadata": {
"filenames": {
"image/svg+xml": "/home/ryan/Documents/Career_docs/cooperrc-gh-pages/_build/jupyter_execute/Julia-learning/day_04_9_0.svg"
}
},
"output_type": "execute_result"
}
],
"source": [
"t = range(0, 1, 11)\n",
"\n",
"plts = plot(bg=:white, xlabel = \"time (years)\", ylabel = \"current value\")\n",
"for p in [1.1, 1.01, 0, -1, -2]\n",
" prob = ODEProblem(f,u0,tspan,p)\n",
" sol = solve(prob)\n",
" t = range(0, 1, 11)\n",
" plot!(t, sol(t), label = string(L\"$\\lambda$\",\"= $p\"))\n",
"end\n",
"plts"
]
},
{
"cell_type": "markdown",
"id": "54480387",
"metadata": {},
"source": [
"## Animating the effect of change\n",
"\n",
"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. \n",
"\n",
"Here I build an animation from exponential decay, $\\lambda = -2$ to exponential growth, $\\lambda = 2$ every 0.1 steps. \n",
"\n",
"First, I set up the plot background, labels, and ylimits. Then, use the `@gif` to loop through each value of $\\lambda$ for parameter, `p`."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "bd0f6bbe",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"â”Œ Info: Saved animation to \n",
"â”‚ fn = /home/ryan/Documents/Career_docs/cooperrc-gh-pages/Julia-learning/tmp.gif\n",
"â”” @ Plots /home/ryan/.julia/packages/Plots/D9pfj/src/animation.jl:114\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
"Plots.AnimatedGif(\"/home/ryan/Documents/Career_docs/cooperrc-gh-pages/Julia-learning/tmp.gif\")"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plts = plot(bg=:white, \n",
" xlabel = \"time (years)\", \n",
" ylabel = \"current value\", \n",
" ylims = (0, 4))\n",
"\n",
"@gif for p in -2:0.1:2\n",
" prob = ODEProblem(f,u0,tspan,p)\n",
" sol = solve(prob)\n",
" t = 0:0.1:1\n",
" plot!(t, sol(t), label = \"\")#string(L\"$\\lambda$\",\"= $p\"))\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "0e16d855",
"metadata": {},
"source": [
"## Wrapping up\n",
"\n",
"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."
]
}
],
"metadata": {
"jupytext": {
"formats": "ipynb,md:myst",
"text_representation": {
"extension": ".md",
"format_name": "myst",
"format_version": 0.13,
"jupytext_version": "1.10.3"
}
},
"kernelspec": {
"display_name": "Julia 1.7.2",
"language": "julia",
"name": "julia-1.7"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.7.2"
},
"source_map": [
13,
32,
35,
44,
50,
62,
64,
76,
82,
90,
101,
111,
123
]
},
"nbformat": 4,
"nbformat_minor": 5
}