Homeworks/5_5_7_Runge-Kutta-Verner_method.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Amirhossein Mahmoudi](https://ammahmoudi.github.io)\n",
    "# Runge-Kutta-Verner method\n",
    "solution to [Numerical Analysis (R. Burden, Aires, A. Burden) [10th Ed.]](https://www.amazon.com/Richard-Burden-Numerical-Analysis-Hardcover/dp/B00SB3UL20) Section 5.5 question 7\n",
    "\n",
    "![7](images/5_5_7.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Runge-Kutta-Verner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "def runge_kutta_verner(f,a,b,alpha,tol,h_max,h_min):\n",
    "    t=a\n",
    "    w=alpha\n",
    "    h=h_max\n",
    "    flag=True\n",
    "    result=[]\n",
    "    result.append([t,w,0,0])\n",
    "    while(flag and t<b):\n",
    "      #  print(\"h:\",h)\n",
    "        k1=h*f(t,w)\n",
    "\n",
    "        k2=h*f(t+h/6.0,w+k1/6.0)\n",
    "\n",
    "        k3=h*f(t+4.0*h/15.0,w+4*k1/75.0+16*k2/75.0)\n",
    "        \n",
    "        k4=h*f(t+2*h/3.0,w+5*k1/6.0-8*k2/3.0+5*k3/2.0)\n",
    "\n",
    "        k5=h*f(t+5*h/6.0,w-165*k1/64.0+55*k2/6.0-425*k3/64.0+85*k4/96.0)\n",
    "\n",
    "        k6=h*f(t+h,w+12*k1/5.0-8*k2+4015*k3/612.0-11*k4/36.0+88*k5/255.0)\n",
    "\n",
    "        k7=h*f(t+h/15.0,w-8263*k1/15000+124*k2/75.0-643*k3/680.0-81*k4/250.0+2484*k5/10625.0)\n",
    "\n",
    "        k8=h*f(t+h,w+3501*k1/1720.0-300*k2/43.0+297275*k3/52632.0-319*k4/2322.0+24068*k5/84065.0+3850*k7/26703.0)\n",
    "\n",
    "        R=abs(-k1/160.0-125*k3/17952.0+k4/144.0-12*k5/1955.0-3*k6/44.0+125*k7/11592.0+43*k8/616.0)/h\n",
    "        if(R<=tol):\n",
    "            t=t+h\n",
    "            w=w+13*k1/160.0+2375*k3/5984.0+5*k4/16.0+12*k5/85.0+3*k6/44.0\n",
    "            result.append([t,w,h,R])\n",
    "\n",
    "        delta=0.871*np.exp(0.20*np.log(tol/R))\n",
    "        if delta<=0.1:\n",
    "            h=0.1*h\n",
    "        elif delta>=4 : h=4.0*h\n",
    "        else :h=delta*h\n",
    "        if h>h_max:h=h_max\n",
    "        if t>=b:flag=False\n",
    "        elif t+h>b:h=b-t\n",
    "        elif h<h_min:\n",
    "            flag=False\n",
    "            print(\"minimum h exceeded.\")\n",
    "            \n",
    "\n",
    "\n",
    "             \n",
    "\n",
    "    return result"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(t,y):\n",
    "    return y/t-y**2/t**2\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.0, 1.0, 0, 0]\n",
      "[1.4208756474753723, 1.0514977474344094, 0.4208756474753723, 9.82603971433408e-07]\n",
      "[1.7887472375470532, 1.1310235710192327, 0.36787159007168085, 8.242650384651363e-08]\n",
      "[2.288747237547053, 1.2520370925446567, 0.5, 2.9124264472324413e-08]\n",
      "[2.788747237547053, 1.3767467259050725, 0.5, 5.224624621824603e-10]\n",
      "[3.288747237547053, 1.5013540148047648, 0.5, 5.352158508054039e-10]\n",
      "[3.788747237547053, 1.624642694007602, 0.5, 2.8299814228138764e-10]\n",
      "[4.0, 1.6762292227022868, 0.211252762452947, 2.1288699755221945e-12]\n"
     ]
    }
   ],
   "source": [
    "a=1.0\n",
    "b=4.0\n",
    "y_a=1.0\n",
    "tol=0.000001\n",
    "h_min=0.05\n",
    "h_max=0.5\n",
    "result=runge_kutta_verner(f,a,b,y_a,tol,h_max,h_min)\n",
    "for r in result:\n",
    "    print(r)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exact answer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.         1.05150869 1.13103364 1.25204675 1.37675632 1.50136369\n",
      " 1.62465253 1.67623914]\n"
     ]
    }
   ],
   "source": [
    "def y_real(t):\n",
    "    return t/(1+np.log(t))\n",
    "\n",
    "exact_answers=y_real(np.array(result)[:,0])\n",
    "print(exact_answers)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Final result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>t</th>\n",
       "      <th>h</th>\n",
       "      <th>R</th>\n",
       "      <th>rkv</th>\n",
       "      <th>exact</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.420876</td>\n",
       "      <td>0.420876</td>\n",
       "      <td>9.826040e-07</td>\n",
       "      <td>1.051498</td>\n",
       "      <td>1.051509</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.788747</td>\n",
       "      <td>0.367872</td>\n",
       "      <td>8.242650e-08</td>\n",
       "      <td>1.131024</td>\n",
       "      <td>1.131034</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2.288747</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>2.912426e-08</td>\n",
       "      <td>1.252037</td>\n",
       "      <td>1.252047</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2.788747</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>5.224625e-10</td>\n",
       "      <td>1.376747</td>\n",
       "      <td>1.376756</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>3.288747</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>5.352159e-10</td>\n",
       "      <td>1.501354</td>\n",
       "      <td>1.501364</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>3.788747</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>2.829981e-10</td>\n",
       "      <td>1.624643</td>\n",
       "      <td>1.624653</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>4.000000</td>\n",
       "      <td>0.211253</td>\n",
       "      <td>2.128870e-12</td>\n",
       "      <td>1.676229</td>\n",
       "      <td>1.676239</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          t         h             R       rkv     exact\n",
       "0  1.000000  0.000000  0.000000e+00  1.000000  1.000000\n",
       "1  1.420876  0.420876  9.826040e-07  1.051498  1.051509\n",
       "2  1.788747  0.367872  8.242650e-08  1.131024  1.131034\n",
       "3  2.288747  0.500000  2.912426e-08  1.252037  1.252047\n",
       "4  2.788747  0.500000  5.224625e-10  1.376747  1.376756\n",
       "5  3.288747  0.500000  5.352159e-10  1.501354  1.501364\n",
       "6  3.788747  0.500000  2.829981e-10  1.624643  1.624653\n",
       "7  4.000000  0.211253  2.128870e-12  1.676229  1.676239"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "table=pd.DataFrame(data=np.c_[result,exact_answers],columns=['t','rkv','h','R','exact'])\n",
    "table=table[['t','h','R','rkv','exact']]\n",
    "table.head(n=10)\n"
   ]
  }
 ],
 "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.10.8"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}