{ "metadata": { "name": "", "signature": "sha256:603c1d9f62c1fe2e9596bb9d4bada2c456d10b7d8c6fd0ac3894ca7b266609a5" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "BVP using finite difference method" ] }, { "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 yexact(x):\n", " return 1.0/(np.exp(x)+np.exp(-x))\n", "\n", "def f(y,dy):\n", " return -y + 2.0*dy**2/y\n", "\n", "def phi(y):\n", " n = len(y) - 1\n", " h = 2.0/n\n", " res = np.zeros(n-1)\n", " for i in range(1,n):\n", " dy = (y[i+1] - y[i-1])/(2.0*h)\n", " res[i-1] = (y[i-1]-2.0*y[i]+y[i+1])/h**2 - f(y[i],dy)\n", " return res\n", " \n", "def dphi(y):\n", " n = len(y) - 1\n", " h = 2.0/n\n", " res = np.zeros((n-1,n-1))\n", " for i in range(1,n):\n", " dy = (y[i+1] - y[i-1])/(2.0*h)\n", " if i > 1:\n", " res[i-1,i-2] = 1.0/h**2 + 4.0*dy/y[i]\n", " res[i-1,i-1] =-2.0/h**2 + 1.0 - 2.0*(dy/y[i])**2\n", " if i < n-1:\n", " res[i-1,i ] = 1.0/h**2 - 4.0*dy/y[i]\n", " return res" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "code", "collapsed": false, "input": [ "n = 100\n", "\n", "# Initial guess for y\n", "y = np.zeros(n+1)\n", "y[:] = 1.0/(exp(1) + exp(-1))\n", "\n", "maxiter = 100\n", "TOL = 1.0e-6\n", "it = 0\n", "while it < maxiter:\n", " b = phi(y)\n", " if np.linalg.norm(b) < TOL:\n", " break\n", " A = dphi(y)\n", " v = np.linalg.solve(A,b)\n", " y[1:n] = y[1:n] - v\n", " it += 1\n", "print \"Number of iterations = %d\" % it\n", "x = np.linspace(-1.0,1.0,n+1)\n", "plt.plot(x,y,'o',x,yexact(x))\n", "plt.legend((\"Numerical\",\"Exact\"));" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Number of iterations = 13\n" ] }, { "metadata": {}, "output_type": "display_data", "svg": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text": [ "" ] } ], "prompt_number": 3 } ], "metadata": {} } ] }