{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Classical iterative methods"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider solving\n",
"$$\n",
"-u_{xx} = f(x), \\qquad x \\in [0,1]\n",
"$$\n",
"with boundary condition\n",
"$$\n",
"u(0) = u(1) = 0\n",
"$$\n",
"Choose\n",
"$$\n",
"f(x) = 1\n",
"$$\n",
"The exact solution is\n",
"$$\n",
"u(x) = \\frac{1}{2}x(1-x)\n",
"$$\n",
"Make a partition of $n$ intervals with spacing and grid points\n",
"$$\n",
"h = \\frac{1}{n}, \\qquad x_i = i h, \\qquad i=0,1,\\ldots,n\n",
"$$\n",
"The finite difference scheme is\n",
"$$\n",
"-\\frac{u_{i-1} - 2 u_i + u_{i+1}}{h^2} = f_i, \\qquad i=1,2,\\ldots,n-1\n",
"$$"
]
},
{
"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": [
"The next function computes the residual norm defined as\n",
"$$\n",
"\\frac{ \\sum_{i=1}^{n-1} \\left[ \\frac{u_{i-1} - 2 u_i + u_{i+1}}{h^2} + f_i \\right]^2 }{ \\sum_{i=1}^{n-1}f_i^2 }\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def residual(h,u,f):\n",
" n = len(u) - 1\n",
" rnorm, fnorm = 0.0, 0.0\n",
" for i in range(1,n):\n",
" r = f[i] + (u[i-1]-2*u[i]+u[i+1])/h**2\n",
" rnorm += r**2; fnorm += f[i]**2\n",
" return np.sqrt(rnorm/fnorm)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We first define some problem parameters."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"xmin, xmax = 0.0, 1.0\n",
"n = 50\n",
"\n",
"h = (xmax - xmin)/n\n",
"x = np.linspace(xmin, xmax, n+1)\n",
"f = np.ones(n+1)\n",
"ue= 0.5*x*(1-x)\n",
"\n",
"TOL = 1.0e-6\n",
"itmax = 10000"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Gauss Jacobi method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For $m=1,2,\\ldots$\n",
"$$\n",
"u_i^m = \\frac{1}{2}[ h^2 f_i + u_{i-1}^{m-1} + u_{i+1}^{m-1}], \\qquad i=1,2,\\ldots,n-1\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of iterations = 6946\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"uold = np.zeros(n+1) # at m-1\n",
"unew = np.zeros(n+1) # at m\n",
"\n",
"for it in range(itmax):\n",
" uold[:] = unew\n",
" for i in range(1,n):\n",
" unew[i] = 0.5*(h**2*f[i] + uold[i-1] + uold[i+1])\n",
" change = residual(h,unew,f)\n",
" if change < TOL:\n",
" break\n",
"\n",
"print(\"Number of iterations = %d\" % it)\n",
"\n",
"plt.plot(x,ue,x,unew)\n",
"plt.legend((\"Exact\",\"Numerical\"));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Gauss-Seidel method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For $m=1,2,\\ldots$\n",
"$$\n",
"u_i^m = \\frac{1}{2}[ h^2 f_i + u_{i-1}^{m} + u_{i+1}^{m-1}], \\qquad i=1,2,\\ldots,n-1\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of iterations = 3474\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"u = np.zeros(n+1)\n",
"\n",
"for it in range(itmax):\n",
" for i in range(1,n):\n",
" u[i] = 0.5*(h**2*f[i] + u[i-1] + u[i+1])\n",
" change = residual(h,u,f)\n",
" if change < TOL:\n",
" break\n",
"\n",
"print(\"Number of iterations = %d\" % it)\n",
"\n",
"plt.plot(x,ue,x,u)\n",
"plt.legend((\"Exact\",\"Numerical\"));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## SOR method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For $m=1,2,\\ldots$\n",
"$$\n",
"z_i = \\frac{1}{2}[ h^2 f_i + u_{i-1}^{m} + u_{i+1}^{m-1}], \\quad u_i^m = (1-\\omega)u_i^{m-1} + \\omega z_i, \\qquad i=1,2,\\ldots,n-1\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of iterations = 146\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"omega = 2.0/(1 + np.sin(np.pi*h))\n",
"\n",
"u = np.zeros(n+1)\n",
"\n",
"for it in range(itmax):\n",
" for i in range(1,n):\n",
" z = 0.5*(h**2*f[i] + u[i-1] + u[i+1])\n",
" u[i] = (1-omega)*u[i] + omega*z\n",
" change = residual(h,u,f)\n",
" if change < TOL:\n",
" break\n",
"\n",
"print(\"Number of iterations = %d\" % it)\n",
"\n",
"plt.plot(x,ue,x,u)\n",
"plt.legend((\"Exact\",\"Numerical\"));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## SSOR method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For $m=1,2,\\ldots$\n",
"$$\n",
"z_i = \\frac{1}{2}[ h^2 f_i + u_{i-1} + u_{i+1}], \\quad u_i = (1-\\omega)u_i + \\omega z_i, \\qquad i=1,2,\\ldots,n-1\n",
"$$\n",
"$$\n",
"z_i = \\frac{1}{2}[ h^2 f_i + u_{i-1} + u_{i+1}], \\quad u_i = (1-\\omega)u_i + \\omega z_i, \\qquad i=n-1,n-2, \\ldots,1\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of iterations = 226\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"omega = 2.0/(1 + np.sin(np.pi*h))\n",
"\n",
"u = np.zeros(n+1)\n",
"\n",
"for it in range(itmax):\n",
" for i in range(1,n): # forward loop\n",
" z = 0.5*(h**2*f[i] + u[i-1] + u[i+1])\n",
" u[i] = (1-omega)*u[i] + omega*z\n",
" for i in range(n-1,0,-1): # backward loop\n",
" z = 0.5*(h**2*f[i] + u[i-1] + u[i+1])\n",
" u[i] = (1-omega)*u[i] + omega*z\n",
" change = residual(h,u,f)\n",
" if change < TOL:\n",
" break\n",
"\n",
"print(\"Number of iterations = %d\" % it)\n",
"\n",
"plt.plot(x,ue,x,u)\n",
"plt.legend((\"Exact\",\"Numerical\"));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"SSOR is not better than optimal SOR. But it is useful as a preconditioner for CG."
]
}
],
"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.10"
}
},
"nbformat": 4,
"nbformat_minor": 4
}