{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[<img src=\"https://gitlab.irit.fr/toc/etu-n7/controle-optimal/-/raw/master/ressources/Logo-toulouse-inp-N7.png\" alt=\"N7\" height=\"80\"/>](https://gitlab.irit.fr/toc/etu-n7/controle-optimal)\n",
"<img src=\"https://gitlab.irit.fr/toc/ens-n7/texCoursN7/-/raw/main/logo-insa.png\" alt=\"INSA\" height=\"80\" style=\"margin-left:50px\"/>\n",
"\n",
"# Méthodes de Runge-Kutta\n",
"\n",
"- Date : 2024-2025\n",
"- Durée approximative : 1h15\n",
"\n",
"<div style=\"width:95%;\n",
" margin:10px;\n",
" padding:8px;\n",
" background-color:#afa;\n",
" border:2px solid #bbffbb;\n",
" border-radius:10px;\n",
" font-weight:bold;\n",
" font-size:1.5em;\n",
" text-align:center;\">\n",
"Introduction\n",
"</div>\n",
"\n",
"**Méthode à un pas explicite.** On désire calculer une approximation de la solution, notée $x(\\cdot, t_0, x_0)$, sur l'intervalle \n",
"$I := [t_0, t_f]$ du problème de Cauchy\n",
"\\begin{equation*}\n",
" \\dot{x}(t) = f(t, x(t)), \\quad x(t_0) = x_0,\n",
"\\end{equation*}\n",
"où $f \\colon \\mathbb{R} \\times \\mathbb{R}^n \\to \\mathbb{R}^n$ est une application suffisamment régulière. \n",
"On considère pour cela une subdivision de l'intervalle $I$ : \n",
"$$\n",
" t_0 < t_1 < \\cdots < t_N := t_f\n",
"$$\n",
"\n",
"On note $h_i := t_{i+1} - t_i$ pour $i \\in \\llbracket 0, N-1 \\rrbracket$, les **pas** de la subdivision et \n",
"$h_{\\max} := \\max_i(h_i)$ le pas le plus long. L'idée est de calculer une approximation de \n",
"la solution en les points de discrétisation. On souhaite donc approcher $x(t_i, t_0, x_0)$ pour \n",
"$i \\in \\llbracket 1, N \\rrbracket$. \n",
"\n",
"On appelle *méthode à un pas explicite*, toute méthode calculant un N-uplet $(x_1, \\ldots, x_N) \\in {(\\mathbb{R}^n)}^N$ tel que la valeur $x_{i+1}$ est calculée en fonction de $t_i$, $h_i$ et $x_i$ de la forme :\n",
"\n",
"$$\n",
" x_{i+1} = x_i + h_i\\, \\Phi(t_i, x_i, h_i).\n",
"$$\n",
"\n",
"Pour approcher la solution, nous utiliserons une méthode à un pas. Si on note $x_h^* := (x_1^*, \\ldots, x_N^*)$ une solution au problème discrétisé (obtenue par la méthode à un pas considérée), alors on souhaite que \n",
"$$\n",
" x_i^* \\approx x(t_i, t_0, x_0) \\quad \\text{pour tout } i \\in \\llbracket 1, N \\rrbracket.\n",
"$$\n",
"\n",
"**Ordre de convergence.** On considère à partir de maintenant un pas de temps $h$ **uniforme**. Une méthode à un pas explicite est **convergente** si pour toute solution $x(\\cdot, t_0, x_0)$, la suite ${(x_i)}_i$ définie par $x_{i+1} = x_i + h\\, \\Phi(t_i, x_i, h)$ vérifie \n",
"$$\n",
" \\max_{1 \\le i \\le N}\\, \\|{x(t_i, t_0, x_0) - x_i}\\| \\to 0 \n",
" \\quad\\text{quand}\\quad h \\to 0.\n",
"$$\n",
"\n",
"Si la convergence est d'ordre $p$ alors il existe une constante $C \\ge 0$ telle que\n",
"\n",
"$$\n",
" E := \\max_{1 \\le i \\le N}\\, \\|{x(t_i, t_0, x_0) - x_i}\\| \\le C\\, h^p.\n",
"$$\n",
"\n",
"Faisons l'hypothèse que $E = M\\, h^p$ pour un certain $M \\ge 0$. En passant au logarithme, on obtient\n",
"\n",
"$$\n",
" \\log(E) = \\log(M) + p\\, \\log(h).\n",
"$$\n",
"\n",
"Nous en déduisons que si on trace $\\log(E)$ en fonction de $\\log(h)$, on doit obtenir une droite de pente $p$. C'est ce que nous allons vérifier dans ce TP.\n",
"\n",
"**Objectif.** Nous allons dans ce TP, implémenter quelques méthodes (à un pas) de [Runge-Kutta](https://fr.wikipedia.org/wiki/Méthodes_de_Runge-Kutta) explicites et \n",
"étudier leur convergence. On appelle *méthode de Runge-Kutta explicite à $s$ étages*, la méthode définie par le schéma\n",
"\n",
"$$\n",
" \\begin{aligned}\n",
" k_1 &= f(t_0, x_0) \\\\\n",
" k_2 &= f(t_0 + c_2 h, x_0 + h a_{21} k_1) \\\\\n",
" \\vdots & \\\\[-1em]\n",
" k_s &= f(t_0 + c_s h, x_0 + h \\sum_{j=1}^{s-1} a_{sj} k_j) \\\\\n",
" x_1 &= x_0 + h \\sum_{i=1}^{s} b_i k_i,\n",
" \\end{aligned}\n",
"$$\n",
"où les coefficients $c_i$, $a_{ij}$ et $b_i$ sont des constantes qui définissent précisément le schéma. \n",
"On supposera toujours dans la suite que $c_1 = 0$ et $c_i = \\sum_{j=1}^{i-1} a_{ij}$, pour $i = 2, \\ldots, s$."
]
},
{
"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, 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, 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": [
"<div style=\"width:95%;\n",
" margin:10px;\n",
" padding:8px;\n",
" background-color:#afa;\n",
" border:2px solid #bbffbb;\n",
" border-radius:10px;\n",
" font-weight:bold;\n",
" font-size:1.5em;\n",
" text-align:center;\">\n",
"L'exemple d'étude\n",
"</div>\n",
"\n",
"On s'intéresse dans ce TP au problème de Cauchy\n",
"\n",
"$$\n",
" \\dot{x}(t) = (1-2t) x(t), \\quad x(0) = x_0 = 1\n",
"$$\n",
"\n",
"sur l'intervalle $[0, 3]$."
]
},
{
"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": "markdown",
"metadata": {},
"source": [
"### Exercice 1\n",
"\n",
"* Calculer la solution analytique de l'exemple ci-dessus et la coder dans la fonction `sol` ci-dessous."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Solution analytique\n",
"function sol(t)\n",
" return 0 # TO UPDATE\n",
"end;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<div style=\"width:95%;\n",
" margin:10px;\n",
" padding:8px;\n",
" background-color:#afa;\n",
" border:2px solid #bbffbb;\n",
" border-radius:10px;\n",
" font-weight:bold;\n",
" font-size:1.5em;\n",
" text-align:center;\">\n",
"La méthode d'Euler\n",
"</div>\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. Implémenter la méthode d'Euler dans la fonction `euler` ci-dessous.\n",
"2. 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": [
"# 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": [
"# 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=(800, 400))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"3. 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=(800, 400))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<div style=\"width:95%;\n",
" margin:10px;\n",
" padding:8px;\n",
" background-color:#afa;\n",
" border:2px solid #bbffbb;\n",
" border-radius:10px;\n",
" font-weight:bold;\n",
" font-size:1.5em;\n",
" text-align:center;\">\n",
"La méthode de Runge \n",
"</div>\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=(800, 400))"
]
},
{
"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=(800, 400))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<div style=\"width:95%;\n",
" margin:10px;\n",
" padding:8px;\n",
" background-color:#afa;\n",
" border:2px solid #bbffbb;\n",
" border-radius:10px;\n",
" font-weight:bold;\n",
" font-size:1.5em;\n",
" text-align:center;\">\n",
"La méthode de Heun \n",
"</div>\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=(800, 400))"
]
},
{
"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=(800, 400))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<div style=\"width:95%;\n",
" margin:10px;\n",
" padding:8px;\n",
" background-color:#afa;\n",
" border:2px solid #bbffbb;\n",
" border-radius:10px;\n",
" font-weight:bold;\n",
" font-size:1.5em;\n",
" text-align:center;\">\n",
"La méthode de Runge-Kutta d'ordre 4\n",
"</div>\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=(800, 400))"
]
},
{
"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=(800, 400))"
]
}
],
"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
}