{
"metadata": {
"name": "",
"signature": "sha256:0145620980b16e919fbd0b1f9c9d4a8e5eb7eea6fea6a9e6005ebd5832752bf6"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "heading",
"level": 1,
"metadata": {},
"source": [
"Conjugate gradient method"
]
},
{
"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",
"\\begin{eqnarray*}\n",
"u_0 &=& 0 \\\\\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 &=& 0\n",
"\\end{eqnarray*}\n",
"We have a matrix equation\n",
"$$\n",
"Au = f\n",
"$$"
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Algorithm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Set initial guess $u_0 = 0$, $r_0 = f - A u_0 = f$, $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": 1
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This function computes the matrix-vector product."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def ax(h,u):\n",
" n = len(u) - 1\n",
" r = np.zeros(n+1)\n",
" for i in range(1,n):\n",
" r[i] = -(u[i-1]-2*u[i]+u[i+1])/h**2\n",
" return r"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"xmin, xmax = 0.0, 1.0\n",
"n = 100\n",
"h = (xmax - xmin)/n\n",
"\n",
"x = np.linspace(0.0, 1.0, n+1)\n",
"f = np.ones(n+1)\n",
"ue= 0.5*x*(1.0-x)\n",
"\n",
"TOL = 1.0e-6\n",
"itmax = 100\n",
"\n",
"u = np.zeros(n+1)\n",
"p = np.zeros(n+1)\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] = 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",
" print it, res_new\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 / p.dot(ap)\n",
" u += alpha * p\n",
" res -= alpha * ap\n",
" res_old = res_new\n",
"\n",
"print \"Number of iterations = %d\" % it\n",
"plt.plot(x,ue,x,u)\n",
"plt.legend((\"Exact\",\"Numerical\"));"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0 9.94987437107\n",
"1 69.2928567747\n",
"2 67.8785680462\n",
"3 66.4642761188\n",
"4 65.049980784\n",
"5 63.6356818145\n",
"6 62.2213789625\n",
"7 60.8070719571\n",
"8 59.3927605016\n",
"9 57.9784442703\n",
"10 56.5641229049\n",
"11 55.1497960105\n",
"12 53.7354631505\n",
"13 52.3211238411\n",
"14 50.9067775448\n",
"15 49.4924236626\n",
"16 48.078061525\n",
"17 46.6636903813\n",
"18 45.249309387\n",
"19 43.8349175886\n",
"20 42.4205139054\n",
"21 41.0060971076\n",
"22 39.5916657897\n",
"23 38.1772183376\n",
"24 36.7627528893\n",
"25 35.3482672843\n",
"26 33.9337590019\n",
"27 32.519225083\n",
"28 31.10466203\n",
"29 29.6900656786\n",
"30 28.2754310312\n",
"31 26.8607520371\n",
"32 25.4460213\n",
"33 24.0312296814\n",
"34 22.6163657558\n",
"35 21.2014150471\n",
"36 19.7863589374\n",
"37 18.3711730709\n",
"38 16.9558249578\n",
"39 15.5402702679\n",
"40 14.1244468918\n",
"41 12.7082650271\n",
"42 11.2915897906\n",
"43 9.87420882907\n",
"44 8.45576726264\n",
"45 7.03562363974\n",
"46 5.61248608016\n",
"47 4.18330013267\n",
"48 2.73861278753\n",
"49 1.22474487139\n",
"50 1.22558586276e-14\n",
"Number of iterations = 50\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"svg": [
"\n",
"\n",
"\n",
"\n"
],
"text": [
""
]
}
],
"prompt_number": 3
}
],
"metadata": {}
}
]
}