{
"metadata": {
"name": "",
"signature": "sha256:eedd71d7e758bb5e91f7dc45682a009f33a1000fcff7ae29f52a0cadb4647bf0"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "heading",
"level": 1,
"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",
"collapsed": false,
"input": [
"import numpy as np\n",
"from matplotlib import pyplot as plt"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def residual(h,u,f):\n",
" n = len(u) - 1\n",
" r = np.zeros(n)\n",
" for i in range(1,n):\n",
" r[i] = f[i] + (u[i-1]-2*u[i]+u[i+1])/h**2\n",
" return np.linalg.norm(r)/np.linalg.norm(f)\n",
"\n",
"xmin, xmax = 0.0, 1.0\n",
"n = 50\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-x)\n",
"\n",
"TOL = 1.0e-6\n",
"itmax = 10000"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "heading",
"level": 2,
"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",
"collapsed": false,
"input": [
"uold = np.zeros(n+1)\n",
"unew = np.zeros(n+1)\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\"));"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Number of iterations = 6936\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"svg": [
"\n",
"\n",
"\n",
"\n"
],
"text": [
""
]
}
],
"prompt_number": 3
},
{
"cell_type": "heading",
"level": 2,
"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",
"collapsed": false,
"input": [
"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\"));"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Number of iterations = 3469\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"svg": [
"\n",
"\n",
"\n",
"\n"
],
"text": [
""
]
}
],
"prompt_number": 4
},
{
"cell_type": "heading",
"level": 2,
"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",
"collapsed": false,
"input": [
"omega = 1.5\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\"));"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Number of iterations = 1150\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"svg": [
"\n",
"\n",
"\n",
"\n"
],
"text": [
""
]
}
],
"prompt_number": 5
}
],
"metadata": {}
}
]
}