diff --git a/tp/runge-kutta.ipynb b/tp/runge-kutta.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..2d2a98cb8a3e7c644670a346d69b6b3ff8b540c9
--- /dev/null
+++ b/tp/runge-kutta.ipynb
@@ -0,0 +1,607 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "<img src=\"https://gitlab.irit.fr/toc/ens-n7/texCoursN7/-/raw/main/LOGO_INP_N7.png\" alt=\"N7\" height=\"80\"/>\n",
+ "\n",
+ "<img src=\"https://gitlab.irit.fr/toc/ens-n7/texCoursN7/-/raw/main/logo-insa.png\" alt=\"INSA\" height=\"80\"/>\n",
+ "\n",
+ "# Méthodes de Runge-Kutta\n",
+ "\n",
+ "- Date : 2023-2024\n",
+ "- Durée approximative : 1h15\n",
+ "\n",
+ "## Introduction\n",
+ "\n",
+ "Nous allons dans ce TP, implémenter quelques méthodes de Runge-Kutta et étudier leur convergence. On considère un pas de temps $h$ uniforme."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# activate local project\n",
+ "using Pkg\n",
+ "Pkg.activate(\".\")\n",
+ "\n",
+ "# load packages\n",
+ "using Plots\n",
+ "using Plots.PlotMeasures\n",
+ "using Polynomials\n",
+ "\n",
+ "#\n",
+ "px = PlotMeasures.px;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Fonctons auxiliaires \n",
+ "\n",
+ "function method_infos(method)\n",
+ "\n",
+ " if method == :euler \n",
+ " method_func = euler\n",
+ " method_name = \"Euler\"\n",
+ " methode_stages = 1\n",
+ " elseif method == :runge \n",
+ " method_func = runge\n",
+ " method_name = \"Runge\"\n",
+ " methode_stages = 2\n",
+ " elseif method == :heun\n",
+ " method_func = heun\n",
+ " method_name = \"Heun\"\n",
+ " methode_stages = 3\n",
+ " elseif method == :rk4\n",
+ " method_func = rk4\n",
+ " method_name = \"RK4\"\n",
+ " methode_stages = 4\n",
+ " else \n",
+ " error(\"Méthode d'intégration non reconnue\")\n",
+ " end\n",
+ "\n",
+ " return method_func, method_name, methode_stages\n",
+ "\n",
+ "end\n",
+ "\n",
+ "function convergence(method, f, x0, tspan, sol, Nspan)\n",
+ "\n",
+ " # Récupération des informations sur la méthode\n",
+ " method_func, method_name, methode_stages = method_infos(method)\n",
+ " \n",
+ " # Ecriture des choix de paramètres\n",
+ " println(\"Méthode d'intégration : \", method_name)\n",
+ "\n",
+ " plts = [] # Liste des graphiques\n",
+ "\n",
+ " for N ∈ Nspan\n",
+ "\n",
+ " # Solution approchée\n",
+ " ts, xs = method_func(f, x0, tspan, N) # Appel de la méthode d'intégraton\n",
+ "\n",
+ " # Affichage de la solution approchée\n",
+ " plt = plot(ts, xs, label=method_name, marker=:circle)\n",
+ "\n",
+ " # Affichage de la solution analytique\n",
+ " plot!(plt, sol, label=\"Analytic\")\n",
+ "\n",
+ " # Mise en forme du graphique\n",
+ " plot!(plt, xlims=(tspan[1]-0.1, tspan[2]+0.1), ylims=(-0.2, 1.65), \n",
+ " title=\"h=$(round((tspan[2]-tspan[1])/N, digits=2))\", titlefont = font(12,\"Calibri\"),\n",
+ " xlabel=\"t\", ylabel=\"x(t)\", left_margin=15px, bottom_margin=15px, top_margin=10px, right_margin=15px)\n",
+ "\n",
+ " # Ajout du graphique à la liste des graphiques\n",
+ " push!(plts, plt)\n",
+ "\n",
+ " end\n",
+ "\n",
+ " return plts\n",
+ "\n",
+ "end\n",
+ "\n",
+ "xlims_ = (1e-4, 1e0)\n",
+ "ylims_ = (1e-13, 1e1)\n",
+ "xlims_nfe_ = (1e1, 1e4)\n",
+ "ylims_nfe_ = (1e-13, 1e1)\n",
+ "\n",
+ "function ordre(method, f, x0, tspan, sol, hspan, nfespan; \n",
+ " xlims=xlims_, ylims=ylims_, xlims_nfe=xlims_nfe_, ylims_nfe=ylims_nfe_)\n",
+ "\n",
+ " plts = [] # Liste des graphiques\n",
+ "\n",
+ " #\n",
+ " plt1 = plot(xaxis=:log, yaxis=:log); push!(plts, plt1)\n",
+ " plt2 = plot(xaxis=:log, yaxis=:log); push!(plts, plt2)\n",
+ "\n",
+ " #\n",
+ " plts = ordre(plts, method, f, x0, tspan, sol, hspan, nfespan, \n",
+ " xlims=xlims, ylims=ylims, xlims_nfe=xlims_nfe, ylims_nfe=ylims_nfe)\n",
+ "\n",
+ " # Mise en forme des graphiques\n",
+ " plot!(plts[1], titlefont = font(12, \"Calibri\"), legend=:topleft,\n",
+ " xlabel=\"h\", ylabel=\"Error\", left_margin=15px, bottom_margin=15px, top_margin=10px, right_margin=15px)\n",
+ "\n",
+ " plot!(plts[2], titlefont = font(12, \"Calibri\"), legend=:topright,\n",
+ " xlabel=\"Calls to f\", ylabel=\"Error\", left_margin=15px, bottom_margin=15px, top_margin=10px, right_margin=15px)\n",
+ "\n",
+ " return plts\n",
+ " \n",
+ "end\n",
+ "\n",
+ "function ordre(plts_in, method, f, x0, tspan, sol, hspan, nfespan; \n",
+ " xlims=xlims_, ylims=ylims_, xlims_nfe=xlims_nfe_, ylims_nfe=ylims_nfe_)\n",
+ "\n",
+ " # Copie du graphique d'entrée\n",
+ " plts = deepcopy(plts_in)\n",
+ "\n",
+ " # Récupération des informations sur la méthode\n",
+ " method_func, method_name, methode_stages = method_infos(method)\n",
+ " \n",
+ " # Ecriture des choix de paramètres\n",
+ " println(\"Méthode d'intégration : \", method_name)\n",
+ "\n",
+ " # Les différents nombre de pas de temps\n",
+ " Nspan = round.(Int, (tspan[2]-tspan[1]) ./ hspan)\n",
+ "\n",
+ " # Calcul de l'erreur\n",
+ " err = []\n",
+ "\n",
+ " for N ∈ Nspan\n",
+ "\n",
+ " # Solution approchée\n",
+ " ts, xs = method_func(f, x0, tspan, N)\n",
+ "\n",
+ " # On calcule l'erreur en norme infinie\n",
+ " push!(err, maximum(abs.(xs .- sol.(ts))))\n",
+ " \n",
+ " end\n",
+ "\n",
+ " # calcul par régression linéaire de la pente de la droite et de l'ordonnée à l'origine\n",
+ " reg = fit(log10.(hspan), log10.(err), 1)\n",
+ " K = 10^reg[0]\n",
+ " p = reg[1]\n",
+ " println(\"\\nconstante du grand O : K = $(round(K, digits=5))\")\n",
+ " println(\"ordre de convergence : p = $(round(p, digits=5))\")\n",
+ "\n",
+ " # Affichage de l'erreur en fonction du pas de temps: on enlève la constante K\n",
+ " plot!(plts[1], hspan, err ./ K, xaxis=:log, yaxis=:log, label=\"$method_name\", marker=:circle)\n",
+ "\n",
+ " # Affichage de la droite de régression\n",
+ " #plot!(plt, hspan, hspan .^ p, label=\"$method_name Regression\", linestyle=:dash)\n",
+ " \n",
+ " # Mise en forme du graphique\n",
+ " plot!(plts[1], xlims=xlims, ylims=ylims)\n",
+ "\n",
+ " # Calcul de l'erreur en fonction du nombre d'appels à la fonction f\n",
+ " err = []\n",
+ "\n",
+ " for Nfe ∈ Nfespan\n",
+ "\n",
+ " #\n",
+ " N = round(Int, Nfe / methode_stages)\n",
+ "\n",
+ " # Solution approchée\n",
+ " ts, xs = method_func(f, x0, tspan, N)\n",
+ "\n",
+ " # On calcule l'erreur en norme infinie\n",
+ " push!(err, maximum(abs.(xs .- sol.(ts))))\n",
+ " \n",
+ " end\n",
+ "\n",
+ " # Affchage de l'erreur en fonction du nombre d'appels à la fonction f\n",
+ " # Nfespan = methode_stages .* Nspan\n",
+ " plot!(plts[2], Nfespan, err ./ K, xaxis=:log, yaxis=:log, label=\"$method_name\", marker=:circle)\n",
+ "\n",
+ " # Mise en forme du graphique\n",
+ " plot!(plts[2], xlims=xlims_nfe, ylims=ylims_nfe)\n",
+ "\n",
+ " return plts\n",
+ " \n",
+ "end;\n",
+ "\n",
+ "hspan_ = 10 .^ range(-4, stop=0, length=20)\n",
+ "Nfespan_ = 10 .^ range(0, stop=4, length=20);"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## La méthode d'Euler\n",
+ "\n",
+ "On rappelle que la méthode d'Euler est donnée par :\n",
+ "\n",
+ "$$\n",
+ "\\left\\{\n",
+ "\\begin{array}{l}\n",
+ "x_{n+1} = x_n + h f(t_n, x_n) \\\\\n",
+ "x_0 = x(t_0)\n",
+ "\\end{array}\n",
+ "\\right.\n",
+ "$$\n",
+ "\n",
+ "### Exercice 1\n",
+ "\n",
+ "1. Calculer la solution analytique et la coder dans la fonction `sol` ci-dessous.\n",
+ "2. Implémenter la méthode d'Euler dans la fonction `euler` ci-dessous.\n",
+ "3. Vérifier la convergence de la méthode d'Euler sur l'exemple donné dans la cellule d'après."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Solution analytique\n",
+ "function sol(t)\n",
+ " return 0 # TO UPDATE\n",
+ "end;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Implémentation de la méthode d'Euler\n",
+ "function euler(f, x0, tspan, N)\n",
+ " t0, tf = tspan\n",
+ " h = (tf - t0) / N\n",
+ " t = t0\n",
+ " x = x0\n",
+ " ts = [t0]\n",
+ " xs = [x0]\n",
+ " for i in 1:N\n",
+ " x = x # TO UPDATE\n",
+ " t = t # TO UPDATE\n",
+ " push!(ts, t)\n",
+ " push!(xs, x)\n",
+ " end\n",
+ " return ts, xs\n",
+ "end;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Définition du problème de Cauchy\n",
+ "f(t, x) = (1-2t)x # Second membre f(t, x)\n",
+ "x0 = 1.0 # Condition initiale\n",
+ "tspan = (0.0, 3.0); # Intervalle de temps"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Convergence de la méthode sur l'exemple\n",
+ "method = :euler\n",
+ "Nspan = [5, 10, 20]\n",
+ "\n",
+ "# Calcul des graphiques\n",
+ "plts = convergence(method, f, x0, tspan, sol, Nspan)\n",
+ "\n",
+ "# Affichage des graphiques\n",
+ "plot(plts..., layout=(1, length(Nspan)), size=(900, 500))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Exercice 2 \n",
+ "\n",
+ "1. Afficher l'erreur globale de la méthode d'Euler en fonction de $h$ et vérifier que l'erreur est bien en $O(h)$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Ordre de convergence de la méthode\n",
+ "method = :euler\n",
+ "hspan = hspan_[1e-7 .≤ hspan_ .≤ 1]\n",
+ "Nfespan = Nfespan_[10 .≤ Nfespan_ .≤ 1e7]\n",
+ "\n",
+ "# Calcul du graphique\n",
+ "plt_order_euler = ordre(method, f, x0, tspan, sol, hspan, Nfespan)\n",
+ "\n",
+ "# Affichage du graphique\n",
+ "plot(plt_order_euler..., layout=(1, 2), size=(900, 500))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## La méthode de Runge \n",
+ "\n",
+ "On rappelle que la méthode de Runge est donnée par :\n",
+ "\n",
+ "$$\n",
+ "\\left\\{\n",
+ "\\begin{array}{l}\n",
+ "\\displaystyle x_{n+1} = x_n + h\\, f\\left(t_{n} + \\frac{h}{2}, ~x_n + \\frac{h}{2} f(t_n, x_n)\\right) \\\\\n",
+ "x_0 = x(t_0)\n",
+ "\\end{array}\n",
+ "\\right.\n",
+ "$$\n",
+ "\n",
+ "### Exercice 3\n",
+ "\n",
+ "1. Implémenter la méthode de Runge dans la fonction `runge` ci-dessous.\n",
+ "2. Vérifier la convergence de la méthode de Runge sur l'exemple donné dans la cellule d'après.\n",
+ "3. Calculer l'erreur globale de la méthode de Runge en fonction de $h$ et vérifier que l'erreur est bien en $O(h^2)$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "function runge(f, x0, tspan, N)\n",
+ " t0, tf = tspan\n",
+ " h = (tf - t0) / N\n",
+ " t = t0\n",
+ " x = x0\n",
+ " ts = [t0]\n",
+ " xs = [x0]\n",
+ " for i in 1:N\n",
+ " k1 = f(t, x)\n",
+ " k2 = 0 # TO UPDATE\n",
+ " x = x # TO UPDATE\n",
+ " t = t # TO UPDATE\n",
+ " push!(ts, t)\n",
+ " push!(xs, x)\n",
+ " end\n",
+ " return ts, xs\n",
+ "end;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Convergence de la méthode sur l'exemple\n",
+ "method = :runge\n",
+ "Nspan = [5, 10, 20]\n",
+ "\n",
+ "# Calcul des graphiques\n",
+ "plts = convergence(method, f, x0, tspan, sol, Nspan)\n",
+ "\n",
+ "# Affichage des graphiques\n",
+ "plot(plts..., layout=(1, length(Nspan)), size=(900, 500))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Ordre de convergence de la méthode\n",
+ "method = :runge\n",
+ "hspan = hspan_[1e-5 .≤ hspan_ .≤ 1]\n",
+ "Nfespan = Nfespan_[10 .≤ Nfespan_ .≤ 1e6]\n",
+ "\n",
+ "# Calcul du graphique\n",
+ "plt_order_runge = ordre(plt_order_euler, method, f, x0, tspan, sol, hspan, Nfespan)\n",
+ "\n",
+ "# Affichage du graphique\n",
+ "plot(plt_order_runge..., layout=(1, 2), size=(900, 500))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## La méthode de Heun \n",
+ "\n",
+ "On rappelle que la méthode de Heun est donnée par :\n",
+ "\n",
+ "$$\n",
+ "\\left\\{\n",
+ "\\begin{array}{l}\n",
+ " \\displaystyle k_1 = f(t_n, x_n) \\\\[1em]\n",
+ " \\displaystyle k_2 = f(t_n + \\frac{h}{3}, x_n + \\frac{h}{3} k_1) \\\\[1em]\n",
+ " \\displaystyle k_3 = f(t_n + \\frac{2}{3}h, x_n + \\frac{2}{3}h\\, k_2) \\\\[1em]\n",
+ " \\displaystyle x_{n+1} = x_n + h\\left(\\frac{1}{4} k_1 + \\frac{3}{4} k_3\\right), \\quad \n",
+ " \\displaystyle x_0 = x(t_0)\n",
+ "\\end{array}\n",
+ "\\right.\n",
+ "$$\n",
+ "\n",
+ "### Exercice 4\n",
+ "\n",
+ "1. Implémenter la méthode de Heun dans la fonction `heun` ci-dessous.\n",
+ "2. Vérifier la convergence de la méthode de Heun sur l'exemple donné dans la cellule d'après.\n",
+ "3. Calculer l'erreur globale de la méthode de Heun en fonction de $h$ et vérifier que l'erreur est bien en $O(h^3)$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "function heun(f, x0, tspan, N)\n",
+ " t0, tf = tspan\n",
+ " h = (tf - t0) / N\n",
+ " t = t0\n",
+ " x = x0\n",
+ " ts = [t0]\n",
+ " xs = [x0]\n",
+ " for i in 1:N\n",
+ " k1 = f(t, x)\n",
+ " k2 = 0 # TO UPDATE\n",
+ " k3 = 0 # TO UPDATE \n",
+ " x = x # TO UPDATE\n",
+ " t = t # TO UPDATE\n",
+ " push!(ts, t)\n",
+ " push!(xs, x)\n",
+ " end\n",
+ " return ts, xs\n",
+ "end;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Convergence de la méthode sur l'exemple\n",
+ "method = :heun\n",
+ "Nspan = [5, 10, 20]\n",
+ "\n",
+ "# Calcul des graphiques\n",
+ "plts = convergence(method, f, x0, tspan, sol, Nspan)\n",
+ "\n",
+ "# Affichage des graphiques\n",
+ "plot(plts..., layout=(1, length(Nspan)), size=(900, 500))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Ordre de convergence de la méthode\n",
+ "method = :heun\n",
+ "hspan = hspan_[1e-4 .≤ hspan_ .≤ 1]\n",
+ "Nfespan = Nfespan_[10 .≤ Nfespan_ .≤ 1e5]\n",
+ "\n",
+ "# Calcul du graphique\n",
+ "plt_order_heun = ordre(plt_order_runge, method, f, x0, tspan, sol, hspan, Nfespan)\n",
+ "\n",
+ "# Affichage du graphique\n",
+ "plot(plt_order_heun..., layout=(1, 2), size=(900, 500))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## La méthode de Runge-Kutta d'ordre 4\n",
+ "\n",
+ "On rappelle que la méthode de Runge-Kutta d'ordre 4 est donnée par :\n",
+ "\n",
+ "$$\n",
+ "\\left\\{\n",
+ "\\begin{array}{l}\n",
+ " \\displaystyle k_1 = f(t_n, x_n) \\\\[1em]\n",
+ " \\displaystyle k_2 = f(t_n + \\frac{h}{2}, x_n + \\frac{h}{2} k_1) \\\\[1em]\n",
+ " \\displaystyle k_3 = f(t_n + \\frac{h}{2}, x_n + \\frac{h}{2} k_2) \\\\[1em]\n",
+ " \\displaystyle k_4 = f(t_n + h, x_n + h\\, k_3) \\\\[1em]\n",
+ " \\displaystyle x_{n+1} = x_n + \\frac{h}{6} (k_1 + 2 k_2 + 2 k_3 + k_4), \\quad \n",
+ " \\displaystyle x_0 = x(t_0)\n",
+ "\\end{array}\n",
+ "\\right.\n",
+ "$$\n",
+ "\n",
+ "### Exercice 5\n",
+ "\n",
+ "1. Implémenter la méthode de Runge-Kutta d'ordre 4 dans la fonction `rk4` ci-dessous.\n",
+ "2. Vérifier la convergence de la méthode de Runge-Kutta d'ordre 4 sur l'exemple donné dans la cellule d'après.\n",
+ "3. Calculer l'erreur globale de la méthode de Runge-Kutta d'ordre 4 en fonction de $h$ et vérifier que l'erreur est bien en $O(h^4)$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "function rk4(f, x0, tspan, N)\n",
+ " t0, tf = tspan\n",
+ " h = (tf - t0) / N\n",
+ " t = t0\n",
+ " x = x0\n",
+ " ts = [t0]\n",
+ " xs = [x0]\n",
+ " for i in 1:N\n",
+ " k1 = f(t, x)\n",
+ " k2 = 0 # TO UPDATE\n",
+ " k3 = 0 # TO UPDATE\n",
+ " k4 = 0 # TO UPDATE\n",
+ " x = x # TO UPDATE\n",
+ " t = t # TO UPDATE\n",
+ " push!(ts, t)\n",
+ " push!(xs, x)\n",
+ " end\n",
+ " return ts, xs\n",
+ "end;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Convergence de la méthode sur l'exemple\n",
+ "method = :rk4\n",
+ "Nspan = [5, 10, 20]\n",
+ "\n",
+ "# Calcul des graphiques\n",
+ "plts = convergence(method, f, x0, tspan, sol, Nspan)\n",
+ "\n",
+ "# Affichage des graphiques\n",
+ "plot(plts..., layout=(1, length(Nspan)), size=(900, 500))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Ordre de convergence de la méthode\n",
+ "method = :rk4\n",
+ "hspan = hspan_[1e-3 .≤ hspan_ .≤ 1]\n",
+ "Nfespan = Nfespan_[10 .≤ Nfespan_ .≤ 1e4]\n",
+ "\n",
+ "# Calcul du graphique\n",
+ "plt_order_rk4 = ordre(plt_order_heun, method, f, x0, tspan, sol, hspan, Nfespan)\n",
+ "\n",
+ "# Affichage du graphique\n",
+ "plot(plt_order_rk4..., layout=(1, 2), size=(900, 500))"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Julia 1.9.0",
+ "language": "julia",
+ "name": "julia-1.9"
+ },
+ "language_info": {
+ "file_extension": ".jl",
+ "mimetype": "application/julia",
+ "name": "julia",
+ "version": "1.9.0"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}