{
"cells": [
{
"cell_type": "markdown",
"id": "4159be0e",
"metadata": {
"hideCode": false,
"hidePrompt": false,
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 2022-04-20 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](http://www.cs.cornell.edu/cv/ResearchPDF/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": 85,
"id": "c1c87f6b",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 85,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f1(t, y; k=5) = -k * (y .- cos(t))\n",
"\n",
"thist, yhist = ode_euler(f1, [1.], tfinal=10, h=.2)\n",
"scatter(thist, yhist[1,:])\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": 83,
"id": "d9a8f791",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 83,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f2(t, y) = [0 1; -1 0] * y\n",
"\n",
"thist, yhist = ode_euler(f2, [0., 1], h=.5, tfinal=10)\n",
"scatter(thist, yhist')\n",
"plot!([cos, sin])"
]
},
{
"cell_type": "code",
"execution_count": 82,
"id": "95eb6628",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Vector{ComplexF64}:\n",
" 0.0 - 1.0im\n",
" 0.0 + 1.0im"
]
},
"execution_count": 82,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eigvals([0 1; -1 0])"
]
},
{
"cell_type": "markdown",
"id": "53e73e58",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Runge-Kutta 4"
]
},
{
"cell_type": "code",
"execution_count": 89,
"id": "747257b9",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 89,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"thist, yhist = ode_rk_explicit(f2, [0., 1], h=.5, tfinal=50)\n",
"scatter(thist, yhist')\n",
"plot!([cos, sin], size=(800, 500))"
]
},
{
"cell_type": "markdown",
"id": "b755e16c",
"metadata": {
"cell_style": "split"
},
"source": [
"* Apparently it is possible to integrate this system using large time steps.\n",
"* This method evaluates $f(y)$ four times per stepso the cost is about equal when the step size $h$ is 4x larger than forward Euler."
]
},
{
"cell_type": "markdown",
"id": "6201b60c",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Linear Stability Analysis\n",
"\n",
"Why did Euler diverge (even if slowly) while RK4 solved this problem accurately?\n",
"And why do both methods diverge if the step size is too large?\n",
"We can understand the convergence of methods by analyzing the test problem\n",
"$$ y' = \\lambda y $$\n",
"for different values of $\\lambda$ in the complex plane.\n",
"One step of the Euler method with step size $h$ maps\n",
"$$ y \\mapsto y + h \\lambda y = \\underbrace{(1 + h \\lambda)}_{R(h \\lambda)} y .$$\n",
"When does this map cause solutions to \"blow up\" and when is it stable?"
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "295efcde",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"plot_stability (generic function with 1 method)"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function plot_stability(Rz, method; xlim=(-3, 2), ylim=(-1.5, 1.5))\n",
" x = xlim[1]:.02:xlim[2]\n",
" y = ylim[1]:.02:ylim[2]\n",
" plot(title=\"Stability: $method\", aspect_ratio=:equal, xlim=xlim, ylim=ylim)\n",
" heatmap!(x, y, (x, y) -> abs(Rz(x + 1im*y)), c=:bwr, clims=(0, 2))\n",
" contour!(x, y, (x, y) -> abs(Rz(x + 1im*y)), color=:black, linewidth=2, levels=[1.])\n",
" plot!(x->0, color=:black, linewidth=1, label=:none)\n",
" plot!([0, 0], [ylim...], color=:black, linewidth=1, label=:none)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 31,
"id": "f233ed2c",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_stability(z -> 1 + z, \"Forward Eulor\")"
]
},
{
"cell_type": "markdown",
"id": "708a2cb7",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Stability for RK4"
]
},
{
"cell_type": "code",
"execution_count": 91,
"id": "82563216",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 91,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_stability(z -> rk_stability(4z, rk4), \"RK4\")\n"
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "653d9295",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_stability(z -> rk_stability(2z, heun), \"Heun's method\")"
]
},
{
"cell_type": "markdown",
"id": "eabf6b5d",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Implicit methods"
]
},
{
"cell_type": "markdown",
"id": "290a4108",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"Recall that forward Euler is the step\n",
"$$ \\tilde y(h) = y(0) + h f(0, y(0)) . $$\n",
"This can be evaluated **explicitly**; all the terms on the right hand side are known so the approximation $\\tilde y(h)$ is computed merely by evaluating the right hand side.\n",
"Let's consider an alternative, **backward Euler** (or \"implicit Euler\"),\n",
"$$ \\tilde y(h) = y(0) + h f(h, \\tilde y(h)) . $$\n",
"This is a (generally) nonlinear equation for $\\tilde y(h)$.\n",
"For the test equation $y' = \\lambda y$, the backward Euler method is\n",
"$$ \\tilde y(h) = y(0) + h \\lambda \\tilde y(h) $$\n",
"or\n",
"$$ \\tilde y(h) = \\underbrace{\\frac{1}{1 - h \\lambda}}_{R(h\\lambda)} y(0) . $$"
]
},
{
"cell_type": "code",
"execution_count": 34,
"id": "c9e9d867",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_stability(z -> 1/(1-z), \"Backward Euler\")"
]
},
{
"cell_type": "markdown",
"id": "f3a71738",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Computing with implicit methods\n",
"\n",
"$$ \\tilde y(h) = y(0) + h f\\big(\\tilde y(h) \\big) $$"
]
},
{
"cell_type": "markdown",
"id": "0ece15d3",
"metadata": {
"cell_style": "split"
},
"source": [
"* Linear solve for linear problem\n",
"* Nonlinear (often Newton) solve for nonlinear problem\n",
"* Need Jacobian or finite differencing"
]
},
{
"cell_type": "code",
"execution_count": 35,
"id": "c479bcf4",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_stability(z -> Rz_theta(z, .5), \"Midpoint method\")"
]
},
{
"cell_type": "markdown",
"id": "6631b25b",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# The $\\theta$ method\n"
]
},
{
"cell_type": "markdown",
"id": "4d01c9a3",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"Forward and backward Euler are bookends of the family known as $\\theta$ methods.\n",
"\n",
"$$ \\tilde u(h) = u(0) + h f\\Big(\\theta h, \\theta\\tilde u(h) + (1-\\theta)u(0) \\Big) $$\n",
"\n",
"which, for linear problems, is solved as\n",
"\n",
"$$ (I - h \\theta A) u(h) = \\Big(I + h (1-\\theta) A \\Big) u(0) . $$\n",
"\n",
"$\\theta=0$ is explicit Euler, $\\theta=1$ is implicit Euler, and $\\theta=1/2$ are the midpoint or trapezoid rules (equivalent for linear problems).\n",
"The stability function is\n",
"$$ R(z) = \\frac{1 + (1-\\theta)z}{1 - \\theta z}. $$"
]
},
{
"cell_type": "code",
"execution_count": 123,
"id": "a035183c",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 123,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Rz_theta(z, theta) = (1 + (1-theta)*z) / (1 - theta*z)\n",
"theta = .6\n",
"plot_stability(z -> Rz_theta(z, theta), \"\\$\\\\theta=$theta\\$\", \n",
" xlim=(-20, 1))"
]
},
{
"cell_type": "markdown",
"id": "85a81669",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# $\\theta$ method for the oscillator"
]
},
{
"cell_type": "code",
"execution_count": 98,
"id": "f9ac0b5e",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"ode_theta_linear (generic function with 1 method)"
]
},
"execution_count": 98,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function ode_theta_linear(A, u0; forcing=zero, tfinal=1, h=0.1, theta=.5)\n",
" u = copy(u0)\n",
" t = 0.\n",
" thist = [t]\n",
" uhist = [u0]\n",
" while t < tfinal\n",
" tnext = min(t+h, tfinal)\n",
" h = tnext - t\n",
" rhs = (I + h*(1-theta)*A) * u .+ h*forcing(t+h*theta)\n",
" u = (I - h*theta*A) \\ rhs\n",
" t = tnext\n",
" push!(thist, t)\n",
" push!(uhist, u)\n",
" end\n",
" thist, hcat(uhist...)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 111,
"id": "8ecc9503",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 111,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Test on oscillator\n",
"A = [0 1; -1 0]\n",
"theta = 1\n",
"thist, uhist = ode_theta_linear(A, [0., 1], h=.2, theta=theta, tfinal=20)\n",
"scatter(thist, uhist')\n",
"plot!([cos, sin])"
]
},
{
"cell_type": "markdown",
"id": "3a83dbf3",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# $\\theta$ method for the cosine decay"
]
},
{
"cell_type": "code",
"execution_count": 109,
"id": "fdf09c98",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 109,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"k = 50\n",
"theta = .5\n",
"thist, uhist = ode_theta_linear(-k, [.2], forcing=t -> k*cos(t), tfinal=5, h=.5, theta=theta)\n",
"scatter(thist, uhist[1,:], title=\"\\$\\\\theta = $theta\\$\")\n",
"plot!(cos, size=(800, 500))"
]
},
{
"cell_type": "markdown",
"id": "c536eb9d",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Stability classes and the $\\theta$ method\n"
]
},
{
"cell_type": "markdown",
"id": "b984e88a",
"metadata": {
"cell_style": "split"
},
"source": [
"## Definition: $A$-stability\n",
"A method is $A$-stable if the stability region\n",
"$$ \\{ z : |R(z)| \\le 1 \\} $$\n",
"contains the entire left half plane $$ \\Re[z] \\le 0 .$$\n",
"This means that the method can take arbitrarily large time steps without becoming unstable (diverging) for any problem that is indeed physically stable."
]
},
{
"cell_type": "markdown",
"id": "232b9cdf",
"metadata": {
"cell_style": "split"
},
"source": [
"## Definition: $L$-stability\n",
"A time integrator with stability function $R(z)$ is $L$-stable if\n",
"$$ \\lim_{z\\to\\infty} R(z) = 0 .$$\n",
"For the $\\theta$ method, we have\n",
"$$ \\lim_{z\\to \\infty} \\frac{1 + (1-\\theta)z}{1 - \\theta z} = \\frac{1-\\theta}{\\theta} . $$\n",
"Evidently only $\\theta=1$ is $L$-stable."
]
},
{
"cell_type": "markdown",
"id": "f8c8e7fa",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Heat equation as linear ODE\n",
"\n",
"* How do different $\\theta \\in [0, 1]$ compare in terms of stability?\n",
"* Are there artifacts even when the solution is stable?"
]
},
{
"cell_type": "code",
"execution_count": 48,
"id": "2ea7e708",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"5×5 SparseMatrixCSC{Float64, Int64} with 15 stored entries:\n",
" -12.5 6.25 â‹… â‹… 6.25\n",
" 6.25 -12.5 6.25 â‹… â‹… \n",
" â‹… 6.25 -12.5 6.25 â‹… \n",
" â‹… â‹… 6.25 -12.5 6.25\n",
" 6.25 â‹… â‹… 6.25 -12.5"
]
},
"execution_count": 48,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using SparseArrays\n",
"\n",
"function heat_matrix(n)\n",
" dx = 2 / n\n",
" rows = [1]\n",
" cols = [1]\n",
" vals = [0.]\n",
" wrap(j) = (j + n - 1) % n + 1\n",
" for i in 1:n\n",
" append!(rows, [i, i, i])\n",
" append!(cols, wrap.([i-1, i, i+1]))\n",
" append!(vals, [1, -2, 1] ./ dx^2)\n",
" end\n",
" sparse(rows, cols, vals)\n",
"end\n",
"heat_matrix(5)"
]
},
{
"cell_type": "code",
"execution_count": 49,
"id": "8b832a79",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1.398294 seconds (4.36 M allocations: 229.282 MiB, 5.36% gc time, 99.92% compilation time)\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 200\n",
"A = heat_matrix(n)\n",
"x = LinRange(-1, 1, n+1)[1:end-1]\n",
"u0 = exp.(-200 * x .^ 2)\n",
"@time thist, uhist = ode_theta_linear(A, u0, h=.1, theta=1, tfinal=1);\n",
"nsteps = size(uhist, 2)\n",
"plot(x, uhist[:, 1:5])"
]
},
{
"cell_type": "markdown",
"id": "0ac0a049",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Advection as a linear ODE"
]
},
{
"cell_type": "code",
"execution_count": 52,
"id": "f73784a3",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [
{
"data": {
"text/plain": [
"5×5 SparseMatrixCSC{Float64, Int64} with 11 stored entries:\n",
" 0.0 -1.25 â‹… â‹… 1.25\n",
" 1.25 â‹… -1.25 â‹… â‹… \n",
" â‹… 1.25 â‹… -1.25 â‹… \n",
" â‹… â‹… 1.25 â‹… -1.25\n",
" -1.25 â‹… â‹… 1.25 â‹… "
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function advect_matrix(n; upwind=false)\n",
" dx = 2 / n\n",
" rows = [1]\n",
" cols = [1]\n",
" vals = [0.]\n",
" wrap(j) = (j + n - 1) % n + 1\n",
" for i in 1:n\n",
" append!(rows, [i, i])\n",
" if upwind\n",
" append!(cols, wrap.([i-1, i]))\n",
" append!(vals, [1., -1] ./ dx)\n",
" else\n",
" append!(cols, wrap.([i-1, i+1]))\n",
" append!(vals, [1., -1] ./ 2dx)\n",
" end\n",
" end\n",
" sparse(rows, cols, vals)\n",
"end\n",
"advect_matrix(5)"
]
},
{
"cell_type": "code",
"execution_count": 54,
"id": "0545bfb7",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.127159 seconds (179.07 k allocations: 11.513 MiB, 13.49% gc time, 83.18% compilation time)\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 50\n",
"A = advect_matrix(n, upwind=false)\n",
"x = LinRange(-1, 1, n+1)[1:end-1]\n",
"u0 = exp.(-9 * x .^ 2)\n",
"@time thist, uhist = ode_theta_linear(A, u0, h=.04, theta=1, tfinal=1.);\n",
"nsteps = size(uhist, 2)\n",
"plot(x, uhist[:, 1:(nsteps÷8):end])"
]
},
{
"cell_type": "markdown",
"id": "4c61e8fb",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Spectrum of operators"
]
},
{
"cell_type": "code",
"execution_count": 53,
"id": "3f736202",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"theta=.5\n",
"h = .1\n",
"plot_stability(z -> Rz_theta(z, theta), \"\\$\\\\theta=$theta, h=$h\\$\")\n",
"ev = eigvals(Matrix(h*advect_matrix(20, upwind=true)))\n",
"scatter!(real(ev), imag(ev))"
]
},
{
"cell_type": "code",
"execution_count": 51,
"id": "833e0984",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 51,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"theta=.5\n",
"h = .2 / 4\n",
"plot_stability(z -> Rz_theta(z, theta), \"\\$\\\\theta=$theta, h=$h\\$\")\n",
"ev = eigvals(Matrix(h*heat_matrix(20)))\n",
"scatter!(real(ev), imag(ev))"
]
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"celltoolbar": "Slideshow",
"hide_code_all_hidden": false,
"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"
},
"rise": {
"enable_chalkboard": true
}
},
"nbformat": 4,
"nbformat_minor": 5
}