04 differential equations - 01 introduction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[![Binder](https://mybinder.org/badge_logo.svg)](https://notebooks.gesis.org/binder/v2/gh/jolin-io/fall-in-love-with-julia/main?filepath=04%20differential%20equations%20-%2001%20introduction.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DifferentialEquations.jl\r\n",
    "\r\n",
    "DifferentialEquations is a suite for numerically solving differential equations written in Julia and available for use in Julia, Python, and R. \r\n",
    "\r\n",
    "The well-optimized DifferentialEquations solvers benchmark as the some of the fastest implementations, using classic algorithms and ones from recent research which routinely outperform the \"standard\" C/Fortran methods, and include algorithms optimized for high-precision and HPC applications.\r\n",
    "\r\n",
    "Equations included within DifferentialEquations.jl are for example:\r\n",
    "- Discrete equations (function maps, discrete stochastic (Gillespie/Markov) simulations)\r\n",
    "- Ordinary differential equations (ODEs)\r\n",
    "- Split and Partitioned ODEs (Symplectic integrators, IMEX Methods)\r\n",
    "- Stochastic ordinary differential equations (SODEs or SDEs)\r\n",
    "- Stochastic differential-algebraic equations (SDAEs)\r\n",
    "- Random differential equations (RODEs or RDEs)\r\n",
    "- Differential algebraic equations (DAEs)\r\n",
    "- Delay differential equations (DDEs)\r\n",
    "- Neutral, retarded, and algebraic delay differential equations (NDDEs, RDDEs, and DDAEs)\r\n",
    "- Stochastic delay differential equations (SDDEs)\r\n",
    "- Experimental support for stochastic neutral, retarded, and algebraic delay differential equations (SNDDEs, SRDDEs, and SDDAEs)\r\n",
    "- Mixed discrete and continuous equations (Hybrid Equations, Jump Diffusions)\r\n",
    "- (Stochastic) partial differential equations ((S)PDEs) (with both finite difference and finite element methods)\r\n",
    "\r\n",
    "For more information look at the documentation https://diffeq.sciml.ai/stable/.\r\n",
    "\r\n",
    "The basics are best explained using Ordinary Differential Equations (ODEs)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\r\n",
    "# Example 1: Solving Scalar Equations\r\n",
    "\r\n",
    "In this example we will solve the equation\r\n",
    "\r\n",
    "$$\\frac{du}{dt} = f(u,p,t)$$\r\n",
    "\r\n",
    "on the time interval $t\\in[0,1]$ where $f(u,p,t)=αu$. We know by Calculus that the solution to this equation is $u(t)=u₀exp⁡(αt)$.\r\n",
    "\r\n",
    "The general workflow is to define a problem, solve the problem, and then analyze the solution. The full code for solving this problem is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using DifferentialEquations\r\n",
    "f(u,p,t) = 1.01*u\r\n",
    "u0 = 1/2\r\n",
    "tspan = (0.0,1.0)\r\n",
    "prob = ODEProblem(f,u0,tspan)\r\n",
    "sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots\r\n",
    "plotly()\r\n",
    "plot(sol, linewidth=5, title=\"Solution to the linear ODE with a thick line\",\r\n",
    "     xaxis=\"Time (t)\", yaxis=\"u(t) (in μm)\", label=\"My Thick Line!\") # legend=false\r\n",
    "plot!(sol.t, t->0.5*exp(1.01t), lw=3, ls=:dash, label=\"True Solution!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol.t[5], sol[5]  # access solution at time points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol(0.13)  # automated interpolation of the solution for any time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 2: Solving Systems of Equations\r\n",
    "\r\n",
    "In this example we will solve the Lorenz equations:\r\n",
    "\r\n",
    "$$\\begin{aligned} \\frac{dx}{dt} &= σ(y-x) \\\\ \\frac{dy}{dt} &= x(ρ-z) - y \\\\ \\frac{dz}{dt} &= xy - βz \\\\ \\end{aligned}$$\r\n",
    "\r\n",
    "Defining your ODE function to be in-place updating can have performance benefits. What this means is that, instead of writing a function which outputs its solution, you write a function which updates a vector that is designated to hold the solution. By doing this, DifferentialEquations.jl's solver packages are able to reduce the amount of array allocations and achieve better performance.\r\n",
    "\r\n",
    "The way we do this is we simply write the output to the 1st input of the function. For example, our Lorenz equation problem would be defined by the function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function lorenz!(du,u,p,t)\r\n",
    " du[1] = 10.0*(u[2]-u[1])\r\n",
    " du[2] = u[1]*(28.0-u[3]) - u[2]\r\n",
    " du[3] = u[1]*u[2] - (8/3)*u[3]\r\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and then we can use this function in a problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u0 = [1.0;0.0;0.0]\r\n",
    "tspan = (0.0,100.0)\r\n",
    "prob = ODEProblem(lorenz!,u0,tspan)\r\n",
    "sol = solve(prob)\r\n",
    "plot(sol,vars=(1,2,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining Parameterized Functions\r\n",
    "\r\n",
    "In many cases you may want to explicitly have parameters associated with your differential equations. This can be used by things like [parameter estimation routines](https://diffeq.sciml.ai/stable/analysis/parameter_estimation/#parameter_estimation). In this case, you use the p values via the syntax:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function parameterized_lorenz!(du,u,p,t)\r\n",
    "  x,y,z = u\r\n",
    "  σ,ρ,β = p\r\n",
    "  du[1] = dx = σ*(y-x)\r\n",
    "  du[2] = dy = x*(ρ-z) - y\r\n",
    "  du[3] = dz = x*y - β*z\r\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and then we add the parameters to the ODEProblem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u0 = [1.0,0.0,0.0]\r\n",
    "tspan = (0.0,100.0)\r\n",
    "p = [10.0,28.0,8/3]\r\n",
    "prob = ODEProblem(parameterized_lorenz!,u0,tspan,p)\r\n",
    "sol = solve(prob)\r\n",
    "plot(sol,vars=(1,2,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the type for the parameters p can be anything: you can use arrays, static arrays, named tuples, etc. to enclose your parameters in a way that is sensible for your problem.\r\n",
    "\r\n",
    "Since the parameters exist within the function, functions defined in this manner can also be used for sensitivity analysis, parameter estimation routines, and bifurcation plotting. This makes DifferentialEquations.jl a full-stop solution for differential equation analysis which also achieves high performance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 3: Solving Nonhomogeneous Equations using Parameterized Functions\r\n",
    "\r\n",
    "Parameterized functions can also be used for building nonhomogeneous ordinary differential equations (these are also referred to as ODEs with nonzero right-hand sides). They are frequently used as models for dynamical systems with external (in general time-varying) inputs. As an example, consider a model of a pendulum consisting of a slender rod of length $l$ and mass $m$:\r\n",
    "\r\n",
    "$$\\begin{aligned} \\frac{\\mathrm{d}\\theta(t)}{\\mathrm{d}t} &= \\omega(t)\\\\ \\frac{\\mathrm{d}\\omega(t)}{\\mathrm{d}t} &= - \\frac{3}{2}\\frac{g}{l}\\sin\\theta(t) + \\frac{3}{ml^2}M(t) \\end{aligned}$$\r\n",
    "\r\n",
    "where $θ$ and $ω$ are the angular deviation of the pendulum from the vertical (hanging) orientation and the angular rate, respectively, $M$ is an external torque (developed, say, by a wind or a motor), and finally, $g$ stands for gravitional acceleration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = 1.0                             # length [m]\r\n",
    "m = 1.0                             # mass[m]\r\n",
    "g = 9.81                            # gravitational acceleration [m/s²]\r\n",
    "\r\n",
    "function pendulum!(du,u,p,t)\r\n",
    "    θ, ω = u\r\n",
    "    du[1] = dθ = ω                  # θ'(t) = ω(t)\r\n",
    "    du[2] = dω = -3g/(2l)*sin(θ) + 3/(m*l^2)*p(t)  # ω'(t) = -3g/(2l) sin θ(t) + 3/(ml^2)M(t)\r\n",
    "end\r\n",
    "\r\n",
    "θ₀ = 0.01                           # initial angular deflection [rad]\r\n",
    "ω₀ = 0.0                            # initial angular velocity [rad/s]\r\n",
    "u₀ = [θ₀, ω₀]                       # initial state vector\r\n",
    "tspan = (0.0,10.0)                  # time interval\r\n",
    "\r\n",
    "M = t->0.1sin(t)                    # external torque [Nm]\r\n",
    "\r\n",
    "prob = ODEProblem(pendulum!,u₀,tspan,M)\r\n",
    "sol = solve(prob)\r\n",
    "\r\n",
    "plot(sol,linewidth=2,xaxis=\"t\",label=[\"θ [rad]\" \"ω [rad/s]\"],layout=(2,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how the external **time-varying** torque ``M`` is introduced as a **parameter** in the `pendulum!` function. Indeed, as a general principle the parameters can be any type; here we specify M as time-varying by representing it by a function, which is expressed by appending the dependence on time `(t)` to the name of the parameter.\r\n",
    "\r\n",
    "Note also that, in contrast with the time-varying parameter, the (vector of) state variables `u`, which is generally also time-varying, is always used without the explicit dependence on time `(t)`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 4: Using Other Types for Systems of Equations\r\n",
    "\r\n",
    "DifferentialEquations.jl can handle many different dependent variable types (generally, anything with a linear index should work!). So instead of solving a vector equation, let's let ``u`` be a matrix! To do this, we simply need to have ``u0`` be a matrix, and define ``f`` such that it takes in a matrix and outputs a matrix. We can define a matrix of linear ODEs as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A  = [1. 0  0 -5\r\n",
    "      4 -2  4 -3\r\n",
    "     -4  0  0  1\r\n",
    "      5 -2  2  3]\r\n",
    "u0 = rand(4,2)\r\n",
    "tspan = (0.0,1.0)\r\n",
    "f(u,p,t) = A*u\r\n",
    "prob = ODEProblem(f,u0,tspan)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here our ODE is on a 4x2 matrix, and the ODE is the linear system defined by multiplication by A. To solve the ODE, we do the same steps as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = solve(prob)\r\n",
    "plot(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can instead use the in-place form by using Julia's in-place matrix multiplication function ``mul!``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\r\n",
    "f(du,u,p,t) = mul!(du,A,u)\r\n",
    "prob = ODEProblem(f,u0,tspan)\r\n",
    "sol = solve(prob)\r\n",
    "plot(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Additionally, we can use non-traditional array types as well. For example, ``StaticArrays.jl`` offers immutable arrays which are stack-allocated, meaning that their usage does not require any (slow) heap-allocations that arrays normally have. This means that they can be used to solve the same problem as above, with the only change being the type for the initial condition and constants:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using StaticArrays\r\n",
    "A  = @SMatrix [ 1.0  0.0 0.0 -5.0\r\n",
    "                4.0 -2.0 4.0 -3.0\r\n",
    "               -4.0  0.0 0.0  1.0\r\n",
    "                5.0 -2.0 2.0  3.0 ]\r\n",
    "u0 = @SMatrix rand(4,2)\r\n",
    "tspan = (0.0,1.0)\r\n",
    "f′(u,p,t) = A*u  # we need a new function here, as otherwise ODEProblem will think we can use the mutation version of f defined above. Give it a try\r\n",
    "prob = ODEProblem(f′,u0,tspan)\r\n",
    "sol = solve(prob)\r\n",
    "plot(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the analysis tools generalize over to systems of equations as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol[4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "still returns the solution at the fourth timestep. It also indexes into the array as well. The last value is the timestep, and the beginning values are for the component. This means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol[5,3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "is the value of the 5th component (by linear indexing) at the 3rd timepoint, or"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol[2,1,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "is the timeseries for the component which is the 2nd row and 1st column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CartesianIndices(A)[5], LinearIndices(A)[2,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Going Beyond ODEs: How to Use the Documentation\r\n",
    "\r\n",
    "Not everything can be covered in the tutorials. Instead, this tutorial will end by pointing you in the directions for the next steps.\r\n",
    "Common API for Defining, Solving, and Plotting\r\n",
    "\r\n",
    "One feature of DifferentialEquations.jl is that this pattern for solving equations is conserved across the different types of differential equations. Every equation has a problem type, a solution type, and the same solution handling (+ plotting) setup. Thus the solver and plotting commands in the **Basics** section applies to all sorts of equations, like stochastic differential equations and delay differential equations. Each of these different problem types are defined in the **Problem Types** section of the docs. Every associated solver algorithm is detailed in the **Solver Algorithms** section, sorted by problem type. The same steps for ODEs can then be used for the analysis of the solution.\r\n",
    "Additional Features and Analysis Tools\r\n",
    "\r\n",
    "In many cases, the common workflow only starts with solving the differential equation. Many common setups have built-in solutions in DifferentialEquations.jl. For example, check out the features for:\r\n",
    "\r\n",
    "* [Handling, parallelizing, and analyzing large Ensemble experiments](https://diffeq.sciml.ai/stable/features/ensemble/#ensemble)\r\n",
    "* [Saving the output to tabular formats like DataFrames and CSVs](https://diffeq.sciml.ai/stable/features/io/#io)\r\n",
    "* [Event handling](https://diffeq.sciml.ai/stable/features/callback_functions/#callbacks)\r\n",
    "* [Parameter estimation (inverse problems)](https://diffeq.sciml.ai/stable/analysis/parameter_estimation/#parameter_estimation)\r\n",
    "* [Quantification of numerical uncertainty and error](https://diffeq.sciml.ai/stable/analysis/uncertainty_quantification/#uncertainty_quantification)\r\n",
    "\r\n",
    "Many more are defined in the relevant sections of the docs. Please explore the rest of the documentation, including tutorials for getting started with other types of equations. In addition, to get help, please either [file an issue at the main repository](https://github.com/JuliaDiffEq/DifferentialEquations.jl) or [come have an informal discussion at our Gitter chatroom](https://gitter.im/JuliaDiffEq/Lobby)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.1",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "name": "julia",
   "version": ""
  },
  "orig_nbformat": 2
 },
 "nbformat": 4,
 "nbformat_minor": 2
}