{
"metadata": {
"name": "",
"signature": "sha256:015b45719dbd9ed337c062ab041572a9a273a4432ecf82275576cd665e0d3e3d"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "heading",
"level": 1,
"metadata": {},
"source": [
"Solution of 2-D Poisson equation: CG method"
]
},
{
"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",
"$$\n",
"This set of equations if of the form\n",
"$$\n",
"A U = b\n",
"$$\n",
"where\n",
"$$\n",
"A \\in \\Re^{(n-2)^2 \\times (n-2)^2}, \\qquad U,b \\in \\Re^{(n-2)^2}\n",
"$$"
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"CG Algorithm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Set initial guess $U_0 = 0$, $r_0 = b - A U_0 = b$, $p_0 = 0$\n",
"* For $k=0,1,\\ldots$\n",
" * If $\\| r_k \\| < TOL \\cdot \\|f\\|$, then stop\n",
" * If $k=0$, $\\beta_1 = 0$\n",
" * If $k > 0$, $\\beta_{k+1} = \\frac{r_k^\\top r_k}{r_{k-1}^\\top r_{k-1}}$\n",
" * $p_{k+1} = r_k + \\beta_{k+1} p_k$\n",
" * $\\alpha_{k+1} = \\frac{r_k^\\top r_k}{p_{k+1}^\\top A p_{k+1}}$\n",
" * $U_{k+1} = U_k + \\alpha_{k+1} p_{k+1}$\n",
" * $r_{k+1} = r_k - \\alpha_{k+1} p_{k+1}$"
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Code"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import numpy as np\n",
"from matplotlib import pyplot as plt"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 10
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This function computes matrix-vector product."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def ax(h,u):\n",
" n = u.shape[0]\n",
" res = np.zeros((n,n))\n",
" for i in range(1,n-1):\n",
" for j in range(1,n-1):\n",
" res[i,j] = -(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",
" return res"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 11
},
{
"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\n",
"\n",
"u = np.zeros((n,n))\n",
"p = np.zeros((n,n))\n",
"res = np.array(f)\n",
"\n",
"# First and last grid point, solution is fixed to zero.\n",
"# Hence we make residual zero, in which case solution\n",
"# will not change at these points.\n",
"res[0,:] = 0.0\n",
"res[n-1,:] = 0.0\n",
"res[:,0] = 0.0\n",
"res[:,n-1] = 0.0\n",
"\n",
"f_norm = np.linalg.norm(f)\n",
"res_old, res_new = 0.0, 0.0\n",
"for it in range(itmax):\n",
" res_new = np.linalg.norm(res)\n",
" if res_new < TOL * f_norm:\n",
" break\n",
" if it==0:\n",
" beta = 0.0\n",
" else:\n",
" beta = res_new**2 / res_old**2\n",
" p = res + beta * p\n",
" ap= ax(h,p)\n",
" alpha = res_new**2 / np.sum(p*ap)\n",
" u += alpha * p\n",
" res -= alpha * ap\n",
" res_old = res_new\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 = 31\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"svg": [
"\n",
"\n",
"\n",
"\n"
],
"text": [
""
]
}
],
"prompt_number": 13
}
],
"metadata": {}
}
]
}