{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.integrate import odeint, quad\n",
    "from scipy.optimize import newton"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Question 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = np.arange(1, 12 + 1)\n",
    "m = np.power(2, k)\n",
    "\n",
    "\n",
    "def f(x):\n",
    "    return 1 / np.sqrt(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "a, b = 1, 2\n",
    "\n",
    "\n",
    "def compute_I(f, a, b, m):\n",
    "    h_list = (b - a) / m\n",
    "    I = []\n",
    "    errors = []\n",
    "    sol_exact = quad(f, a, b)[0]\n",
    "\n",
    "    for h in h_list:\n",
    "        t = np.arange(a, b, h)\n",
    "        y = np.array(\n",
    "            [3 / 4 * h * f(t[i] + h / 3) + h / 4 * f(t[i] + h) for i in range(len(t))]\n",
    "        )\n",
    "        I_approx = np.sum(y)\n",
    "        I.append(I_approx)\n",
    "        errors.append(np.abs(I_approx - sol_exact))\n",
    "\n",
    "    return I, h_list, errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I1 = [0.8281067086445473, 0.8283842224576623, 0.8284216302273131, 0.8284264316763478, 0.8284270377897814, 0.8284271138586767, 0.828427123384198, 0.8284271245758775, 0.8284271247248971, 0.8284271247435282, 0.8284271247458573, 0.8284271247461485]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 59728 characters]",
      "text/plain": [
       "<Figure size 1000x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "I1, h_list, errors1 = compute_I(f, a, b, m)\n",
    "\n",
    "print(f\"I1 = {I1}\")\n",
    "\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.plot(np.log(h_list), np.log(errors1), \"o-\", label=\"Erreur numérique\")\n",
    "plt.plot(np.log(h_list), 2 * np.log(h_list), \"--\", label=\"Ordre 2\")\n",
    "plt.plot(np.log(h_list), 4 * np.log(h_list), \"--\", label=\"Ordre 4\")\n",
    "plt.xlabel(\"log(h)\")\n",
    "plt.ylabel(\"log(Erreur)\")\n",
    "plt.title(\"Convergence de la méthode d'intégration\")\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 54464 characters]",
      "text/plain": [
       "<Figure size 1000x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a, b = 0, 1\n",
    "I2, h_list, errors2 = compute_I(f, a, b, m)\n",
    "\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.plot(np.log(h_list), np.log(errors2), label=\"Approximation de l'intégrale\")\n",
    "plt.plot(np.log(h_list), np.log(h_list), \"--\", label=\"h\")\n",
    "plt.plot(np.log(h_list), 2 * np.log(h_list), \"--\", label=\"h^2\")\n",
    "plt.xlabel(\"h\")\n",
    "plt.ylabel(\"Approximation de l'intégrale\")\n",
    "plt.title(\"Approximation de l'intégrale par la méthode de Simpson\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dans le premier cas, l'ordre de convergence est **2** alors que dans le second cas, l'ordre de convergence est **1**. Cela est du au fait que lorsque l'on se rapproche trop de 0, donc quand l'intégrale est entre 0 et 1 et que h devient trop petit, les approximations sont trop approximatives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Question 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RKI(f, y0, vt, tol=1e-6, itmax=20):\n",
    "    N, T = len(vt), vt[-1]\n",
    "    yn = np.zeros((len(y0), N))\n",
    "    yn[:, 0] = y0\n",
    "    h = T / N\n",
    "\n",
    "    for n in range(N - 1):\n",
    "        p1 = f(vt[n], yn[:, n])\n",
    "        def F1(p2):\n",
    "            return f(vt[n] + h / 3, yn[:, n] + h / 6 * (p1 + p2)) - p2\n",
    "        p2 = newton(F1, yn[:, n], fprime=None, tol=tol, maxiter=itmax)\n",
    "        def F2(yn1):\n",
    "            return yn[:, n] + h / 4 * (3 * p2 + f(vt[n + 1], yn1)) - yn1\n",
    "        yn[:, n + 1] = newton(F2, yn[:, n], fprime=None, tol=tol, maxiter=itmax)\n",
    "    return yn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Question 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 35696 characters]",
      "text/plain": [
       "<Figure size 1000x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error with h=0.05: 0.00515308673377779\n"
     ]
    }
   ],
   "source": [
    "a, b = [0, 2]\n",
    "\n",
    "\n",
    "def f(t, y):\n",
    "    return t * np.power(y, 3) - t * y\n",
    "\n",
    "\n",
    "y0 = [0.5]\n",
    "\n",
    "\n",
    "def sol_exact(t):\n",
    "    return 1 / np.sqrt(1 + 3 * np.exp(np.power(t, 2)))\n",
    "\n",
    "\n",
    "x_fine = np.linspace(a, b, 1000)\n",
    "y_fine = sol_exact(x_fine)\n",
    "\n",
    "h = 0.05\n",
    "vt = np.arange(a, b, h)\n",
    "y = RKI(f, y0, vt)\n",
    "\n",
    "y_exact_interp = np.interp(vt, x_fine, y_fine)\n",
    "\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.plot(x_fine, y_fine, label=\"Solution exacte\")\n",
    "plt.scatter(vt, y, label=\"Solution numérique\", color=\"red\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "error = np.max(np.abs(y - y_exact_interp))\n",
    "print(f\"Error with h={h}: {error}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 52564 characters]",
      "text/plain": [
       "<Figure size 1000x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Convergence order: 1.03\n"
     ]
    }
   ],
   "source": [
    "k = np.arange(1, 10 + 1)\n",
    "h_list = 1 / np.power(2, k)\n",
    "\n",
    "errors = []\n",
    "for h in h_list:\n",
    "    vt = np.arange(a, b, h)\n",
    "    y = RKI(f, y0, vt)\n",
    "    y_exact = sol_exact(vt)\n",
    "    errors.append(np.max(np.abs(y - y_exact)))\n",
    "\n",
    "log_h = np.log(h_list)\n",
    "log_errors = np.log(errors)\n",
    "order = np.polyfit(log_h, log_errors, 1)[0]\n",
    "\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.plot(log_h, log_errors, \"o-\", label=f\"Erreur (ordre {order:.2f})\")\n",
    "plt.plot(log_h, log_h, \"--\", label=\"h\")\n",
    "plt.plot(log_h, 2 * log_h, \"--\", label=\"h^2\")\n",
    "plt.plot(log_h, 4 * log_h, \"--\", label=\"h^4\")\n",
    "plt.xlabel(\"log(h)\")\n",
    "plt.ylabel(\"log(error)\")\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.show()\n",
    "\n",
    "print(f\"Convergence order: {order:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### L'ordre de convergence de ce schéma est **1**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Question 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 169964 characters]",
      "text/plain": [
       "<Figure size 1500x500 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def F(t, Y):\n",
    "    x, y, z = Y\n",
    "    return np.array(\n",
    "        [\n",
    "            1 + np.power(x, 2) * y - (z + 1) * x,\n",
    "            x * z - np.power(x, 2) * y,\n",
    "            -x * z + 1.45,\n",
    "        ]\n",
    "    )\n",
    "\n",
    "\n",
    "h = 0.025\n",
    "y0 = np.array([1, 1, 1])\n",
    "t = np.arange(0, 40, h)\n",
    "\n",
    "sol_exact = odeint(F, y0, t, tfirst=True)\n",
    "y = RKI(F, y0, t)\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 5))\n",
    "\n",
    "ax1.scatter(y[0], y[1], label=\"Solution numérique\", color=\"red\")\n",
    "ax1.plot(sol_exact[:, 0], sol_exact[:, 1], label=\"Solution exacte\", color=\"blue\")\n",
    "ax1.legend()\n",
    "ax1.set_title(\"x vs y\")\n",
    "\n",
    "ax2.scatter(y[1], y[2], label=\"Solution numérique\", color=\"red\")\n",
    "ax2.plot(sol_exact[:, 1], sol_exact[:, 2], label=\"Solution exacte\", color=\"blue\")\n",
    "ax2.legend()\n",
    "ax2.set_title(\"y vs z\")\n",
    "\n",
    "ax3.scatter(y[0], y[2], label=\"Solution numérique\", color=\"red\")\n",
    "ax3.plot(sol_exact[:, 0], sol_exact[:, 2], label=\"Solution exacte\", color=\"blue\")\n",
    "ax3.legend()\n",
    "ax3.set_title(\"x vs z\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 35452 characters]",
      "text/plain": [
       "<Figure size 1000x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def R(z):\n",
    "    return (1 + 3 / 4 * z * (1 + z / 6) / (1 - z / 6)) / (1 - z / 4)\n",
    "\n",
    "\n",
    "x = np.linspace(-15, 5, 100)\n",
    "y = np.linspace(-7.5, 7.5, 100)\n",
    "X, Y = np.meshgrid(x, y)\n",
    "Z = R(X + 1j * Y)\n",
    "plt.figure(figsize=(10, 7))\n",
    "plt.contour(X, Y, np.abs(Z), levels=[1], cmap=\"rainbow\")\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "studies",
   "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.13.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
