{
"cells": [
{
"cell_type": "markdown",
"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": "markdown",
"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": "markdown",
"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": "markdown",
"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": "markdown",
"metadata": {},
"source": [
"## Now we start coding"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"import numpy as np\n",
"from matplotlib import pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now code the problem specific data."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"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))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The function below implements the two step RK scheme."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"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/(np.exp(1)+np.exp(-1))\n",
" dphi= u[n][2]\n",
" return phi,dphi,u"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us see the solution corresponding to the initial guess."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"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.grid(True)\n",
"plt.legend((\"Numerical\",\"Exact\"));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The initial guess is far from the solution. We will not start the shooting method."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1 2.712498e-01 1.113544e-01\n",
" 2 2.534755e-01 1.223885e-01\n",
" 3 2.472677e-01 2.633095e-02\n",
" 4 2.467661e-01 1.840413e-03\n",
" 5 2.467633e-01 1.033511e-05\n",
"Root s = 2.467633e-01\n",
"phi(s) = 3.296041e-10\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"n = 100\n",
"s = 0.2\n",
"maxiter = 100\n",
"TOL = 1.0e-8\n",
"it = 0\n",
"x = np.linspace(-1.0,1.0,n+1)\n",
"plt.figure(figsize=(8,5))\n",
"plt.plot(xe,ye)\n",
"leg = ['Exact']\n",
"while it < maxiter:\n",
" p, dp, u = solvephi(n,s)\n",
" plt.plot(x,u[:,0])\n",
" leg.append(str(it))\n",
" if np.abs(p) < TOL:\n",
" break\n",
" s = s - p/dp\n",
" it += 1\n",
" print('%5d %16.6e %16.6e'%(it,s,np.abs(p)))\n",
"plt.legend(leg); plt.grid(True)\n",
"print(\"Root s = %e\" % s)\n",
"print(\"phi(s) = %e\" % p)\n",
"plt.figure(figsize=(8,5))\n",
"plt.plot(x,u[:,0],xe,ye)\n",
"plt.grid(True)\n",
"plt.title('Final solution')\n",
"plt.legend((\"Numerical\",\"Exact\"));"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.14"
}
},
"nbformat": 4,
"nbformat_minor": 1
}