{
 "cells": [
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T12:51:01.523548Z",
     "start_time": "2025-03-04T12:51:01.520087Z"
    }
   },
   "cell_type": "code",
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import odeint"
   ],
   "id": "fe01d49f9c6740ce",
   "outputs": [],
   "execution_count": 259
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Exercise 2.18 : Numerical precision",
   "id": "4a91fd48ddb1c4fc"
  },
  {
   "metadata": {
    "collapsed": true,
    "ExecuteTime": {
     "end_time": "2025-03-04T12:51:01.559159Z",
     "start_time": "2025-03-04T12:51:01.555615Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def f(t, x):\n",
    "    return 2 * np.power(np.abs(x), 1 / 2)\n",
    "\n",
    "\n",
    "def f_exact(x, lam):\n",
    "    return np.power(x - lam, 2) * (x >= lam)"
   ],
   "id": "initial_id",
   "outputs": [],
   "execution_count": 260
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T12:51:01.565527Z",
     "start_time": "2025-03-04T12:51:01.562740Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def euler_explicit(N, T, y0, f, type=int, args=()):\n",
    "    h = T / N\n",
    "    t = np.linspace(0, T, N + 1)\n",
    "    y = np.zeros(N + 1, dtype=type)\n",
    "    y[0] = y0\n",
    "    for n in range(N):\n",
    "        y[n + 1] = y[n] + h * f(t[n], y[n], *args)\n",
    "    return t, y"
   ],
   "id": "f447a3bdff49e850",
   "outputs": [],
   "execution_count": 261
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T12:51:01.639877Z",
     "start_time": "2025-03-04T12:51:01.582392Z"
    }
   },
   "cell_type": "code",
   "source": [
    "T = 10\n",
    "N = 1000\n",
    "u0 = 0\n",
    "\n",
    "t, sol_appr = euler_explicit(N, T, u0, f)\n",
    "x = np.linspace(u0, T, 100)\n",
    "\n",
    "plt.scatter(t, sol_appr, label=\"Approximation\")\n",
    "plt.plot(x, f_exact(x, T), label=\"Exact solution\", color=\"red\")\n",
    "plt.legend()\n",
    "plt.show()"
   ],
   "id": "957b7a23b6076f77",
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ],
      "image/png": "[encoded data, 16688 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "execution_count": 262
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T12:51:01.657056Z",
     "start_time": "2025-03-04T12:51:01.654570Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def euler_implicit(N, T, y0):\n",
    "    h = T / N\n",
    "    t = np.linspace(0, T, N + 1)\n",
    "    y = np.zeros(N + 1)\n",
    "    y[0] = y0\n",
    "    for n in range(N):\n",
    "        y[n + 1] = np.power(h + np.sqrt(h ** 2 + y[n]), 2)\n",
    "    return t, y"
   ],
   "id": "57527c2e5eb6d17d",
   "outputs": [],
   "execution_count": 263
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T12:51:01.739113Z",
     "start_time": "2025-03-04T12:51:01.672072Z"
    }
   },
   "cell_type": "code",
   "source": [
    "T = 10\n",
    "N = 1000\n",
    "u0 = 0\n",
    "\n",
    "t, sol_appr = euler_implicit(N, T, u0)\n",
    "x = np.linspace(0, T, 100)\n",
    "\n",
    "plt.scatter(t, sol_appr, label=\"Approximation\")\n",
    "plt.plot(x, f_exact(x, T), label=\"Exact solution\", color=\"red\")\n",
    "plt.plot(x, x ** 2, label=\"Square function\", color=\"green\")\n",
    "plt.legend()\n",
    "plt.show()"
   ],
   "id": "2ff29a8872a51c89",
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ],
      "image/png": "[encoded data, 34248 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "execution_count": 264
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Exercise 2.19 : SIR model, schema order",
   "id": "ecea2cd935dea93a"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T12:51:01.757489Z",
     "start_time": "2025-03-04T12:51:01.754373Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def F(y, t, a, r):\n",
    "    S, I, R = y\n",
    "    dS = -r * S * I\n",
    "    dI = r * S * I - a * I\n",
    "    dR = a * I\n",
    "    return np.array([dS, dI, dR])"
   ],
   "id": "61588e94b56ea779",
   "outputs": [],
   "execution_count": 265
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T12:51:01.776543Z",
     "start_time": "2025-03-04T12:51:01.773086Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def euler_explicit_sir(N, T, u0, f, args):\n",
    "    t = np.linspace(0, T, N + 1)\n",
    "    y = np.zeros((N + 1, len(y0)))\n",
    "    y[0] = y0\n",
    "    h = T / N\n",
    "    for i in range(N):\n",
    "        y[i + 1] = y[i] + h * F(y[i], t[i], *args)\n",
    "    return t, y"
   ],
   "id": "a32f0facb958a1a0",
   "outputs": [],
   "execution_count": 266
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T12:51:01.794953Z",
     "start_time": "2025-03-04T12:51:01.791398Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def heun_sir(N, T, u0, f, args):\n",
    "    t = np.linspace(0, T, N + 1)\n",
    "    y = np.zeros((N + 1, len(y0)))\n",
    "    y[0] = y0\n",
    "    h = T / N\n",
    "    for i in range(N):\n",
    "        k1 = h * F(y[i], t[i], *args)\n",
    "        k2 = h * F(y[i] + k1, t[i + 1], *args)\n",
    "        y[i + 1] = y[i] + 0.5 * (k1 + k2)\n",
    "    return t, y"
   ],
   "id": "ca061c36b09f7200",
   "outputs": [],
   "execution_count": 267
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T12:51:01.814042Z",
     "start_time": "2025-03-04T12:51:01.810178Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def runge_kutta_order_4_sir(N, T, u0, f, args):\n",
    "    t = np.linspace(0, T, N + 1)\n",
    "    y = np.zeros((N + 1, len(y0)))\n",
    "    y[0] = y0\n",
    "    h = T / N\n",
    "    for i in range(N):\n",
    "        k1 = h * F(y[i], t[i], *args)\n",
    "        k2 = h * F(y[i] + 0.5 * k1, t[i] + 0.5 * h, *args)\n",
    "        k3 = h * F(y[i] + 0.5 * k2, t[i] + 0.5 * h, *args)\n",
    "        k4 = h * F(y[i] + k3, t[i + 1], *args)\n",
    "        y[i + 1] = y[i] + (k1 + 2 * k2 + 2 * k3 + k4) / 6\n",
    "    return t, y"
   ],
   "id": "afc8571b14a6b464",
   "outputs": [],
   "execution_count": 268
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T12:51:01.947339Z",
     "start_time": "2025-03-04T12:51:01.876657Z"
    }
   },
   "cell_type": "code",
   "source": [
    "t0, T = 52, 60\n",
    "N = 150\n",
    "\n",
    "y0 = np.array([10e6, 10, 0], dtype=np.float64)  # Initial conditions: S, I, R\n",
    "y0 = y0 / np.sum(y0)  # Normalize initial conditions\n",
    "a = 0.33  # Recovery rate\n",
    "r = 0.5  # Infection rate\n",
    "t = np.linspace(t0, T, N)\n",
    "\n",
    "sol = odeint(F, y0, t, args=(a, r))\n",
    "\n",
    "plt.plot(t, sol[:, 0], label='S(t)')\n",
    "plt.plot(t, sol[:, 1], label='I(t)')\n",
    "plt.plot(t, sol[:, 2], label='R(t)')\n",
    "plt.legend()\n",
    "plt.show()"
   ],
   "id": "3ff4bbd4286d2b69",
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ],
      "image/png": "[encoded data, 14796 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "execution_count": 269
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T12:51:02.453901Z",
     "start_time": "2025-03-04T12:51:01.967378Z"
    }
   },
   "cell_type": "code",
   "source": [
    "errors_EE = []\n",
    "errors_H = []\n",
    "errors_RK4 = []\n",
    "\n",
    "h_list = [0.05, 0.01, 0.1, 0.5, 1, 2, 4]\n",
    "\n",
    "\n",
    "def calculate_errors(sol_exact, sol_appr):\n",
    "    return np.max(\n",
    "        np.power(np.abs(sol_appr - sol_exact), 2)[np.isfinite(np.power(np.abs(sol_appr - sol_exact), 2))]\n",
    "    )\n",
    "\n",
    "\n",
    "for h in h_list:\n",
    "    N = int(T / h)\n",
    "    t = np.linspace(0, T, N + 1)\n",
    "\n",
    "    _, sol_appr_EE = euler_explicit_sir(N, T, y0, F, args=(a, r))\n",
    "    _, sol_appr_H = heun_sir(N, T, y0, F, args=(a, r))\n",
    "    _, sol_appr_RK4 = runge_kutta_order_4_sir(N, T, y0, F, args=(a, r))\n",
    "    sol_exact = odeint(F, y0, t, args=(a, r))\n",
    "\n",
    "    errors_EE.append(calculate_errors(sol_exact, sol_appr_EE))\n",
    "    errors_H.append(calculate_errors(sol_exact, sol_appr_H))\n",
    "    errors_RK4.append(calculate_errors(sol_exact, sol_appr_RK4))\n",
    "\n",
    "plt.plot(errors_EE, label=\"Euler Explicit\")\n",
    "plt.plot(errors_H, label=\"Heun\")\n",
    "plt.plot(errors_RK4, label=\"Runge Kutta order 4\")\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.legend()\n",
    "plt.show()"
   ],
   "id": "c147c62e0a55bf28",
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ],
      "image/png": "[encoded data, 36560 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "execution_count": 270
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## Exercise 2.20 : Stability of Euler Explicit",
   "id": "ee44f24a438fb7b6"
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T12:51:02.474489Z",
     "start_time": "2025-03-04T12:51:02.471385Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def euler_implicit(N, T, u0, lam, type=int):\n",
    "    h = T / N\n",
    "    t = np.linspace(0, T, N + 1)\n",
    "    u = np.zeros(N + 1, dtype=type)\n",
    "    u[0] = u0\n",
    "    for n in range(N):\n",
    "        u[n + 1] = u[n] / (1 - lam * h)\n",
    "    return t, u"
   ],
   "id": "4484c71f50b8e691",
   "outputs": [],
   "execution_count": 271
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T12:51:02.644858Z",
     "start_time": "2025-03-04T12:51:02.493576Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def f(t, x, lam):\n",
    "    return x * lam\n",
    "\n",
    "\n",
    "T = 10 * 2 * np.pi\n",
    "lam = 1j\n",
    "N = 400\n",
    "x0 = 1\n",
    "\n",
    "t, x_appr_EE = euler_explicit(N, T, x0, f, type=complex, args=(lam,))\n",
    "_, x_appr_EI = euler_implicit(N, T, x0, lam, type=complex)\n",
    "x_exact = x0 * np.exp(lam * t)\n",
    "\n",
    "# Plot the real parts\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.plot(t, np.real(x_appr_EI), label='Numerical Solution by EI')\n",
    "plt.plot(t, np.real(x_appr_EE), label='Numerical Solution by EE')\n",
    "plt.plot(t, np.real(x_exact), label='Exact Solution', linestyle='--')\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Real Part')\n",
    "plt.legend()\n",
    "plt.title('Real Part of the Solution')\n",
    "\n",
    "# Plot the imaginary parts\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.plot(t, np.imag(x_appr_EI), label='Numerical Solution by EI')\n",
    "plt.plot(t, np.imag(x_appr_EE), label='Numerical Solution by EE')\n",
    "plt.plot(t, np.imag(x_exact), label='Exact Solution', linestyle='--')\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Imaginary Part')\n",
    "plt.legend()\n",
    "plt.title('Imaginary Part of the Solution')\n",
    "\n",
    "plt.show()"
   ],
   "id": "1b8a3146e7060c52",
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 1200x600 with 2 Axes>"
      ],
      "image/png": "[encoded data, 88148 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "execution_count": 272
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-03-04T12:51:02.662942Z",
     "start_time": "2025-03-04T12:51:02.661576Z"
    }
   },
   "cell_type": "code",
   "source": "",
   "id": "aeb67c5158f41a0b",
   "outputs": [],
   "execution_count": null
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
