{
"metadata": {
"name": "",
"signature": "sha256:a0162e732ce49cf3f32f8764f264cf7b63a04d287536e3e1b8417b87e16b0d78"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "heading",
"level": 1,
"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",
"collapsed": false,
"input": [
"import numpy as np\n",
"from matplotlib import pyplot as plt"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 6
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def residual(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",
" r = -(u[i-1,j]-2*u[i,j]+u[i+1,j])/h**2 \\\n",
" - (u[i,j-1]-2*u[i,j]+u[i,j+1])/h**2 \\\n",
" - f[i,j]\n",
" res += r**2\n",
" return np.sqrt(res)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 7
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"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)\n",
"f = np.ones((n,n))\n",
"\n",
"TOL = 1.0e-6\n",
"itmax = 2000"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 8
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Gauss-Jacobi"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"uold = np.zeros((n,n))\n",
"unew = np.zeros((n,n))\n",
"res0 = residual(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(unew,f)\n",
" if res < TOL * res0:\n",
" break\n",
"\n",
"print \"Number of iterations = %d\" % it\n",
"cs = plt.contourf(X,Y,unew)\n",
"plt.colorbar(cs);"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Number of iterations = 1102\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"svg": [
"\n",
"\n",
"\n",
"\n"
],
"text": [
""
]
}
],
"prompt_number": 9
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Gauss-Seidel"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"u = np.zeros((n,n))\n",
"res0 = residual(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(u,f)\n",
" if res < TOL * res0:\n",
" break\n",
"\n",
"print \"Number of iterations = %d\" % it\n",
"cs = plt.contourf(X,Y,u)\n",
"plt.colorbar(cs);"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Number of iterations = 552\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"svg": [
"\n",
"\n",
"\n",
"\n"
],
"text": [
""
]
}
],
"prompt_number": 10
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"SOR"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"omega= 1.5\n",
"u = np.zeros((n,n))\n",
"res0 = residual(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(u,f)\n",
" if res < TOL * res0:\n",
" break\n",
"\n",
"print \"Number of iterations = %d\" % it\n",
"cs = plt.contourf(X,Y,u)\n",
"plt.colorbar(cs);"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Number of iterations = 177\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"svg": [
"\n",
"\n",
"\n",
"\n"
],
"text": [
""
]
}
],
"prompt_number": 11
}
],
"metadata": {}
}
]
}