{
"cells": [
{
"cell_type": "markdown",
"id": "4159be0e",
"metadata": {
"hideCode": false,
"hidePrompt": false,
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 2023-04-21 Stiff equations\n",
"\n",
"## Last time\n",
"\n",
"* Notes on integration\n",
"* Ordinary differential equations (ODE)\n",
"* Explicit methods for solving ODE\n",
"* Stability\n",
"\n",
"## Today\n",
"\n",
"* Implicit and explicit methods\n",
"* Stiff equations\n",
"* PDE examples"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "eb781a7a",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [
{
"data": {
"text/plain": [
"ode_rk_explicit (generic function with 1 method)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using LinearAlgebra\n",
"using Plots\n",
"default(linewidth=4, legendfontsize=12)\n",
"\n",
"struct RKTable\n",
" A::Matrix\n",
" b::Vector\n",
" c::Vector\n",
" function RKTable(A, b)\n",
" s = length(b)\n",
" A = reshape(A, s, s)\n",
" c = vec(sum(A, dims=2))\n",
" new(A, b, c)\n",
" end\n",
"end\n",
"\n",
"function rk_stability(z, rk)\n",
" s = length(rk.b)\n",
" 1 + z * rk.b' * ((I - z*rk.A) \\ ones(s))\n",
"end\n",
"\n",
"rk4 = RKTable([0 0 0 0; .5 0 0 0; 0 .5 0 0; 0 0 1 0], [1, 2, 2, 1] / 6)\n",
"heun = RKTable([0 0; 1 0], [.5, .5])\n",
"Rz_theta(z, theta) = (1 + (1 - theta)*z) / (1 - theta*z)\n",
"\n",
"function ode_rk_explicit(f, u0; tfinal=1, h=0.1, table=rk4)\n",
" u = copy(u0)\n",
" t = 0.\n",
" n, s = length(u), length(table.c)\n",
" fY = zeros(n, s)\n",
" thist = [t]\n",
" uhist = [u0]\n",
" while t < tfinal\n",
" tnext = min(t+h, tfinal)\n",
" h = tnext - t\n",
" for i in 1:s\n",
" ti = t + h * table.c[i]\n",
" Yi = u + h * sum(fY[:,1:i-1] * table.A[i,1:i-1], dims=2)\n",
" fY[:,i] = f(ti, Yi)\n",
" end\n",
" u += h * fY * table.b\n",
" t = tnext\n",
" push!(thist, t)\n",
" push!(uhist, u)\n",
" end\n",
" thist, hcat(uhist...)\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "1bcd691a",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Ordinary Differential Equations\n",
"\n",
"Given initial condition $y_0 = y(t=0)$, find $y(t)$ for $t > 0$ that satisfies\n",
"\n",
"$$ y' \\equiv \\dot y \\equiv \\frac{\\partial y}{\\partial t} = f(t, y) $$"
]
},
{
"cell_type": "markdown",
"id": "f0fb11e0",
"metadata": {
"cell_style": "split"
},
"source": [
"| Application | $y$ | $f$ |\n",
"| --- | --- | --- |\n",
"| Orbital dynamics | position, momentum | conservation of momentum|\n",
"| Chemical reactions | concentration | conservation of atoms |\n",
"| Epidemiology | infected/recovered population | transmission and recovery |"
]
},
{
"cell_type": "markdown",
"id": "e8089acc",
"metadata": {
"cell_style": "split"
},
"source": [
"* $y$ can be a scalar or a vector"
]
},
{
"cell_type": "markdown",
"id": "1858f407",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Solving differential equations\n"
]
},
{
"cell_type": "markdown",
"id": "87c81fd0",
"metadata": {
"cell_style": "split"
},
"source": [
"## Linear equations\n",
"\n",
"$$ y' = A(t) y + \\text{source}(t)$$\n",
"\n",
"* Autonomous if $A(t) = A$ and source independent of $t$\n",
"\n",
"* Suppose $y$ and $a = A$ are scalars: $y(t) = e^{at} y_0$"
]
},
{
"cell_type": "markdown",
"id": "1d2ffb98",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## Can do the same for systems\n",
"\n",
"$$ y(t) = e^{A t} y_0 $$\n",
"\n",
"### What does it mean to exponentiate a matrix?\n",
"\n",
"Taylor series!\n",
"\n",
"$$ e^A = 1 + A + \\frac{A^2}{2} + \\frac{A^3}{3!} + \\dotsb $$\n",
"and there are many [practical ways to compute it](https://www.cs.cornell.edu/cv/files/2021/10/19ways.pdf).\n"
]
},
{
"cell_type": "markdown",
"id": "31469d4e",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## Question\n",
"Suppose that the diagonalization $A = X \\Lambda X^{-1}$ exists and derive a finite expression for the matrix exponential using the scalar `exp` function."
]
},
{
"cell_type": "markdown",
"id": "ad344c00",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Forward Euler Method\n",
"The simplest method for solving $y'(t) = f(t,y)$ is\n",
"to use numerical differentiation to write\n",
"$$ y' \\approx \\frac{y(h) - y(0)}{h} $$\n",
"which yields the solution estimate\n",
"$$ \\tilde y(h) = y(0) + h f(0, y(0)) $$\n",
"where $h$ is the step size.\n",
"Let's try this on a scalar problem\n",
"$$ y' = -k (y - \\cos t) $$\n",
"where $k$ is a parameter controlling the rate at which the solution $u(t)$ is pulled toward the curve $\\cos t$."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "cde76654",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"ode_euler (generic function with 1 method)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function ode_euler(f, y0; tfinal=10., h=0.1)\n",
" y = copy(y0)\n",
" t = 0.\n",
" thist = [t]\n",
" yhist = [y0]\n",
" while t < tfinal\n",
" tnext = min(t+h, tfinal)\n",
" h = tnext - t\n",
" y += h * f(t, y)\n",
" t = tnext\n",
" push!(thist, t)\n",
" push!(yhist, y)\n",
" end\n",
" thist, hcat(yhist...)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 117,
"id": "c1c87f6b",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 117,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f1(t, y; k=10) = -k * (y .- cos(t))\n",
"\n",
"thist, yhist = ode_euler(f1, [1.], tfinal=20, h=.2)\n",
"plot(thist, yhist[1,:], marker=:circle)\n",
"plot!(cos)"
]
},
{
"cell_type": "markdown",
"id": "eedad8e0",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Forward Euler on a linear system\n",
"\n",
"$$ \\begin{bmatrix} y_1 \\\\ y_2 \\end{bmatrix}' = \\begin{bmatrix} 0 & 1 \\\\ -1 & 0 \\end{bmatrix} \\begin{bmatrix} y_1 \\\\ y_2 \\end{bmatrix}$$"
]
},
{
"cell_type": "code",
"execution_count": 118,
"id": "d9a8f791",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"