{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Solution of 2-D Poisson equation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We consider the Poisson equation\n",
"$$\n",
" -(u_{xx} + u_{yy}) = f, \\quad\\textrm{in}\\quad \\Omega = [0,1] \\times [0,1]\n",
"$$\n",
"with boundary condition\n",
"$$\n",
"u = 0, \\quad\\textrm{on}\\quad \\partial\\Omega\n",
"$$\n",
"We will take\n",
"$$\n",
"f = 1\n",
"$$\n",
"Make a uniform grid with $n$ points in each direction, with spacing\n",
"$$\n",
"h = \\frac{1}{n-1}\n",
"$$\n",
"The grid points are\n",
"$$\n",
"x_i = ih, \\quad 0 \\le i \\le n-1\n",
"$$\n",
"$$\n",
"y_j = jh, \\quad 0 \\le j \\le n-1\n",
"$$\n",
"The finite difference approximation at $(i,j)$ is\n",
"$$\n",
"-\\frac{u_{i-1,j} - 2u_{i,j} + u_{i+1,j}}{\\Delta x^2} - \\frac{u_{i,j-1} - 2u_{i,j} + u_{i,j+1}}{\\Delta y^2} = f_{i,j}, \\qquad 1 \\le i \\le n-2, \\quad 1 \\le j \\le n-2\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 9,
"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 computers the $L^2$ norm of the residual which is used to test for convergence."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"def residual(h,u,f):\n",
" n = u.shape[0]\n",
" res = 0.0\n",
" for i in range(1,n-1):\n",
" for j in range(1,n-1):\n",
" uxx = (u[i-1,j]-2*u[i,j]+u[i+1,j])/h**2\n",
" uyy = (u[i,j-1]-2*u[i,j]+u[i,j+1])/h**2\n",
" res += (uxx + uyy + f[i,j])**2\n",
" return np.sqrt(res)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, define some problem parameters."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"xmin,xmax = 0.0,1.0\n",
"n = 21\n",
"h = (xmax - xmin)/(n-1)\n",
"x = np.linspace(xmin,xmax,n)\n",
"X,Y = np.meshgrid(x,x,indexing='ij')\n",
"plt.plot(X,Y,'.')\n",
"plt.title('Grid points')\n",
"plt.axis('equal')\n",
"f = np.ones((n,n))\n",
"TOL = 1.0e-6\n",
"itmax = 2000"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Gauss-Jacobi"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of iterations = 1102\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"uold = np.zeros((n,n))\n",
"unew = np.zeros((n,n))\n",
"res0 = residual(h,unew,f)\n",
"\n",
"for it in range(itmax):\n",
" for i in range(1,n-1):\n",
" for j in range(1,n-1):\n",
" unew[i,j] = 0.25*(h**2 * f[i,j]\n",
" + (uold[i-1,j] + uold[i+1,j] + uold[i,j-1] + uold[i,j+1]))\n",
" uold[:,:] = unew\n",
" res = residual(h,unew,f)\n",
" if res < TOL * res0:\n",
" break\n",
"\n",
"print(\"Number of iterations = %d\" % it)\n",
"plt.contourf(X,Y,unew,cmap=plt.cm.jet)\n",
"plt.axis('equal')\n",
"plt.colorbar();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Gauss-Seidel"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of iterations = 552\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"u = np.zeros((n,n))\n",
"res0 = residual(h,u,f)\n",
"\n",
"for it in range(itmax):\n",
" for i in range(1,n-1):\n",
" for j in range(1,n-1):\n",
" u[i,j] = 0.25*(h**2 * f[i,j]\n",
" + (u[i-1,j] + u[i+1,j] + u[i,j-1] + u[i,j+1]))\n",
" res = residual(h,u,f)\n",
" if res < TOL * res0:\n",
" break\n",
"\n",
"print(\"Number of iterations = %d\" % it)\n",
"plt.contourf(X,Y,u,cmap=plt.cm.jet)\n",
"plt.axis('equal')\n",
"plt.colorbar();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## SOR"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of iterations = 58\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"omega= 2.0/(1 + np.sin(np.pi*h))\n",
"u = np.zeros((n,n))\n",
"res0 = residual(h,u,f)\n",
"\n",
"for it in range(itmax):\n",
" for i in range(1,n-1):\n",
" for j in range(1,n-1):\n",
" z = 0.25*(h**2 * f[i,j]\n",
" + (u[i-1,j] + u[i+1,j] + u[i,j-1] + u[i,j+1]))\n",
" u[i,j] = (1.0-omega)*u[i,j] + omega*z\n",
" res = residual(h,u,f)\n",
" if res < TOL * res0:\n",
" break\n",
"\n",
"print(\"Number of iterations = %d\" % it)\n",
"plt.contourf(X,Y,u,cmap=plt.cm.jet)\n",
"plt.axis('equal')\n",
"plt.colorbar();"
]
}
],
"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.9.10"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 1
}