7_Boundary_Value_Problems.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ★ Boundary Value Problems ★"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import modules\n",
    "import numpy as np\n",
    "import scipy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7.1 Shooting Method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "Apply the shooting method to the boundary value problem\n",
    "\n",
    "$$\n",
    "\\left\\{\\begin{matrix}\\begin{align*}\n",
    " & y'' = 4y \\\\\n",
    " & y(0) = 1 \\\\\n",
    " & y(1) = 3\n",
    "\\end{align*}\\end{matrix}\\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 383,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ode_rkf45(f, a, b, y0, h = 1e-3, tol = 1e-6):\n",
    "    w = y0.astype(np.float64)\n",
    "    t = a\n",
    "    while(t < b):\n",
    "        w_this, t_this = w, t\n",
    "        s1 = f(t, w)\n",
    "        hs1 = h * s1\n",
    "        s2 = f(t + h / 4, w + hs1 / 4)\n",
    "        hs2 = h * s2\n",
    "        s3 = f(t + 3 / 8 * h, w + 3 / 32 * hs1 + 9 / 32 * hs2)\n",
    "        hs3 = h * s3\n",
    "        s4 = f(t + 12 / 13 * h, w + 1932 / 2197 * hs1 - 7200 / 2197 * hs2 + 7296 / 2197 * hs3)\n",
    "        hs4 = h * s4\n",
    "        s5 = f(t + h, w + 439 / 216 * hs1 - 8 * hs2 + 3680 / 513 * hs3 - 845 / 4104 * hs4)\n",
    "        hs5 = s5 * h\n",
    "        s6 = f(t + h / 2, w - 8 / 27 * hs1 + 2 * hs2 - 3544 / 2565  * hs3 + 1859 / 4104 * hs4 - 11 / 40 * hs5)\n",
    "        z = w + h * (16 / 135 * s1 + 6656 / 12825 * s3 + 28561 / 56430 * s4 - 9 / 50 * s5 + 2 / 55 * s6)\n",
    "        w += h * (25 / 216 * s1 + 1408 / 2565 * s3 + 2197 / 4104 * s4 - s5 / 5)\n",
    "        t += h\n",
    "        e = abs(z - w)\n",
    "        if np.all(e / np.abs(w) < tol):\n",
    "            w = z\n",
    "        else:\n",
    "            h = 0.8 * pow(tol * abs(w) / e, 1 / 5) * h\n",
    "        \n",
    "    return w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 384,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.42030605]\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import fsolve\n",
    "\n",
    "f = lambda t, w : np.array([w[1], 4 * w[0]])\n",
    "\n",
    "def F(s):\n",
    "    return (ode_rkf45(f, 0, 1, np.array([1, s])) - 3)[0]\n",
    "\n",
    "x = fsolve(F, 0)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 385,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.42030605]\n"
     ]
    }
   ],
   "source": [
    "from scipy.integrate import ode\n",
    "from scipy.optimize import fsolve\n",
    "\n",
    "def F(s):\n",
    "    a = 0\n",
    "    ya = 1.0\n",
    "    b = 1\n",
    "    step = 100\n",
    "    dt = (b - a) / step\n",
    "    f = lambda t, w : np.array([w[1], 4 * w[0]])\n",
    "    r = ode(f).set_integrator('dopri5').set_initial_value(np.array([ya, s]), 0)\n",
    "    while r.successful() and step > 0:\n",
    "        x = r.integrate(r.t + dt)[0]\n",
    "        step -= 1\n",
    "        \n",
    "    return x - 3\n",
    "\n",
    "x = fsolve(F, 0)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7.2 Finite Difference Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "y'(t) &=  \\frac{y(t+h) - y(t-h)}{2h} - \\frac{h^2}{6}y'''(c) \\\\\n",
    "y''(t) &= \\frac{y(t+h) - 2y(t) + y(t-h)}{h^2} + \\frac{h^2}{12}f''''(c)\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example \n",
    "\n",
    "Solve the BVP\n",
    "\n",
    "$$\n",
    "\\left\\{\\begin{matrix}\\begin{align*}\n",
    "y'' &= 4y \\\\\n",
    "y(0) &= 1 \\\\\n",
    "y(1) &= 3\n",
    "\\end{align*}\\end{matrix}\\right.\n",
    "$$\n",
    "\n",
    "using finite differences\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    " & \\frac{w_{i+1} - 2w_{i} + w_{i-1}}{h^2} - 4w_i = 0 \\\\\n",
    " & \\rightarrow w_{i-1} + (-4h^2 - 2)w_i + w_{i_+1} = 0 \\\\\n",
    " & \\Rightarrow \\\\\n",
    " & 1 + (-4h^2 - 2)w_1 + w_2 = 0 \\\\\n",
    " & w_1 + (-4h^2 - 2)w_2 + w_3 = 0 \\\\\n",
    " & w_2 + (-4h^2 - 2)w_3 + 3 = 0 \\\\\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "-\\frac{9}{4} & 1 & 0 \\\\ \n",
    "1 & -\\frac{9}{4} & 1 \\\\ \n",
    "0 & 1 & -\\frac{9}{4} \n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "w_1 \\\\ \n",
    "w_2 \\\\ \n",
    "w_3  \n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "-1 \\\\ \n",
    "0 \\\\ \n",
    "-3  \n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.02494331  1.30612245  1.9138322 ]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([-9/4, 1, 0, 1, -9/4, 1, 0, 1, -9/4]).reshape(3, 3)\n",
    "b = np.array([-1, 0, -3])\n",
    "x = np.linalg.solve(A, b)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7.3 Collocation and the Finite Element Method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "Apply the Finite Element Method to the BVP\n",
    "\n",
    "$$\n",
    "\\left\\{\\begin{matrix}\\begin{align*}\n",
    "y''  &= 4y \\\\\n",
    "y(0) &= 1  \\\\\n",
    "y(1) &= 3  \n",
    "\\end{align*}\\end{matrix}\\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.01091223,  1.2855407 ,  1.89552762])"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def finite_element_method_sp( a, b, ya, yb, n):\n",
    "    h = (b - a) / (n + 1)\n",
    "    alpha = (8 / 3) * h + 2 / h\n",
    "    beta = (2 / 3) * h - 1 / h\n",
    "    A = np.zeros(n * n).reshape(n, n)\n",
    "    np.fill_diagonal(A, alpha)\n",
    "    dia_range = np.arange(n - 1)\n",
    "    A[dia_range, dia_range + 1] = beta\n",
    "    A[dia_range + 1, dia_range] = beta\n",
    "    b = np.zeros(n)\n",
    "    b[0] = -ya * beta\n",
    "    b[-1] = -yb * beta\n",
    "    x = np.linalg.solve(A, b)\n",
    "    return x\n",
    "    \n",
    "finite_element_method_sp(0, 1, 1, 3, 3)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}