euler_rk_shooting.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0e66d527",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2d03eb8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Euler's method for equations\n",
    "def ode_Euler(f, y0, t):\n",
    "    \n",
    "    y = np.zeros(len(t))\n",
    "    y[0] = y0\n",
    "    for i in range (0, len(t) - 1):\n",
    "        y[i+1] = y[i] + f(y[i], t[i])*(t[i+1] - t[i])\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f895b906",
   "metadata": {},
   "outputs": [],
   "source": [
    "#RK4 for equations\n",
    "def rk4(x0, x, y0, h):\n",
    "    \n",
    "    N = int(((x-x0)/h))\n",
    "    y = y0\n",
    "        \n",
    "    for k in range(1, N+1):\n",
    "        k1 = h * f(x0,y)\n",
    "        k2 = h * f(x0 + 0.5 * h, y + 0.5 * k1)\n",
    "        k3 = h * f(x0 + 0.5 * h, y + 0.5 * k2)\n",
    "        k4 = h * f(x0 + h, y + k3)\n",
    "        y = y + (1/6) * (k1 + 2*k2 + 2*k3 + k4)\n",
    "        x0 += h\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "02a88d2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Euler's method for systems\n",
    "def euler(f, a, b, y, dy, N=20):\n",
    "    \n",
    "    N -= 1\n",
    "    h = (b-a)/N\n",
    "    \n",
    "    dy1 = []\n",
    "    y1 = []\n",
    "    x1 = []\n",
    "    \n",
    "    for x in np.linspace(a+h,b,N):\n",
    "        \n",
    "        dy += h * f(dy,y,x) # integration to get y'\n",
    "        y  += h * dy        # integration to get y\n",
    "        \n",
    "        dy1.append(dy)\n",
    "        y1.append(y)\n",
    "        x1.append(x)\n",
    "    \n",
    "    return y1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8f7ad4d0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.0502970622766254, 1.0985450108165298, 1.1450135214583803, 1.1899457615585987, 1.2335484641729064, 1.2759831154877042, 1.317358484913293, 1.357724693245853, 1.3970689754119887, 1.435313251497881, 1.4723135741015827, 1.507861472669681, 1.541687167569064, 1.5734645793948072, 1.6028180136038062, 1.6293303581225673, 1.6525525931645513, 1.6720143790654451, 1.6872354603375146]\n"
     ]
    }
   ],
   "source": [
    "f = lambda dy,y,x: np.sin(3*x) - y\n",
    "data = euler(f,0,1,1,1)\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6e3bb453",
   "metadata": {},
   "outputs": [],
   "source": [
    "#RK4 method for systems\n",
    "def RK4(f, a, b, y, dy, N=20):\n",
    "    \n",
    "    N -= 1\n",
    "    h = (b - a)/N\n",
    "    \n",
    "    dy1 = []\n",
    "    y1 = []\n",
    "    x1 = []\n",
    "    \n",
    "    z = lambda dy,y,x: dy #z(dy,y,x) = dy\n",
    "\n",
    "    for x in np.linspace(a+h,b,N):\n",
    "        \n",
    "        m1 = h*f(dy,y,x)\n",
    "        k1 = h*z(dy,y,x)\n",
    "        m2 = h*f(dy+m1/2,y+k1/2,x+h/2)\n",
    "        k2 = h*z(dy+m1/2,y+k1/2,x+h/2)\n",
    "        m3 = h*f(dy+m2/2,y+k2/2,x+h/2)\n",
    "        k3 = h*z(dy+m2/2,y+k2/2,x+h/2)\n",
    "        m4 = h*f(dy+m3,y+m3,x+h)\n",
    "        k4 = h*z(dy+m3,y+m3,x+h)\n",
    "        \n",
    "        dy += (m1+2*m2+2*m3+m4)/6\n",
    "        y += (k1+2*k2+2*k3+k4)/6\n",
    "        \n",
    "        dy1.append(dy)\n",
    "        y1.append(y)\n",
    "        x1.append(x)\n",
    "        \n",
    "\n",
    "    return y1[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "75002da9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.7114047635921288\n"
     ]
    }
   ],
   "source": [
    "f = lambda dy,y,x: np.sin(3*x) - y\n",
    "data = RK4(f,0,1,1,1)\n",
    "print(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f3d4c386",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Shooting method\n",
    "def shoot(f, x, y, yb, b):\n",
    "    eta = [1, 0.8]\n",
    "    phi = [RK4(f, x, b, y, eta[0])- yb, RK4(f, x, b, y, eta[1])- yb]\n",
    "    eps = min(phi)\n",
    "    #print(phi)\n",
    "    \n",
    "    while abs(eps) >= 0.0001:\n",
    "        tmp = eta[1] - phi[1]*(eta[1] - eta[0])/(phi[1] - phi[0])\n",
    "        eta[0], eta[1] = eta[1], tmp\n",
    "        #print(phi)\n",
    "        phi = [RK4(f, x, b, y, eta[0])- yb, RK4(f, x, b, y, eta[1])- yb]\n",
    "        eps = min(phi)\n",
    "        #print(eps)\n",
    "    return RK4(f, x, b, y, eta[1]), eta[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "2d30a737",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.000003248004767, -0.20123692993690498)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = lambda dy, y, x: np.exp(x) + np.sin(y)\n",
    "a = shoot(f, 0, 1, 2, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10b052cd",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}