{ "metadata": { "name": "", "signature": "sha256:86cd5003ecced9014042c33f9057ae2c2b256a57d86abe7970a94cf5d7767723" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "Shooting method for BVP" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider the BVP\n", "\n", "$$\n", "y'' = -y + \\frac{2 (y')^2}{y}, \\qquad -1 < x < +1\n", "$$\n", "with boundary conditions\n", "$$\n", "y(-1) = y(+1) = (e + e^{-1})^{-1}\n", "$$\n", "The exact solution is\n", "$$\n", "y(x) = (e^x + e^{-x})^{-1}\n", "$$" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Formulate as an initial value problem" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "y'' = -y + \\frac{2 (y')^2}{y}, \\qquad -1 < x < +1\n", "$$\n", "with initial conditions\n", "$$\n", "y(-1) = (e + e^{-1})^{-1}, \\qquad y'(-1) = s\n", "$$\n", "Let the solution of this IVP be denoted as $y(x;s)$. We have to determine $s$ so that\n", "$$\n", "\\phi(s) = y(1;s) - (e + e^{-1})^{-1} = 0\n", "$$" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Newton method" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To find the root of $\\phi(s)$, we use Newton method with an initial guess $s_0$. Then the Newton method updates the guess by\n", "$$\n", "s_{m+1} = s_m - \\frac{\\phi(s_m)}{\\frac{d}{ds}\\phi(s_m)}, \\qquad m=0,1,2,\\ldots\n", "$$\n", "Define\n", "$$\n", "z_s(x) = \\frac{\\partial}{\\partial s} y(x;s)\n", "$$\n", "Then\n", "$$\n", "z_s'(x) = \\frac{\\partial}{\\partial s} y'(x;s)\n", "$$\n", "The derivative of $\\phi(s)$ is given by\n", "$$\n", "\\frac{d}{ds}\\phi(s) = \\frac{\\partial}{\\partial s} y(1;s) = z_s(1)\n", "$$\n", "We can write an equation for $z_s(x)$\n", "$$\n", "z_s'' = \\left[ -1 - 2 \\left( \\frac{y'}{y} \\right)^2 \\right] z_s + 4 \\frac{y'}{y} z_s'\n", "$$\n", "with initial conditions\n", "$$\n", "z_s(-1) = 0, \\qquad z_s'(-1) = 1\n", "$$" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "First order ODE system" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the vector\n", "$$\n", "u = \\begin{bmatrix}\n", "u_1 \\\\ u_2 \\\\ u_3 \\\\ u_4 \\end{bmatrix} = \\begin{bmatrix}\n", "y \\\\ y' \\\\ z_s \\\\ z_s' \\end{bmatrix}\n", "$$\n", "Then\n", "$$\n", "u_1' = u_2, \\qquad u_2' = -u_1 + \\frac{2 u_2^2}{u_1}, \\qquad u_3' = u_4, \\qquad u_4' = \\left[ -1 - 2 \\left( \\frac{u_2}{u_1} \\right)^2 \\right] u_3 + 4 \\frac{u_2}{u_1} u_4\n", "$$\n", "Hence we get the first order ODE system\n", "$$\n", "u' = f(u) = \\begin{bmatrix}\n", "u_2 \\\\\n", "-u_1 + \\frac{2 u_2^2}{u_1} \\\\\n", "u_4 \\\\\n", "\\left[ -1 - 2 \\left( \\frac{u_2}{u_1} \\right)^2 \\right] u_3 + 4 \\frac{u_2}{u_1} u_4\n", "\\end{bmatrix}\n", "$$\n", "with initial condition\n", "$$\n", "u(-1) = \\begin{bmatrix}\n", "(e+e^{-1})^{-1} \\\\\n", "s \\\\\n", "0 \\\\\n", "1\n", "\\end{bmatrix}\n", "$$\n", "Once we solve this IVP, we get\n", "$$\n", "\\phi(s) = u_1(1) - (e + e^{-1})^{-1}, \\qquad \\frac{d}{ds}\\phi(s) = z_s(1) = u_3(1)\n", "$$" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Now we start coding" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np\n", "from matplotlib import pyplot as plt" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 43 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now code the problem specific data." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def f(u):\n", " rhs = np.zeros(4)\n", " rhs[0] = u[1]\n", " rhs[1] = -u[0] + 2.0*u[1]**2/u[0]\n", " rhs[2] = u[3]\n", " rhs[3] = (-1.0-2.0*(u[1]/u[0])**2)*u[2] + 4.0*u[1]*u[3]/u[0]\n", " return rhs\n", "\n", "def initialcondition(s):\n", " u = np.zeros(4)\n", " u[0] = 1.0/(np.exp(1)+np.exp(-1))\n", " u[1] = s\n", " u[2] = 0.0\n", " u[3] = 1.0\n", " return u\n", "\n", "def yexact(x):\n", " return 1.0/(np.exp(x) + np.exp(-x))" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 44 }, { "cell_type": "code", "collapsed": false, "input": [ "def solvephi(n,s):\n", " h = 2.0/n\n", " u = np.zeros((n+1,4))\n", " u[0] = initialcondition(s)\n", " for i in range(1,n+1):\n", " u1 = u[i-1] + 0.5*h*f(u[i-1])\n", " u[i] = u[i-1] + h*f(u1)\n", " phi = u[n][0] - 1.0/(exp(1)+exp(-1))\n", " dphi= u[n][2]\n", " return phi,dphi,u" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 45 }, { "cell_type": "code", "collapsed": false, "input": [ "n = 100\n", "s = 0.2\n", "p, dp, u = solvephi(n,s)\n", "x = np.linspace(-1.0,1.0,n+1)\n", "xe = np.linspace(-1.0,1.0,1000)\n", "ye = yexact(xe)\n", "plt.plot(x,u[:,0],xe,ye)\n", "plt.legend((\"Numerical\",\"Exact\"))" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 46, "text": [ "" ] }, { "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" ], "text": [ "" ] } ], "prompt_number": 46 }, { "cell_type": "code", "collapsed": false, "input": [ "n = 100\n", "s = 0.2\n", "maxiter = 100\n", "TOL = 1.0e-8\n", "it = 0\n", "while it < maxiter:\n", " p, dp, u = solvephi(n,s)\n", " if np.abs(p) < TOL:\n", " break\n", " s = s - p/dp\n", " it += 1\n", "print \"Root = %e\" % s\n", "x = np.linspace(-1.0,1.0,n+1)\n", "xe = np.linspace(-1.0,1.0,1000)\n", "ye = yexact(xe)\n", "plt.plot(x,u[:,0],xe,ye)\n", "plt.legend((\"Numerical\",\"Exact\"))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Root = 2.467633e-01\n" ] }, { "metadata": {}, "output_type": "pyout", "prompt_number": 49, "text": [ "" ] }, { "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" ], "text": [ "" ] } ], "prompt_number": 49 }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }