{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Thomas tri-diagonal method for 1-D BVP"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider solving\n",
"$$\n",
"-u_{xx} = f(x), \\qquad x \\in [a,b]\n",
"$$\n",
"with boundary condition\n",
"$$\n",
"u(a) = \\alpha, \\qquad u(b) = \\beta\n",
"$$\n",
"\n",
"Make a partition of $n$ intervals with uniform spacing and grid points\n",
"$$\n",
"h = \\frac{b-a}{n}, \\qquad x_i = a + i h, \\qquad i=0,1,\\ldots,n\n",
"$$\n",
"The finite difference scheme is\n",
"\\begin{eqnarray*}\n",
"u_0 &=& \\alpha \\\\\n",
"-\\frac{u_{i-1} - 2 u_i + u_{i+1}}{h^2} &=& f_i, \\qquad i=1,2,\\ldots,n-1 \\\\\n",
"u_n &=& \\beta\n",
"\\end{eqnarray*}\n",
"We can rewrite this as\n",
"\\begin{eqnarray*}\n",
"u_0 &=& \\alpha \\\\\n",
"-u_{i-1} + 2 u_i - u_{i+1} &=& h^2 f_i, \\qquad i=1,2,\\ldots,n-1 \\\\\n",
"u_n &=& \\beta\n",
"\\end{eqnarray*}\n",
"The above set of equations can be written as a matrix equation where the matrix is tridiagonal. \n",
"$$\n",
"\\left[\\begin{array}{cccccc}\n",
"1 & 0 & 0 & 0 & \\ldots & 0\\\\\n",
"-1 & 2 & -1 & 0 & \\ldots & 0 \\\\\n",
"0 & -1 & 2 & -1 & \\ldots & 0\\\\\n",
"\\vdots & & & & & \\vdots \\\\\n",
"0 & \\ldots & 0 & -1 & 2 & -1 \\\\\n",
"0 & \\ldots & \\ldots & 0 & 0 & 1\n",
"\\end{array}\\right] \\begin{bmatrix}\n",
"u_0 \\\\ u_1 \\\\ u_2 \\\\ \\vdots \\\\ u_{n-1} \\\\ u_n \\end{bmatrix} =\n",
"\\begin{bmatrix}\n",
"\\alpha \\\\ h^2 f_1 \\\\ h^2 f_2 \\\\ \\vdots \\\\ h^2 f_{n-1} \\\\ \\beta \\end{bmatrix}\n",
"$$\n",
"We first implement Thomas algorithm but before that let us import some Python stuff."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"import numpy as np\n",
"from matplotlib import pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Thomas Algorithm\n",
"\n",
"Consider a tridiagonal matrix problem\n",
"\n",
"$$\n",
"A x = f\n",
"$$\n",
"\n",
"The general form of a $N \\times N$ tridiagonal matrix is\n",
"\n",
"$$\n",
"A = \\left[\\begin{array}{cccccc}\n",
"a_0 & c_0 & 0 & 0 & \\ldots & 0\\\\\n",
"b_1 & a_1 & c_1 & 0 & \\ldots & 0 \\\\\n",
"0 & b_2 & a_2 & c_2 & \\ldots & 0\\\\\n",
"\\vdots & & & & & \\vdots \\\\\n",
"0 & \\ldots & 0 & b_{N-2} & a_{N-2} & c_{N-2} \\\\\n",
"0 & \\ldots & \\ldots & 0 & b_{N-1} & a_{N-1}\n",
"\\end{array}\\right]\n",
"$$\n",
"\n",
"We write this as the following LU decomposition\n",
"\n",
"$$\n",
"A = \\left[\\begin{array}{cccccc}\n",
"\\alpha_0 & 0 & 0 & 0 & \\ldots & 0\\\\\n",
"b_1 & \\alpha_1 & 0 & 0 & \\ldots & 0 \\\\\n",
"0 & b_2 & \\alpha_2 & 0 & \\ldots & 0\\\\\n",
"\\vdots & & & & & \\vdots \\\\\n",
"0 & \\ldots & 0 & b_{N-2} & \\alpha_{N-2} & 0 \\\\\n",
"0 & \\ldots & \\ldots & 0 & b_{N-1} & \\alpha_{N-1}\n",
"\\end{array}\\right]\n",
"\\left[\\begin{array}{cccccc}\n",
"1 & \\gamma_0 & 0 & 0 & \\ldots & 0\\\\\n",
"0 & 1 & \\gamma_1 & 0 & \\ldots & 0 \\\\\n",
"0 & 0 & 1 & \\gamma_2 & \\ldots & 0\\\\\n",
"\\vdots & & & & & \\vdots \\\\\n",
"0 & \\ldots & 0 & 0 & 1 & \\gamma_{N-2} \\\\\n",
"0 & \\ldots & \\ldots & 0 & 0 & 1\n",
"\\end{array}\\right]\n",
"$$\n",
"\n",
"By matching the lhs with rhs, we can determine the $\\alpha_i$ and $\\gamma_i$\n",
"\n",
"$$\n",
"\\alpha_0 = a_0, \\qquad \\gamma_0 = c_0/\\alpha_0\n",
"$$\n",
"\n",
"$$\n",
"\\alpha_i = a_i - b_i \\gamma_{i-1}, \\quad \\gamma_i = c_i/\\alpha_i \\qquad i=1,2,\\ldots,N-2\n",
"$$\n",
"\n",
"$$\n",
"\\alpha_{N-1} = a_{N-1} - b_{N-1} \\gamma_{N-2}\n",
"$$\n",
"\n",
"Actually there is no need to have separate arrays for $\\alpha$, $\\beta$. We will modify the $a$ and overwrite it with $\\alpha$ and $c$ array will contain $\\gamma$, while the $b$ array is unchanged. The LU decomposition is given by\n",
"\n",
"$$\n",
"c_0 = c_0/a_0\n",
"$$\n",
"\n",
"$$\n",
"a_i = a_i - b_i c_{i-1}, \\quad c_i = c_i/a_i \\qquad i=1,2,\\ldots,N-2\n",
"$$\n",
"\n",
"$$\n",
"a_{N-1} = a_{N-1} - b_{N-1} c_{N-2}\n",
"$$\n",
"\n",
"Now we solve $Lx = f$\n",
"\n",
"$$\n",
"x_0 = f_0/a_0\n",
"$$\n",
"\n",
"$$\n",
"x_i = \\frac{1}{a_i}(f_i - b_i x_{i-1}), \\qquad i=1,2,\\ldots,N-1\n",
"$$\n",
"\n",
"Now we solve $Ux = x$, i.e., we store the result in the same array $x$. \n",
"\n",
"$$\n",
"x_i = x_i - c_i x_{i+1}, \\qquad i=N-2,N-3,\\ldots,0\n",
"$$\n",
"\n",
"We implement two functions; the first one computes the LU decomposition of the above matrix and the second one solves the matrix problem. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# construct LU decomposition\n",
"def tdma1(a,b,c):\n",
" n = len(a)\n",
" c[0] = c[0]/a[0]\n",
" for i in range(1,n):\n",
" a[i] = a[i] - b[i]*c[i-1]\n",
" c[i] = c[i]/a[i]\n",
" a[n-1] = a[n-1] - b[n-1]*c[n-2]\n",
" return a,b,c\n",
"\n",
"# solve\n",
"def tdma2(a,b,c,f):\n",
" n = len(f)\n",
" x = np.empty_like(f)\n",
" # solve L y = f\n",
" x[0] = f[0]/a[0]\n",
" for i in range(1,n):\n",
" x[i] = (f[i] - b[i]*x[i-1])/a[i]\n",
" # solve U x = y\n",
" for i in range(n-2,-1,-1):\n",
" x[i] = x[i] - c[i]*x[i+1]\n",
" return x"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solution of ODE"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Choose\n",
"$$\n",
"f(x) = \\sin(2\\pi x), \\qquad x \\in [a,b] = [0,1]\n",
"$$\n",
"for which the exact solution is\n",
"$$\n",
"u(x) = \\frac{1}{(2\\pi)^2}\\sin(2\\pi x)\n",
"$$\n",
"and $\\alpha = \\beta = 0$. We now solve the problem by calling the above two functions."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Max error norm = 0.00013349155840115118\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"xmin, xmax = 0.0, 1.0\n",
"uexact = lambda x: np.sin(2*np.pi*x)/(2*np.pi)**2\n",
"frhs = lambda x: np.sin(2*np.pi*x)\n",
"\n",
"n = 25 # Number of intervals, so we have n+1 points\n",
"h = (xmax - xmin)/n\n",
"\n",
"x = np.linspace(xmin, xmax, n+1) # Grid\n",
"f = h**2 * frhs(x) # Right hand side\n",
"f[0] = uexact(xmin); f[n] = uexact(xmax);\n",
"\n",
"# Create the three diagonals\n",
"a = 2.0*np.ones(n+1)\n",
"a[0] = 1.0; a[n] = 1.0\n",
"\n",
"b = -1.0*np.ones(n+1)\n",
"b[n] = 0.0\n",
"\n",
"c = -1.0*np.ones(n+1)\n",
"c[0] = 0.0\n",
"\n",
"# Compute LU decomposition and solve the problem\n",
"a,b,c, = tdma1(a,b,c)\n",
"u = tdma2(a,b,c,f)\n",
"\n",
"ue = uexact(x) # Exact solution\n",
"\n",
"print('Max error norm = ',np.abs(u-ue).max())\n",
"\n",
"plt.plot(x,ue,x,u,'o')\n",
"plt.xlabel('x'); plt.ylabel('u')\n",
"plt.grid(True)\n",
"plt.legend((\"Exact\",\"Numerical\"));"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
}
},
"nbformat": 4,
"nbformat_minor": 4
}