{
"cells": [
{
"cell_type": "markdown",
"id": "4159be0e",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 2021-10-01 Runge-Kutta methods\n",
"\n",
"## Last time\n",
"\n",
"* Stability for advection-diffusion\n",
"* Cost scaling\n",
"\n",
"## Today\n",
"\n",
"* Accuracy and stability barriers\n",
"* Runge-Kutta methods"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "009bd630",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [
{
"data": {
"text/plain": [
"5×5 SparseMatrixCSC{Float64, Int64} with 15 stored entries:\n",
" -1.25 -0.625 ⋅ ⋅ 1.875\n",
" 1.875 -1.25 -0.625 ⋅ ⋅ \n",
" ⋅ 1.875 -1.25 -0.625 ⋅ \n",
" ⋅ ⋅ 1.875 -1.25 -0.625\n",
" -0.625 ⋅ ⋅ 1.875 -1.25"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots\n",
"using LinearAlgebra\n",
"using SparseArrays\n",
"using Zygote\n",
"\n",
"function vander(x, k=nothing)\n",
" if k === nothing\n",
" k = length(x)\n",
" end\n",
" V = ones(length(x), k)\n",
" for j = 2:k\n",
" V[:, j] = V[:, j-1] .* x\n",
" end\n",
" V\n",
"end\n",
"\n",
"function fdstencil(source, target, k)\n",
" \"kth derivative stencil from source to target\"\n",
" x = source .- target\n",
" V = vander(x)\n",
" rhs = zero(x)'\n",
" rhs[k+1] = factorial(k)\n",
" rhs / V\n",
"end\n",
"\n",
"function newton(residual, jacobian, u0; maxits=20)\n",
" u = u0\n",
" uhist = [copy(u)]\n",
" normhist = []\n",
" for k in 1:maxits\n",
" f = residual(u)\n",
" push!(normhist, norm(f))\n",
" J = jacobian(u)\n",
" delta_u = - J \\ f\n",
" u += delta_u\n",
" push!(uhist, copy(u))\n",
" end\n",
" uhist, normhist\n",
"end\n",
"\n",
"function plot_stability(Rz, title; xlims=(-6, 2), ylims=(-4, 4))\n",
" x = LinRange(xlims[1], xlims[2], 100)\n",
" y = LinRange(ylims[1], ylims[2], 100)\n",
" heatmap(x, y, (x, y) -> abs(Rz(x + 1im*y)), c=:bwr, clims=(0, 2), aspect_ratio=:equal, title=title)\n",
"end\n",
"Rz_theta(z, theta) = (1 + (1 - theta)*z) / (1 - theta*z)\n",
"plot_stability(z -> Rz_theta(z, .5), \"Midpoint method\")\n",
"\n",
"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\n",
"\n",
"function advdiff_matrix(n; kappa=1, wind=1, upwind=0.)\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+1))\n",
" diffuse = [-1, 2, -1] * kappa / dx^2\n",
" advect_upwind = [-1, 1, 0] * wind / dx\n",
" advect_center = [-1, 0, 1] * wind / 2dx\n",
" stencil = -diffuse - upwind * advect_upwind - (1 - upwind) * advect_center\n",
" append!(vals, stencil)\n",
" end\n",
" sparse(rows, cols, vals)\n",
"end\n",
"advdiff_matrix(5, kappa=.1)"
]
},
{
"cell_type": "markdown",
"id": "b9a21ca1",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Stiffness\n",
"\n",
"> Stiff equations are problems for which explicit methods don't work. (Hairer and Wanner, 2002)\n",
"\n",
"* \"stiff\" dates to Curtiss and Hirschfelder (1952)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "4951cb46",
"metadata": {
"cell_style": "split",
"scrolled": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"k = 5\n",
"thist, uhist = ode_theta_linear(-k, [.2], forcing=t -> k*cos(t), tfinal=5, h=.5, theta=1)\n",
"scatter(thist, uhist[1,:])\n",
"plot!(cos)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "792b5615",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "-"
}
},
"outputs": [
{
"data": {
"text/plain": [
"mms_error (generic function with 1 method)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function mms_error(h; theta=.5, k=10)\n",
" u0 = [.2]\n",
" thist, uhist = ode_theta_linear(-k, u0, forcing=t -> k*cos(t), tfinal=3, h=h, theta=theta)\n",
" T = thist[end]\n",
" u_exact = (u0 .- k^2/(k^2+1)) * exp(-k*T) .+ k*(sin(T) + k*cos(T))/(k^2 + 1)\n",
" uhist[1,end] .- u_exact]\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "05b77814",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"hs = .5 .^ (1:8)\n",
"errors = mms_error.(hs, theta=0, k=10)\n",
"plot(hs, norm.(errors), marker=:auto, xscale=:log10, yscale=:log10)\n",
"plot!(hs, hs, label=\"\\$h\\$\", legend=:topleft)\n",
"plot!(hs, hs.^2, label=\"\\$h^2\\$\")"
]
},
{
"cell_type": "markdown",
"id": "b3552de0",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"\n",
"## Discuss: is advection-diffusion stiff?"
]
},
{
"cell_type": "code",
"execution_count": 74,
"id": "e17e0349",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 74,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"theta=0\n",
"n = 20\n",
"dx = 2 / n\n",
"kappa = .01\n",
"lambda_min = -4 * kappa / dx^2\n",
"cfl = 1\n",
"h = min(-2 / lambda_min, cfl * dx)\n",
"\n",
"plot_stability(z -> Rz_theta(z, theta),\n",
" \"\\$\\\\theta=$theta, h=$h, Pe_{\\\\Delta x} = $(dx/kappa)\\$\", xlims=(-2, 1), ylims=(-2, 2))\n",
"ev = eigvals(Matrix(h*advdiff_matrix(n, kappa=kappa, wind=1, upwind=0)))\n",
"scatter!(real(ev), imag(ev))"
]
},
{
"cell_type": "markdown",
"id": "a4ef244b",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Cost scaling"
]
},
{
"cell_type": "markdown",
"id": "c2b44b7c",
"metadata": {},
"source": [
"* Spatial discretization with error $O((\\Delta x)^p)$\n",
"* Time discretization with error $O((\\Delta t)^q)$"
]
},
{
"cell_type": "markdown",
"id": "bb792add",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Runge-Kutta methods"
]
},
{
"cell_type": "markdown",
"id": "5a6de051",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"The methods we have considered thus far can all be expressed as Runge-Kutta methods, which are expressed in terms of $s$ \"stage\" equations (possibly coupled) and a completion formula. For the ODE\n",
"\n",
"$$ \\dot u = f(t, u) $$\n",
"\n",
"the Runge-Kutta method is\n",
"\n",
"$$\\begin{split}\n",
"Y_i = u(t) + h \\sum_j a_{ij} f(t+c_j h, Y_j) \\\\\n",
"u(t+h) = u(t) + h \\sum_j b_j f(t+c_j h, Y_j)\n",
"\\end{split}$$\n",
"where $c$ is a vector of *abscissa*, $A$ is a table of coefficients, and $b$ is a vector of completion weights."
]
},
{
"cell_type": "markdown",
"id": "374ae5b1",
"metadata": {
"cell_style": "split"
},
"source": [
"These coefficients are typically expressed in a Butcher Table\n",
"$$ \\left[ \\begin{array}{c|c} c & A \\\\ \\hline & b^T \\end{array} \\right] = \\left[ \\begin{array}{c|cc}\n",
"c_0 & a_{00} & a_{01} \\\\\n",
"c_1 & a_{10} & a_{11} \\\\\n",
"\\hline\n",
"& b_0 & b_1\n",
"\\end{array} \\right] . $$\n",
"We will see that, for consistency, the abscissa $c$ are always the row sums of $A$ and that $\\sum_i b_i = 1$.\n",
"\n",
"If the matrix $A$ is strictly lower triangular, then the method is **explicit** (does not require solving equations). "
]
},
{
"cell_type": "markdown",
"id": "20f96c5e",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Past methods as Runge-Kutta"
]
},
{
"cell_type": "markdown",
"id": "c2934702",
"metadata": {
"cell_style": "split"
},
"source": [
"* Forward Euler\n",
"$$ \\left[ \\begin{array}{c|cc}\n",
"0 & 0 \\\\\n",
"\\hline\n",
"& 1\n",
"\\end{array} \\right] ,$$\n",
"* Backward Euler\n",
"$$ \\left[ \\begin{array}{c|c}\n",
"1 & 1 \\\\\n",
"\\hline\n",
"& 1\n",
"\\end{array} \\right] ,$$\n",
"* Midpoint\n",
"$$ \\left[ \\begin{array}{c|c}\n",
"\\frac 1 2 & \\frac 1 2 \\\\\n",
"\\hline\n",
"& 1\n",
"\\end{array} \\right]. $$"
]
},
{
"cell_type": "markdown",
"id": "01446736",
"metadata": {
"cell_style": "split"
},
"source": [
"Indeed, the $\\theta$ method is\n",
"$$ \\left[ \\begin{array}{c|c}\n",
"\\theta & \\theta \\\\\n",
"\\hline\n",
"& 1\n",
"\\end{array} \\right] $$\n",
"and an alternative \"endpoint\" variant of $\\theta$ (a generalization of the trapezoid rule) is\n",
"$$ \\left[ \\begin{array}{c|cc}\n",
"0 & 0 & 0 \\\\\n",
"1 & 1-\\theta & \\theta \\\\\n",
"\\hline\n",
"& 1-\\theta & \\theta\n",
"\\end{array} \\right]. $$"
]
},
{
"cell_type": "markdown",
"id": "f48a123a",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Stability function for Runge-Kutta"
]
},
{
"cell_type": "markdown",
"id": "b6e5ae85",
"metadata": {
"cell_style": "split"
},
"source": [
"To develop an algebraic expression for stability in terms of the Butcher Table, we consider the test equation\n",
"\n",
"$$ \\dot u = \\lambda u $$\n",
"\n",
"and apply the RK method to yield\n",
"\n",
"$$ \\begin{split} Y_i = u(0) + h \\sum_j a_{ij} \\lambda Y_j \\\\\n",
"u(h) = u(0) + h \\sum_j b_j \\lambda Y_j \\end{split} $$"
]
},
{
"cell_type": "markdown",
"id": "ab160b5a",
"metadata": {
"cell_style": "split"
},
"source": [
"or, in matrix form,\n",
"\n",
"$$ \\begin{split} Y = \\mathbb 1 u(0) + h \\lambda A Y \\\\\n",
"u(h) = u(0) + h \\lambda b^T Y \\end{split} $$\n",
"\n",
"where $\\mathbb 1$ is a column vector of length $s$ consisting of all ones.\n",
"This reduces to\n",
"$$ u(h) = \\underbrace{\\Big( 1 + h\\lambda b^T (I - h \\lambda A)^{-1} \\mathbb 1 \\Big)}_{R(h\\lambda)} u(0) . $$"
]
},
{
"cell_type": "markdown",
"id": "e9218609",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Plotting stability\n",
"\n",
"$$ u(h) = \\underbrace{\\Big( 1 + h\\lambda b^T (I - h \\lambda A)^{-1} \\mathbb 1 \\Big)}_{R(h\\lambda)} u(0) $$"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "1352b886",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"rk_stability (generic function with 1 method)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"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"
]
},
{
"cell_type": "code",
"execution_count": 77,
"id": "949b6657",
"metadata": {
"cell_style": "split",
"scrolled": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 77,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"theta = .5\n",
"rk_theta(theta) = RKTable([theta], [1])\n",
"rk_theta_endpoint(theta) = RKTable([0 0;1-theta theta], [1-theta, theta])\n",
"table = rk_theta(theta)\n",
"plot_stability(z -> rk_stability(z, table), \"\\$\\\\theta = $theta\\$\")"
]
},
{
"cell_type": "markdown",
"id": "cbd694a7",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Heun's method and RK4"
]
},
{
"cell_type": "code",
"execution_count": 78,
"id": "539b7e41",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 78,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"heun = RKTable([0 0; 1 0], [.5, .5])\n",
"plot_stability(z -> rk_stability(z, heun), \"Heun's method\")"
]
},
{
"cell_type": "code",
"execution_count": 79,
"id": "6b89a362",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 79,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rk4 = RKTable([0 0 0 0; .5 0 0 0; 0 .5 0 0; 0 0 1 0],\n",
"[1, 2, 2, 1] / 6)\n",
"plot_stability(z -> rk_stability(z, rk4), \"RK4\")"
]
},
{
"cell_type": "markdown",
"id": "af1eba76",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# An explicit RK solver\n",
"\n",
"\\begin{align}\n",
"Y_i &= u(t) + h \\sum_j a_{ij} f(t+c_j h, Y_j) &\n",
"u(t+h) &= u(t) + h \\sum_j b_j f(t+c_j h, Y_j)\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "5914460e",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"ode_rk_explicit (generic function with 1 method)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"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": "code",
"execution_count": 83,
"id": "5825e078",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 83,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"linear_oscillator(t, u) = [0 -1; 1 0] * u\n",
"thist, uhist = ode_rk_explicit(linear_oscillator, [1., 0], tfinal=6, h=2)\n",
"plot(thist, uhist', marker=:auto)\n",
"plot!([sin, cos])"
]
},
{
"cell_type": "markdown",
"id": "287e99d7",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Measuring convergence and accuracy"
]
},
{
"cell_type": "code",
"execution_count": 84,
"id": "69c67678",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"linear_oscillator_exact (generic function with 1 method)"
]
},
"execution_count": 84,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function mms_error(h, f, u_exact; table=rk4)\n",
" u0 = u_exact(0)\n",
" thist, uhist = ode_rk_explicit(f, u0, tfinal=3, h=h, table=table)\n",
" T = thist[end]\n",
" uhist[:,end] - u_exact(T)\n",
"end\n",
"\n",
"hs = .5 .^ (-3:8)\n",
"linear_oscillator_exact(t) = [cos(t), sin(t)]"
]
},
{
"cell_type": "code",
"execution_count": 85,
"id": "8d786062",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 85,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"heun_errors = mms_error.(hs, linear_oscillator, linear_oscillator_exact, table=heun)\n",
"rk4_errors = mms_error.(hs, linear_oscillator, linear_oscillator_exact, table=rk4)\n",
"plot(hs, [norm.(heun_errors) norm.(rk4_errors)], label=[\"heun\" \"rk4\"], marker=:auto)\n",
"plot!(hs, [hs hs.^2 hs.^3 hs.^4], label=[\"\\$h^$p\\$\" for p in [1 2 3 4]], legend=:topleft, xscale=:log10, yscale=:log10)"
]
},
{
"cell_type": "markdown",
"id": "08a09b1f",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Work-precision (or accuracy vs cost)\n"
]
},
{
"cell_type": "code",
"execution_count": 37,
"id": "b327e374",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"heun_nfeval = length(heun.c) ./ hs\n",
"rk4_nfeval = length(rk4.c) ./ hs\n",
"plot(heun_nfeval, norm.(heun_errors), marker=:auto, label=\"heun\")\n",
"plot!(rk4_nfeval, norm.(rk4_errors), marker=:auto, label=\"rk4\")\n",
"plot!(xscale=:log10, yscale=:log10)"
]
},
{
"cell_type": "markdown",
"id": "22f9d20a",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Effective RK stability diagrams"
]
},
{
"cell_type": "code",
"execution_count": 40,
"id": "d0dbec9a",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function rk_eff_stability(z, rk)\n",
" s = length(rk.b)\n",
" z = ]\n",
" 1 + z * rk.b' * ((I - z*rk.A) \\ ones(s))\n",
"end\n",
"plot_stability(z -> rk_eff_stability(z, heun), \"Heun's method\", xlims=(-2, 2), ylims=(-2, 2))"
]
},
{
"cell_type": "code",
"execution_count": 39,
"id": "8cbdd3dc",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_stability(z -> rk_eff_stability(z, rk4), \"RK4\", xlims=(-2, 2), ylims=(-2, 2))"
]
},
{
"cell_type": "markdown",
"id": "aa7d76ad",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Runge-Kutta order conditions"
]
},
{
"cell_type": "markdown",
"id": "d22f0a66",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"We consider the autonomous differential equation\n",
"\n",
"$$ \\dot u = f(u) . $$\n",
"\n",
"Higher derivatives of the exact soultion can be computed using the chain rule, e.g.,\n",
"\n",
"\\begin{align*}\n",
"\\ddot u(t) &= f'(u) \\dot u = f'(u) f(u) \\\\\n",
"\\dddot u(t) &= f''(u) f(u) f(u) + f'(u) f'(u) f(u) . \\\\\n",
"\\end{align*}\n",
"\n",
"Note that if $f(u)$ is linear, $f''(u) = 0$.\n",
"Meanwhile, the numerical solution is a function of the time step $h$,\n",
"\n",
"$$\\begin{split}\n",
"Y_i(h) &= u(0) + h \\sum_j a_{ij} f(Y_j) \\\\\n",
"U(h) &= u(0) + h \\sum_j b_j f(Y_j).\n",
"\\end{split}$$"
]
},
{
"cell_type": "markdown",
"id": "714e9ddf",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"We will take the limit $h\\to 0$ and equate derivatives of the numerical solution. First we differentiate the stage equations,\n",
"\n",
"\\begin{split}\n",
"Y_i(0) &= u(0) \\\\\n",
"\\dot Y_i(0) &= \\sum_j a_{ij} f(Y_j) \\\\\n",
"\\ddot Y_i(0) &= 2 \\sum_j a_{ij} \\dot f(Y_j) \\\\\n",
"&= 2 \\sum_j a_{ij} f'(Y_j) \\dot Y_j \\\\\n",
"&= 2\\sum_{j,k} a_{ij} a_{jk} f'(Y_j) f(Y_k) \\\\\n",
"\\dddot Y_i(0) &= 3 \\sum_j a_{ij} \\ddot f (Y_j) \\\\\n",
"&= 3 \\sum_j a_{ij} \\Big( \\sum_k f''(Y_j) \\dot Y_j \\dot Y_k + f'(Y_j) \\ddot Y_j \\Big) \\\\\n",
"&= 3 \\sum_{j,k,\\ell} a_{ij} a_{jk} \\Big( a_{j\\ell} f''(Y_j) f(Y_k) f(Y_\\ell) + 2 a_{k\\ell} f'(Y_j) f'(Y_k) f(Y_\\ell) \\Big)\n",
"\\end{split}"
]
},
{
"cell_type": "markdown",
"id": "cb9aab61",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": ""
}
},
"source": [
"where we have used Liebnitz's formula for the $m$th derivative,\n",
"$$ (h \\phi(h))^{(m)}|_{h=0} = m \\phi^{(m-1)}(0) .$$\n",
"Similar formulas apply for $\\dot U(0)$, $\\ddot U(0)$, and $\\dddot U(0)$, with $b_j$ in place of $a_{ij}$."
]
},
{
"cell_type": "markdown",
"id": "dc019575",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Order conditions are nonlinear algebraic equations"
]
},
{
"cell_type": "markdown",
"id": "f0afb541",
"metadata": {
"slideshow": {
"slide_type": ""
}
},
"source": [
"Equating terms $\\dot u(0) = \\dot U(0)$ yields\n",
"$$ \\sum_j b_j = 1, $$\n",
"equating $\\ddot u(0) = \\ddot U(0)$ yields\n",
"$$ 2 \\sum_{j,k} b_j a_{jk} = 1 , $$\n",
"and equating $\\dddot u(0) = \\dddot U(0)$ yields the two equations\n",
"\\begin{split}\n",
" 3\\sum_{j,k,\\ell} b_j a_{jk} a_{j\\ell} &= 1 \\\\\n",
" 6 \\sum_{j,k,\\ell} b_j a_{jk} a_{k\\ell} &= 1 .\n",
"\\end{split}"
]
},
{
"cell_type": "markdown",
"id": "966e0bca",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Observations\n",
"\n",
"* These are systems of nonlinear equations for the coefficients $a_{ij}$ and $b_j$. There is no guarantee that they have solutions.\n",
"* The number of equations grows rapidly as the order increases.\n",
"\n",
"| | $u^{(1)}$ | $u^{(2)}$ | $u^{(3)}$ | $u^{(4)}$ | $u^{(5)}$ | $u^{(6)}$ | $u^{(7)}$ | $u^{(8)}$ | $u^{(9)}$ | $u^{(10)}$ |\n",
"| ----------- |---|---|---|---|---|----|----|-----|-----|-----|\n",
"| # terms | 1 | 1 | 2 | 4 | 9 | 20 | 48 | 115 | 286 | 719 |\n",
"| cumulative | 1 | 2 | 4 | 8 | 17 | 37 | 85 | 200 | 486 | 1205 |\n",
"\n",
"* Usually the number of order conditions does not exactly match the number of free parameters, meaning that the remaining parameters can be optimized (usually numerically) for different purposes, such as to minimize the leading error terms or to maximize stability in certain regions of the complex plane. Finding globally optimal solutions can be extremely demanding."
]
},
{
"cell_type": "markdown",
"id": "d24b9d83",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* Rooted trees provide a convenient notation\n",
"\n",
"![Trees and elementary differentials up to order 5. From Hairer, Nørsett, and Wanner.](../img/rktrees.png)"
]
},
{
"cell_type": "markdown",
"id": "d8dd49ba",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Theorem (from Hairer, Nørsett, and Wanner)\n",
"\n",
"A Runge-Kutta method is of order $p$ if and only if\n",
"$$ \\gamma(\\mathcal t) \\sum_{j} b_j \\Phi_j(t) = 1 $$\n",
"for all trees $t$ of order $\\le p$.\n",
"\n",
"For a linear autonomous equation\n",
"$$ \\dot u = A u $$\n",
"we only need one additional order condition per order of accuracy because $f'' = 0$.\n",
"These conditions can also be derived by equating derivatives of the stability function $R(z)$ with the exponential $e^z$.\n",
"For a linear non-autonomous equation\n",
"$$ \\dot u = A(t) u + g(t) $$\n",
"or more generally, an autonomous system with quadratic right hand side,\n",
"$$ \\dot u = B (u \\otimes u) + A u + C $$\n",
"where $B$ is a rank 3 tensor, we have $f''' = 0$, thus limiting the number of order conditions."
]
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Julia 1.6.2",
"language": "julia",
"name": "julia-1.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.6.3"
},
"rise": {
"enable_chalkboard": true
}
},
"nbformat": 4,
"nbformat_minor": 5
}