{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Theory\n",
"\n",
"Some of the rates related formulas are from Leif Andersen's 3-volume book. "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import IPython; IPython.display.HTML('''''')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Automatic Differentiation Example\n",
"\n",
"* 把數學函數全部重新定義一遍，output type 不用 float 而是用一個新的 class（Variable）同時存下數值（data）和導數（grad）\n",
" * 等於把本來的 $\\mathbb R \\mapsto \\mathbb R$ 函數都變成 $\\mathbb R \\mapsto \\mathbb R^2$，同時輸出點值和導數\n",
"* Variable class 超載所有 arithmetic operator 把微分規則寫在裡面（如萊布尼茲公式）\n",
"* 在自定義的數學函數裡手動給他們的導數，例如\n",
"\\begin{align*}\n",
"\\frac{d}{dx}\\sin(x) = \\cos(x),\\qquad\\frac{d}{dx}\\cos(x) = -\\sin(x).\n",
"\\end{align*}\n",
"* Nothing is symbolic"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"data: -1.0717958634011433e-07, grad: 1.9999999999999885"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import math\n",
"\n",
"class Variable:\n",
" \n",
" def __init__(self, data, grad):\n",
" self.data = data\n",
" self.grad = grad\n",
" \n",
" def __add__(self, other):\n",
" res_data = self.data + other.data\n",
" res_grad = self.grad + other.grad\n",
" \n",
" return Variable(res_data, res_grad)\n",
"\n",
" def __mul__(self, other):\n",
" res_data = self.data * other.data\n",
" res_grad = self.data * other.grad + self.grad * other.data # Leibniz rule\n",
"\n",
" return Variable(res_data, res_grad)\n",
"\n",
" def __repr__(self):\n",
" return f'data: {self.data}, grad: {self.grad}'\n",
"\n",
"def sin(x):\n",
" return Variable(math.sin(x), math.cos(x))\n",
"\n",
"def cos(x):\n",
" return Variable(math.cos(x), -math.sin(x))\n",
"\n",
"def f(x):\n",
" return sin(x)*cos(x) + cos(x)*sin(x) # equals sin(2x) with derivative 2 cos(2x)\n",
"\n",
"f(3.1415926) # expecting point value: sin(2*pi) = 0 and derivative 2*cos(2*pi) = 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Crank-Nicolson Summary\n",
"\n",
"To solve $u_t = u_{xx}$, first discretize $x$ to get the system of ODEs \n",
"\\begin{align*}\n",
"\\dot U = \\delta_x^+\\delta_x^- U, \n",
"\\end{align*}\n",
"where the operators are defined as\n",
"\\begin{align*}\n",
"\\delta_x^+ &= \\frac{S-I}{\\Delta x},\\\\\n",
"\\delta_x^- &= \\frac{I-S^{-1}}{\\Delta x},\\\\\n",
"\\delta_x^+\\delta_x^- &= \\frac{S-2I+S^{-1}}{(\\Delta x)^2}, \n",
"\\end{align*}\n",
"and $U$ is a vector function of $t$: \n",
"\\begin{align*}\n",
"U(t) = \\begin{pmatrix}\n",
"U_1(t)\\\\\n",
"U_2(t)\\\\\n",
"U_3(t)\\\\\n",
"\\vdots\\\\\n",
"U_{m-1}(t)\\\\\n",
"\\end{pmatrix}. \n",
"\\end{align*}\n",
"In matrix form, the operators are \n",
"\\begin{align*}\n",
"SU &= \\begin{pmatrix}\n",
"U_2\\\\\n",
"U_3\\\\\n",
"U_4\\\\\n",
"\\vdots\\\\\n",
"U_{m}\\\\\n",
"\\end{pmatrix} = \\begin{pmatrix}\n",
"0 & 1 &&&\\\\\n",
"& 0 & 1 &&\\\\\n",
"&& 0 & 1 &\\\\\n",
"&&& \\ddots \\\\\n",
"&&&& 0\\\\\n",
"\\end{pmatrix}U + \\begin{pmatrix}\n",
"0\\\\\n",
"0\\\\\n",
"0\\\\\n",
"\\vdots\\\\\n",
"\\beta_R\\\\\n",
"\\end{pmatrix}, \\\\\n",
"S^{-1}U &= \\begin{pmatrix}\n",
"U_0\\\\\n",
"U_1\\\\\n",
"U_2\\\\\n",
"\\vdots\\\\\n",
"U_{m-2}\\\\\n",
"\\end{pmatrix} = \\begin{pmatrix}\n",
"0 &&&&\\\\\n",
"1 & 0 &&&\\\\\n",
"& 1 & 0 &&\\\\\n",
"&&& \\ddots \\\\\n",
"&&& 1 & 0\\\\\n",
"\\end{pmatrix}U + \\begin{pmatrix}\n",
"\\beta_L\\\\\n",
"0\\\\\n",
"0\\\\\n",
"\\vdots\\\\\n",
"0\\\\\n",
"\\end{pmatrix}, \n",
"\\end{align*}\n",
"where $\\beta_L = U_0$ and $\\beta_R = U_m$ are Dirichlet boundary conditions. \n",
"Now discretize $t$ and apply the trapezoidal rule\n",
"\\begin{align*}\n",
"U^{n+1} &= U^{n} + \\int_{t_n}^{t_{n+1}} \\dot U\\,dt\\\\\n",
"&= U^{n} + \\int_{t_n}^{t_{n+1}} \\delta_x^+\\delta_x^- U\\,dt\\\\\n",
"&\\approx U^{n} + \\Delta t\\delta_x^+\\delta_x^- \\left(\\frac{U^{n}+U^{n+1}}{2}\\right). \n",
"\\end{align*}\n",
"So the Crank-Nicolson method is \n",
"\\begin{align*}\n",
"\\left(I - \\frac{\\Delta t}{2}\\delta_x^+\\delta_x^- \\right)U^{n+1} = \\left(I + \\frac{\\Delta t}{2}\\delta_x^+\\delta_x^- \\right)U^{n}. \n",
"\\end{align*}\n",
"To write in matrix form, expand the operator to get\n",
"\\begin{align*}\n",
"\\left(I - \\alpha\\left(S-2I+S^{-1}\\right) \\right)U^{n+1} = \\left(I + \\alpha\\left(S-2I+S^{-1}\\right) \\right)U^{n}, \n",
"\\end{align*}\n",
"where\n",
"\\begin{align*}\n",
"\\alpha = \\frac{\\Delta t}{2(\\Delta x)^2}. \n",
"\\end{align*}\n",
"Rearrange to get\n",
"\\begin{align*}\n",
"\\left(-\\alpha S + (1+2\\alpha)I -\\alpha S^{-1} \\right)U^{n+1} = \\left(\\alpha S + (1-2\\alpha)I + \\alpha S^{-1} \\right)U^{n}. \n",
"\\end{align*}\n",
"Thus the scheme in matrix form is\n",
"\\begin{align*}\n",
"&\n",
"\\begin{pmatrix}\n",
"1+2\\alpha & -\\alpha&&&&\\\\\n",
"-\\alpha & 1+2\\alpha & -\\alpha&&&\\\\\n",
"& -\\alpha & 1+2\\alpha & -\\alpha&&\\\\\n",
"&&& \\ddots &\\\\\n",
"&&& -\\alpha & 1+2\\alpha & -\\alpha \\\\\n",
"&&&& -\\alpha & 1+2\\alpha\n",
"\\end{pmatrix}U^{n+1} + \\begin{pmatrix}\n",
"-\\alpha\\beta_L^{n+1}\\\\\n",
"0\\\\\n",
"0\\\\\n",
"\\vdots\\\\\n",
"0\\\\\n",
"-\\alpha\\beta_R^{n+1}\\\\\n",
"\\end{pmatrix} \\\\\n",
"&= \n",
"\\begin{pmatrix}\n",
"1-2\\alpha & \\alpha&&&&\\\\\n",
"\\alpha & 1-2\\alpha & \\alpha&&&\\\\\n",
"& \\alpha & 1-2\\alpha & \\alpha&&\\\\\n",
"&&& \\ddots &\\\\\n",
"&&& \\alpha & 1-2\\alpha & \\alpha \\\\\n",
"&&&& \\alpha & 1-2\\alpha\n",
"\\end{pmatrix}U^{n} + \\begin{pmatrix}\n",
"\\alpha\\beta_L^{n}\\\\\n",
"0\\\\\n",
"0\\\\\n",
"\\vdots\\\\\n",
"0\\\\\n",
"\\alpha\\beta_R^{n}\\\\\n",
"\\end{pmatrix}. \n",
"\\end{align*}\n",
"Note that the solver can only go forward in time: One can only solve the unknown $U^{n+1}$ given $U^{n}$, not the other way around, or otherwise the scheme will become unconditionally unstable. \n",
"\n",
"### Neumann Boundary Conditions\n",
"When Neumann boundary conditions $\\beta_L^\\prime, \\beta_R^\\prime$ are given instead of Dirichlet, we set\n",
"\\begin{align*}\n",
"&-\\frac{3}{2\\Delta x} U_0 + \\frac{4}{2\\Delta x} U_1 -\\frac{1}{2\\Delta x} U_2 = \\beta_L^\\prime, \\\\\n",
"&\\frac{3}{2\\Delta x} U_m - \\frac{4}{2\\Delta x} U_{m-1} +\\frac{1}{2\\Delta x} U_{m-2} = \\beta_R^\\prime, \n",
"\\end{align*}\n",
"where the left hand side are 2nd order approximation of the $u_x$ values at both boundaries, derived from the method of undetermined coefficients. This can be computed by, for example, the [Finite Difference Coefficients Calculator](https://web.media.mit.edu/~crtaylor/calculator.html). Equivalently, we can write\n",
"\\begin{align*}\n",
"U_0 &= -\\frac{2\\Delta x}{3}\\beta_L^\\prime + \\frac{4}{3} U_1 - \\frac{1}{3} U_2, \\\\\n",
"U_m &= \\frac{2\\Delta x}{3}\\beta_R^\\prime + \\frac{4}{3} U_{m-1} - \\frac{1}{3} U_{m-2}. \n",
"\\end{align*}\n",
"Replacing all $\\beta_L$ and $\\beta_R$ in the above linear system by these new expressions of $U_0$ and $U_m$, we obtain the following linear system with Neumann BCs: \n",
"\\begin{align*}\n",
"&\n",
"\\begin{pmatrix}\n",
"1+2\\alpha-\\frac{4}{3}\\alpha & -\\alpha+\\frac{1}{3}\\alpha&&&&\\\\\n",
"-\\alpha & 1+2\\alpha & -\\alpha&&&\\\\\n",
"& -\\alpha & 1+2\\alpha & -\\alpha&&\\\\\n",
"&&& \\ddots &\\\\\n",
"&&& -\\alpha & 1+2\\alpha & -\\alpha \\\\\n",
"&&&& -\\alpha+\\frac{1}{3}\\alpha & 1+2\\alpha-\\frac{4}{3}\\alpha\n",
"\\end{pmatrix}U^{n+1} + \\begin{pmatrix}\n",
"\\frac{2\\alpha\\Delta x}{3}(\\beta_L^\\prime)^{n+1}\\\\\n",
"0\\\\\n",
"0\\\\\n",
"\\vdots\\\\\n",
"0\\\\\n",
"-\\frac{2\\alpha\\Delta x}{3}(\\beta_R^\\prime)^{n+1}\\\\\n",
"\\end{pmatrix} \\\\\n",
"&= \n",
"\\begin{pmatrix}\n",
"1-2\\alpha +\\frac{4}{3}\\alpha & \\alpha-\\frac{1}{3}\\alpha&&&&\\\\\n",
"\\alpha & 1-2\\alpha & \\alpha&&&\\\\\n",
"& \\alpha & 1-2\\alpha & \\alpha&&\\\\\n",
"&&& \\ddots &\\\\\n",
"&&& \\alpha & 1-2\\alpha & \\alpha \\\\\n",
"&&&& \\alpha-\\frac{1}{3}\\alpha & 1-2\\alpha+\\frac{4}{3}\\alpha\n",
"\\end{pmatrix}U^{n} + \\begin{pmatrix}\n",
"-\\frac{2\\alpha\\Delta x}{3}(\\beta_L^\\prime)^{n}\\\\\n",
"0\\\\\n",
"0\\\\\n",
"\\vdots\\\\\n",
"0\\\\\n",
"\\frac{2\\alpha\\Delta x}{3}(\\beta_R^\\prime)^{n}\\\\\n",
"\\end{pmatrix}.\n",
"\\end{align*}\n",
"\n",
"### Nonlinear Term\n",
"\n",
"When there is a nonlinear term in the equation, say the PDE is $u_t = Lu + N(u)$, then the corresponding Crank-Nicolson method is \n",
"\\begin{align*}\n",
"U^{n+1} &= U^{n} + \\int_{t_n}^{t_{n+1}} \\dot U\\,dt\\\\\n",
"&= U^{n} + \\int_{t_n}^{t_{n+1}} L_h U + N(U)\\,dt\\\\\n",
"&\\approx U^{n} + \\Delta t\\left(\\frac{L_h U^{n} + N(U^{n}) + L_h U^{n+1} + N(U^{n+1})}{2}\\right). \n",
"\\end{align*}\n",
"Thus\n",
"\\begin{align*}\n",
"\\left(I - \\frac{\\Delta t}{2}L_h \\right)U^{n+1} - \\frac{\\Delta t}{2} N(U^{n+1}) = \\left(I + \\frac{\\Delta t}{2}L_h \\right)U^{n} + \\frac{\\Delta t}{2} N(U^{n}). \n",
"\\end{align*}\n",
"Now on the right hand side apply the 2nd order approximation (in $t$)\n",
"\\begin{align*}\n",
"N(U^{n+1}) = N(U^{n}) + (U^{n+1} - U^{n})\\circ N^\\prime(U^n), \n",
"\\end{align*}\n",
"where $\\circ$ denotes the element-wise product. \n",
"We obtain \n",
"\\begin{align*}\n",
"&\\left(I - \\frac{\\Delta t}{2}L_h \\right)U^{n+1} - \\frac{\\Delta t}{2} \\left(N(U^{n}) + (U^{n+1} - U^{n})\\circ N^\\prime(U^n)\\right) = \\left(I + \\frac{\\Delta t}{2}L_h \\right)U^{n} + \\frac{\\Delta t}{2} N(U^{n})\\\\\n",
"\\implies\\quad&\\left(I - \\frac{\\Delta t}{2}L_h \\right)U^{n+1} - \\frac{\\Delta t}{2} N^\\prime(U^n)\\circ U^{n+1} = \\left(I + \\frac{\\Delta t}{2}L_h\\right)U^{n} - \\frac{\\Delta t}{2} N^\\prime(U^n)\\circ U^{n} + \\Delta t N(U^n)\\\\\n",
"\\implies\\quad&\\left(I - \\frac{\\Delta t}{2}L_h - \\frac{\\Delta t}{2} N^\\prime(U^n)\\right)U^{n+1} = \\left(I + \\frac{\\Delta t}{2}L_h - \\frac{\\Delta t}{2} N^\\prime(U^n)\\right)U^{n} + \\Delta t N(U^n). \n",
"\\end{align*}\n",
"In most cases this requires reconstruction of the matrix at each time step. \n",
"Note that the above argument is not specific to finite difference method. Same argument works for spectral methods too. \n",
"\n",
"\n",
"### Alternating-Direction Implicit Method (No Cross Term)\n",
"\n",
"To solve $u_t = u_{xx} + u_{yy}$, first discretize $x$ and $y$ to get the system of ODEs \n",
"\\begin{align*}\n",
"\\dot U = (\\delta_x^+\\delta_x^- + \\delta_y^+\\delta_y^-) U, \n",
"\\end{align*}\n",
"where $U$ is a matrix function of $t$: \n",
"\\begin{align*}\n",
"U(t) = \\begin{pmatrix}\n",
"U_{11}(t) & U_{12}(t) & \\cdots & U_{1m_2}(t)\\\\\n",
"U_{21}(t) & U_{22}(t) & \\cdots & U_{2m_2}(t)\\\\\n",
"\\vdots & \\vdots & \\ddots & \\vdots\\\\\n",
"U_{m_11}(t) & U_{m_22}(t) & \\cdots & U_{m_1m_2}(t)\\\\\n",
"\\end{pmatrix}. \n",
"\\end{align*}\n",
"The Crank-Nicolson method is\n",
"\\begin{align*}\n",
"\\left(I - \\frac{\\Delta t}{2}(\\delta_x^+\\delta_x^- + \\delta_y^+\\delta_y^-) \\right)U^{n+1} = \\left(I + \\frac{\\Delta t}{2}(\\delta_x^+\\delta_x^- + \\delta_y^+\\delta_y^-)\\right)U^{n}, \n",
"\\end{align*}\n",
"where the operators can be approximated by\n",
"\\begin{align*}\n",
"\\left(I - \\frac{\\Delta t}{2}(\\delta_x^+\\delta_x^- + \\delta_y^+\\delta_y^-) \\right) &\\approx \\left(I - \\frac{\\Delta t}{2}\\delta_x^+\\delta_x^- \\right)\\left(I - \\frac{\\Delta t}{2}\\delta_y^+\\delta_y^- \\right), \\\\\n",
"\\left(I + \\frac{\\Delta t}{2}(\\delta_x^+\\delta_x^- + \\delta_y^+\\delta_y^-)\\right) &\\approx \\left(I + \\frac{\\Delta t}{2}\\delta_x^+\\delta_x^- \\right)\\left(I + \\frac{\\Delta t}{2}\\delta_y^+\\delta_y^- \\right). \n",
"\\end{align*}\n",
"These approximations have truncation error smaller than $O(\\Delta t^2 + \\Delta x^2 + \\Delta y^2)$. \n",
"The system of equations we need to solve can be written as \n",
"\\begin{align*}\n",
"\\left(I - \\frac{\\Delta t}{2}\\delta_x^+\\delta_x^- \\right)\\left(I - \\frac{\\Delta t}{2}\\delta_y^+\\delta_y^- \\right)U^{n+1} = \\left(I + \\frac{\\Delta t}{2}\\delta_x^+\\delta_x^- \\right)\\left(I + \\frac{\\Delta t}{2}\\delta_y^+\\delta_y^- \\right)U^{n}, \n",
"\\end{align*}\n",
"which can be split into \n",
"\\begin{align*}\n",
"&\\left(I - \\frac{\\Delta t}{2}\\delta_x^+\\delta_x^- \\right)V^{n} = \\left(I + \\frac{\\Delta t}{2}\\delta_y^+\\delta_y^- \\right)U^{n}, \\\\\n",
"&\\left(I - \\frac{\\Delta t}{2}\\delta_y^+\\delta_y^- \\right)U^{n+1} = \\left(I + \\frac{\\Delta t}{2}\\delta_x^+\\delta_x^- \\right)V^{n}, \n",
"\\end{align*}\n",
"where $V^n$ is an intermediate value. This is the Peaceman-Rachford scheme. Although not immediately obvious (as the operators are not commutative), it can be shown that thus obtained solution will be equivalent to the above. Now finding $U^{n+1}$ is a two-step process. Since each step involves operator in only one direction, one can simply apply the Thomas algorithm $m_j$ times, avoiding iterative solvers. \n",
"\n",
"### ADI With Cross Term\n",
"\n",
"To solve $u_t = u_{xx} + u_{yy} + u_{xy}$, first discretize in space to get the system of ODEs\n",
"$$\n",
"\\dot U = (\\delta_x^+\\delta_x^- + \\delta_y^+\\delta_y^- + \\delta_x^0\\delta_y^0) U, \n",
"$$\n",
"where\n",
"\\begin{align*}\n",
"\\delta_x^0 &= \\frac{S_x-S_x^{-1}}{2\\Delta x}, \\\\\n",
"\\delta_y^0 &= \\frac{S_y-S_y^{-1}}{2\\Delta y}. \n",
"\\end{align*}\n",
"No way to split the operator into separate ones in different directions, we treat the mixed derivative fully explicitly and write\n",
"\\begin{align*}\n",
"\\left(I - \\frac{\\Delta t}{2}(\\delta_x^+\\delta_x^- + \\delta_y^+\\delta_y^-) \\right)U^{n+1} = \\left(I + \\frac{\\Delta t}{2}(\\delta_x^+\\delta_x^- + \\delta_y^+\\delta_y^-) + \\Delta t\\delta_x^0\\delta_y^0\\right)U^{n}. \n",
"\\end{align*}\n",
"By splitting the opertor only on the left, we obtain \n",
"\\begin{align*}\n",
"\\left(I - \\frac{\\Delta t}{2}\\delta_x^+\\delta_x^- \\right)\\left(I - \\frac{\\Delta t}{2}\\delta_y^+\\delta_y^- \\right)U^{n+1} = \\left(I + \\frac{\\Delta t}{2}(\\delta_x^+\\delta_x^- + \\delta_y^+\\delta_y^-) + \\Delta t\\delta_x^0\\delta_y^0\\right)U^{n}. \n",
"\\end{align*}\n",
"A variant of this scheme is the Douglas-Rachford scheme\n",
"\\begin{align*}\n",
"&\\left(I - \\frac{\\Delta t}{2}\\delta_x^+\\delta_x^- \\right)V^n = \\left(I + \\frac{\\Delta t}{2}\\delta_x^+\\delta_x^- + \\Delta t\\delta_y^+\\delta_y^- + \\Delta t\\delta_x^0\\delta_y^0\\right)U^{n}, \\\\\n",
"&\\left(I - \\frac{\\Delta t}{2}\\delta_y^+\\delta_y^- \\right)U^{n+1} = V^n - \\frac{\\Delta t}{2}\\delta_y^+\\delta_y^- U^n. \n",
"\\end{align*}\n",
"Since the cross term is not integrated with the trapezoidal rule, the convergence order in time of the Douglas-Rachford reduces to 1. To recover 2nd order convergence in time, apply a predictor-corrector scheme: \n",
"\n",
"* Predictor: \n",
"\\begin{align*}\n",
"\\left(I - \\frac{\\Delta t}{2}\\delta_x^+\\delta_x^- \\right)V_1^{n} &= \\left(I + \\frac{\\Delta t}{2}\\delta_x^+\\delta_x^- + \\Delta t\\delta_y^+\\delta_y^- + \\Delta t\\delta_x^0\\delta_y^0\\right)U^{n}, \\\\\n",
"\\left(I - \\frac{\\Delta t}{2}\\delta_y^+\\delta_y^- \\right)V_2^{n} &= V_1^n - \\frac{\\Delta t}{2}\\delta_y^+\\delta_y^- U^n. \n",
"\\end{align*}\n",
"\n",
"* Corrector: \n",
"\\begin{align*}\n",
"&\\left(I - \\frac{\\Delta t}{2}\\delta_x^+\\delta_x^- \\right)V_3^{n} = \\left(I + \\frac{\\Delta t}{2}\\delta_x^+\\delta_x^- + \\Delta t \\delta_y^+\\delta_y^- + \\frac{\\Delta t}{2}\\delta_x^0\\delta_y^0\\right)U^{n} + \\frac{\\Delta t}{2} V_2^{n}, \\\\\n",
"&\\left(I - \\frac{\\Delta t}{2}\\delta_y^+\\delta_y^- \\right)U^{n+1} = V_3^n - \\frac{\\Delta t}{2}\\delta_y^+\\delta_y^- U^n. \n",
"\\end{align*}\n",
"\n",
"\n",
"Here the Douglas-Rachford scheme is first applied and the result $V_2^{n}$ is a predictor of the solution at time $t_{n+1}$, to be corrected. This scheme has 2nd order convergence both in time and space. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Spectral Method Summary (Nodal Continuous Galerkin)\n",
"\n",
"To solve $u_t = u_{xx}$ in $x\\in[-1, 1]$, first write it in weak form and integrate by part in $x$: \n",
"\\begin{align*}\n",
"\\int_{-1}^1 u_t(t, x)\\phi(x)\\,dx &= \\int_{-1}^{1} u_{xx}(t, x)\\phi(x)\\,dx\\\\\n",
"&= u_{x}\\phi(x)\\Big|_{-1}^1 - \\int_{-1}^1 u_{x}\\phi^\\prime(x)\\,dx\\\\\n",
"&= - \\int_{-1}^1 u_{x}\\phi^\\prime(x)\\,dx, \n",
"\\end{align*}\n",
"where the last equation holds as long as we choose a function $\\phi(x)$ that vanishes at the boundaries. \n",
"Replacing all functions by Lagrange interpolating polynomials and integrals by Gauss-Lobatto quadratures \n",
"\\begin{align*}\n",
"u(t, x) &\\approx \\sum_{p=0}^N U_p(t)\\ell_p(x),\\\\\n",
"\\phi(x) &\\approx \\sum_{i=0}^N \\phi_i\\ell_i(x),\\\\\n",
"\\int_{-1}^1 f(x)\\,dx &\\approx \\sum_{n=0}^N f(x_n)w_n,\n",
"\\end{align*}\n",
"we get\n",
"\\begin{align*}\n",
"\\sum_{p=0}^N\\sum_{i=0}^N\\sum_{n=0}^N \\dot U_p(t)\\ell_p(x_n)\\phi_i\\ell_i(x_n)w_n =& -\\sum_{p=0}^N\\sum_{i=0}^N\\sum_{n=0}^N U_p(t)\\ell_p^\\prime(x_n)\\phi_i\\ell_i^\\prime(x_n)w_n. \n",
"\\end{align*}\n",
"Note that doing so we have introduced an error in space and the error is of order $O\\left(C^{-N}\\right)$ as the functions are smooth. \n",
"There are many zero terms in the above equation as the Lagrange basis polynomials have the property\n",
"\\begin{align*}\n",
"\\ell_p(x_n) =\n",
"\\begin{cases}\n",
"1 &\\mbox{ if }p=n\\\\\n",
"0 &\\mbox{ if }p\\ne n\n",
"\\end{cases}.\n",
"\\end{align*}\n",
"Get rid of the zero terms to get \n",
"\\begin{align*}\n",
"\\sum_{i=0}^N \\dot U_i(t)\\phi_iw_i =& -\\sum_{p=0}^N\\sum_{i=0}^N\\sum_{n=0}^N U_p(t)\\ell_p^\\prime(x_n)\\phi_i\\ell_i^\\prime(x_n)w_n\\\\\n",
"=& -\\sum_{p=0}^N\\sum_{i=0}^N\\sum_{n=0}^N U_p(t)D_{np}\\phi_iD_{ni}w_n, \n",
"\\end{align*}\n",
"where $D_{ij} = \\ell_j^\\prime(x_i)$. \n",
"As this holds for arbitrary $\\phi_i$'s, the coefficients of $\\phi_i$ on both sides should equal: \n",
"\\begin{align*}\n",
"\\dot U_i(t) w_i =& -\\sum_{p=0}^N\\sum_{n=0}^N U_p(t)D_{np}D_{ni}w_n\\qquad \\forall i=1, 2,\n",
"\\ldots, N-1. \n",
"\\end{align*}\n",
"Next we write the equation in matrix form using the identity $(Ax)_i = \\sum_j A_{ij}x_j$. We obtain \n",
"\\begin{align*}\n",
"W\\dot{\\vec{U}} = -D^TWD\\vec U, \n",
"\\end{align*}\n",
"where\n",
"\\begin{align*}\n",
"&\\vec U = \n",
"\\begin{pmatrix}\n",
"U_0(t)\\\\\n",
"U_1(t)\\\\\n",
"U_2(t)\\\\\n",
"\\vdots\\\\\n",
"U_N(t)\n",
"\\end{pmatrix}, W = \n",
"\\begin{pmatrix}\n",
"w_0 & & & & 0 \\\\\n",
"& w_1 & & & \\\\\n",
"& & w_2 & & \\\\\n",
"& & & \\ddots & \\\\\n",
"0 & & & & w_N\n",
"\\end{pmatrix}, \\\\\\\\\n",
"&D = \n",
"\\begin{pmatrix}\n",
"D_{00} & D_{01} & D_{02} & \\cdots & D_{0N} \\\\\n",
"D_{10} & D_{11} & D_{12} & \\cdots & D_{1N} \\\\\n",
"D_{20} & D_{21} & D_{22} & \\cdots & D_{2N} \\\\\n",
"\\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
"D_{N0} & D_{N1} & D_{N2} & \\cdots & D_{NN}\n",
"\\end{pmatrix} = \\begin{pmatrix}\n",
"\\ell_0^\\prime(\\xi_0) & \\ell_1^\\prime(\\xi_0) & \\ell_2^\\prime(\\xi_0) & \\cdots & \\ell_N^\\prime(\\xi_0) \\\\\n",
"\\ell_0^\\prime(\\xi_1) & \\ell_1^\\prime(\\xi_1) & \\ell_2^\\prime(\\xi_1) & \\cdots & \\ell_N^\\prime(\\xi_1) \\\\\n",
"\\ell_0^\\prime(\\xi_2) & \\ell_1^\\prime(\\xi_2) & \\ell_2^\\prime(\\xi_2) & \\cdots & \\ell_N^\\prime(\\xi_2) \\\\\n",
"\\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
"\\ell_0^\\prime(\\xi_N) & \\ell_1^\\prime(\\xi_N) & \\ell_2^\\prime(\\xi_N) & \\cdots & \\ell_N^\\prime(\\xi_N)\n",
"\\end{pmatrix}, \n",
"\\end{align*}\n",
"or equivalently\n",
"\\begin{align*}\n",
"\\dot{\\vec{U}} = -W^{-1}D^TWD\\vec U = L\\vec U, \n",
"\\end{align*}\n",
"where $L = -W^{-1}D^TWD$. \n",
"Now discretize $t$ and apply the trapezoidal rule,\n",
"\\begin{align*}\n",
"\\vec U(t_{n+1}) &= \\vec U(t_n) + \\int_{t_n}^{t_{n+1}} \\dot{\\vec U}(t)\\,dt\\\\\n",
"&= \\vec U(t_n) + \\int_{t_n}^{t_{n+1}} L\\vec U(t)\\,dt\\\\\n",
"&\\approx \\vec U(t_n) + \\frac{\\Delta t L}{2}\\left(\\vec U(t_n) + \\vec U(t_{n+1}) \\right).\n",
"\\end{align*}\n",
"Rearrange the terms to obtain the linear system we need to solve: \n",
"\\begin{align*}\n",
"\\left(I - \\frac{\\Delta t}{2}L\\right)\\vec U(t_{n+1}) =\n",
"\\left(I + \\frac{\\Delta t}{2}L\\right)\\vec U(t_n).\n",
"\\end{align*}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Spectral Method in Multidimensions (Nodal Continuous Galerkin)\n",
"\n",
"The goal here is to solve $u_t = \\mathcal L u(\\vec x)$ on $\\vec x\\in\\Omega = [-1, 1]\\times[-1, 1]$. For simplicity, think of $\\mathcal L$ as for example the 2D Laplacian $\\mathcal L u(x, y) = u_{xx} + u_{yy}$, but the below argument holds for any operators in the form\n",
"\\begin{align*}\n",
"\\mathcal L = a\\frac{\\partial^2}{\\partial x^2} + b\\frac{\\partial^2}{\\partial x\\partial y} + c\\frac{\\partial^2}{\\partial y^2} + d\\frac{\\partial}{\\partial x} + e\\frac{\\partial}{\\partial y} + f. \n",
"\\end{align*}\n",
"First write the PDE in weak form: \n",
"\\begin{align*}\n",
"\\int_\\Omega u_t(t, \\vec x)\\phi(\\vec x)\\,d\\vec x &= \\int_\\Omega \\mathcal L u(t, \\vec x)\\phi(\\vec x)\\,d\\vec x. \n",
"\\end{align*}\n",
"Now replace all functions by Lagrange interpolating polynomials and integrals by Gauss-Lobatto quadratures \n",
"\\begin{align*}\n",
"u(t, \\vec x) &\\approx \\sum_{\\mathbf p} U_\\mathbf p(t)\\ell_\\mathbf p(\\vec x),\\\\\n",
"\\phi(\\vec x) &\\approx \\sum_{\\mathbf i} \\phi_\\mathbf i\\ell_\\mathbf i(\\vec x),\\\\\n",
"\\int_\\Omega f(\\vec x)\\,d\\vec x &\\approx \\sum_{\\mathbf n} f(\\vec x_\\mathbf n)w_\\mathbf n,\n",
"\\end{align*}\n",
"where the boldface indices are in multidimensions, for example $\\mathbf i$ can be $(0, 0), (0, 1), (1, 1), \\ldots$. \n",
"We obtain\n",
"\\begin{align*}\n",
"\\sum_{\\mathbf p}\\sum_{\\mathbf i}\\sum_{\\mathbf n} \\dot U_\\mathbf p(t)\\ell_\\mathbf p(\\vec x_\\mathbf n)\\phi_\\mathbf i\\ell_\\mathbf i(\\vec x_\\mathbf n)w_\\mathbf n =& \\sum_{\\mathbf p}\\sum_{\\mathbf i}\\sum_{\\mathbf n} U_\\mathbf p(t)\\mathcal L\\ell_\\mathbf p(\\vec x_\\mathbf n)\\phi_\\mathbf i\\ell_\\mathbf i(\\vec x_\\mathbf n)w_\\mathbf n. \n",
"\\end{align*}\n",
"Note that doing so we have introduced an error in space and the error is of order $O\\left(C^{-N}\\right)$ if the functions are analytic, where $N$ is total number of nodes in one dimension. \n",
"There are many zero terms in the above equation as the Lagrange basis polynomials have the property\n",
"\\begin{align*}\n",
"\\ell_\\mathbf p(\\vec x_\\mathbf n) =\n",
"\\begin{cases}\n",
"1 &\\mbox{ if }\\mathbf p=\\mathbf n\\\\\n",
"0 &\\mbox{ if }\\mathbf p\\ne \\mathbf n\n",
"\\end{cases}.\n",
"\\end{align*}\n",
"Get rid of the zero terms to get \n",
"\\begin{align*}\n",
"\\sum_{\\mathbf i} \\dot U_\\mathbf i(t)\\phi_\\mathbf iw_\\mathbf i \n",
"=& \\sum_{\\mathbf p}\\sum_{\\mathbf i} U_\\mathbf p(t)\\mathcal L\\ell_\\mathbf p(\\vec x_\\mathbf i)\\phi_\\mathbf iw_\\mathbf i\\\\\n",
"=& \\sum_{\\mathbf p}\\sum_{\\mathbf i} U_\\mathbf p(t)L_{\\mathbf i\\mathbf p}\\phi_\\mathbf iw_\\mathbf i, \n",
"\\end{align*}\n",
"where $L_{\\mathbf i\\mathbf j} = \\mathcal L\\ell_\\mathbf j(\\vec x_\\mathbf i)$. \n",
"As this holds for arbitrary $\\phi_\\mathbf i$'s, the coefficients of $\\phi_\\mathbf iw_\\mathbf i$ on both sides should equal: \n",
"\\begin{align*}\n",
"\\dot U_\\mathbf i(t) =& \\sum_{\\mathbf p}U_\\mathbf p(t)L_{\\mathbf i\\mathbf p}\\qquad \\forall \\mathbf i. \n",
"\\end{align*}\n",
"Next we write the equation in matrix form using the identity $(Ax)_\\mathbf i = \\sum_{\\mathbf j} A_{\\mathbf i\\mathbf j}x_\\mathbf j$. We obtain a clean\n",
"\\begin{align*}\n",
"\\dot{\\vec{U}} = L\\vec U. \n",
"\\end{align*}\n",
"\n",
"Now discretize $t$ and apply the trapezoidal rule,\n",
"\\begin{align*}\n",
"\\vec U(t_{n+1}) &= \\vec U(t_n) + \\int_{t_\\mathbf n}^{t_{n+1}} \\dot{\\vec U}(t)\\,dt\\\\\n",
"&= \\vec U(t_n) + \\int_{t_n}^{t_{n+1}} L\\vec U(t)\\,dt\\\\\n",
"&\\approx \\vec U(t_n) + \\frac{\\Delta t L}{2}\\left(\\vec U(t_n) + \\vec U(t_{n+1}) \\right).\n",
"\\end{align*}\n",
"Rearrange the terms to obtain the linear system we need to solve: \n",
"\\begin{align*}\n",
"\\left(I - \\frac{\\Delta t}{2}L\\right)\\vec U(t_{n+1}) =\n",
"\\left(I + \\frac{\\Delta t}{2}L\\right)\\vec U(t_n).\n",
"\\end{align*}\n",
"\n",
"### Dirichlet Boundary Conditions\n",
"\n",
"Given Dirichlet boundary conditions, in the system of ODEs \n",
"\\begin{align*}\n",
"\\dot U_\\mathbf i(t) =& \\sum_{\\mathbf p}U_\\mathbf p(t)L_{\\mathbf i\\mathbf p}\\qquad \\forall \\mathbf i, \n",
"\\end{align*}\n",
"$U_\\mathbf i(t)$ values at the boundary are known so no need to solve ODEs for those. \n",
"Denote by $\\mathbf i_{\\text{int}}$ and $\\mathbf p_{\\text{int}}$ indices for interior points of $\\Omega$, and denote by $\\mathbf i_{\\text{bdy}}$ and $\\mathbf p_{\\text{bdy}}$ indices for points at the boundary. We are only interested in \n",
"\\begin{align*}\n",
"\\dot U_{\\mathbf i_{\\text{int}}}(t) =& \\sum_{\\mathbf p}U_\\mathbf p(t)L_{\\mathbf i_{\\text{int}}\\mathbf p}\\\\\n",
"=& \\sum_{\\mathbf p_{\\text{int}}}U_{\\mathbf p_{\\text{int}}}(t)L_{\\mathbf i_{\\text{int}}\\mathbf p_{\\text{int}}} + \\sum_{\\mathbf p_{\\text{bdy}}}U_{\\mathbf p_{\\text{bdy}}}(t)L_{\\mathbf i_{\\text{int}}\\mathbf p_{\\text{bdy}}}, \n",
"\\end{align*}\n",
"where the 2nd term in the last equation is simply a known scalar as the $U_{\\mathbf p_{\\text{bdy}}}(t)$ values are given. \n",
"In matrix form, the system of ODEs is\n",
"\\begin{align*}\n",
"\\dot{\\vec U}_{\\text{int}} =& L_{\\text{int}} \\vec U_{\\text{int}} + L_{\\text{int, bdy}} \\vec U_{\\text{bdy}}, \n",
"\\end{align*}\n",
"where \n",
"${\\vec U}_{\\text{int}}$ is a vector of all $U_{\\mathbf i_{\\text{int}}}(t)$ at interior points, \n",
"${\\vec U}_{\\text{bdy}}$ a vector of all $U_{\\mathbf i_{\\text{bdy}}}(t)$ at the boundary, \n",
"and $L_{\\text{int, bdy}}$ is a (tall) matrix of all $L_{{\\mathbf i_{\\text{int}}}{\\mathbf p_{\\text{bdy}}}} = \\mathcal L\\ell_{\\mathbf p_{\\text{bdy}}}(\\vec x_{\\mathbf i_{\\text{int}}})$, that's the Laplacian of all boundary Lagrange basis functions at interior points. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sparse Grid"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"tags": []
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | q | \n",
" 2 | \n",
" 3 | \n",
" 4 | \n",
" 5 | \n",
" 6 | \n",
" 7 | \n",
" 8 | \n",
" 9 | \n",
" 10 | \n",
"
\n",
" \n",
" \n",
" \n",
" | d | \n",
" 2 | \n",
" 2 | \n",
" 2 | \n",
" 2 | \n",
" 2 | \n",
" 2 | \n",
" 2 | \n",
" 2 | \n",
" 2 | \n",
"
\n",
" \n",
" | number of grid points on one boundary | \n",
" 3 | \n",
" 5 | \n",
" 9 | \n",
" 17 | \n",
" 33 | \n",
" 65 | \n",
" 129 | \n",
" 257 | \n",
" 513 | \n",
"
\n",
" \n",
" | number of full grid points | \n",
" 9 | \n",
" 25 | \n",
" 81 | \n",
" 289 | \n",
" 1089 | \n",
" 4225 | \n",
" 16641 | \n",
" 66049 | \n",
" 263169 | \n",
"
\n",
" \n",
" | number of sparse grid points | \n",
" 9 | \n",
" 21 | \n",
" 49 | \n",
" 113 | \n",
" 257 | \n",
" 577 | \n",
" 1281 | \n",
" 2817 | \n",
" 6145 | \n",
"
\n",
" \n",
" | number of full grid interior points | \n",
" 1 | \n",
" 9 | \n",
" 49 | \n",
" 225 | \n",
" 961 | \n",
" 3969 | \n",
" 16129 | \n",
" 65025 | \n",
" 261121 | \n",
"
\n",
" \n",
" | number of sparse grid interior points | \n",
" 1 | \n",
" 5 | \n",
" 17 | \n",
" 49 | \n",
" 129 | \n",
" 321 | \n",
" 769 | \n",
" 1793 | \n",
" 4097 | \n",
"
\n",
" \n",
" | ratio number of interior points (sparse/full) x 100 (in %) | \n",
" 100 | \n",
" 55 | \n",
" 34 | \n",
" 21 | \n",
" 13 | \n",
" 8 | \n",
" 4 | \n",
" 2 | \n",
" 1 | \n",
"
\n",
" \n",
"
\n",
"