Skip to content
Snippets Groups Projects
Commit f0c568e2 authored by Olivier Cots's avatar Olivier Cots
Browse files

add RK TP

parent 8dd38eef
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 607 additions and 0 deletions
tp/runge-kutta.ipynb 0 → 100644
+ 607
− 0
View file @ f0c568e2
%% Cell type:markdown id: tags:
<img src="https://gitlab.irit.fr/toc/ens-n7/texCoursN7/-/raw/main/LOGO_INP_N7.png" alt="N7" height="80"/>
<img src="https://gitlab.irit.fr/toc/ens-n7/texCoursN7/-/raw/main/logo-insa.png" alt="INSA" height="80"/>
# Méthodes de Runge-Kutta
- Date : 2023-2024
- Durée approximative : 1h15
## Introduction
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 id: tags:
``` julia
# activate local project
using Pkg
Pkg.activate(".")
# load packages
using Plots
using Plots.PlotMeasures
using Polynomials
#
px = PlotMeasures.px;
```
%% Cell type:code id: tags:
``` julia
# Fonctons auxiliaires
function method_infos(method)
if method == :euler
method_func = euler
method_name = "Euler"
methode_stages = 1
elseif method == :runge
method_func = runge
method_name = "Runge"
methode_stages = 2
elseif method == :heun
method_func = heun
method_name = "Heun"
methode_stages = 3
elseif method == :rk4
method_func = rk4
method_name = "RK4"
methode_stages = 4
else
error("Méthode d'intégration non reconnue")
end
return method_func, method_name, methode_stages
end
function convergence(method, f, x0, tspan, sol, Nspan)
# Récupération des informations sur la méthode
method_func, method_name, methode_stages = method_infos(method)
# Ecriture des choix de paramètres
println("Méthode d'intégration : ", method_name)
plts = [] # Liste des graphiques
for N ∈ Nspan
# Solution approchée
ts, xs = method_func(f, x0, tspan, N) # Appel de la méthode d'intégraton
# Affichage de la solution approchée
plt = plot(ts, xs, label=method_name, marker=:circle)
# Affichage de la solution analytique
plot!(plt, sol, label="Analytic")
# Mise en forme du graphique
plot!(plt, xlims=(tspan[1]-0.1, tspan[2]+0.1), ylims=(-0.2, 1.65),
title="h=$(round((tspan[2]-tspan[1])/N, digits=2))", titlefont = font(12,"Calibri"),
xlabel="t", ylabel="x(t)", left_margin=15px, bottom_margin=15px, top_margin=10px, right_margin=15px)
# Ajout du graphique à la liste des graphiques
push!(plts, plt)
end
return plts
end
xlims_ = (1e-4, 1e0)
ylims_ = (1e-13, 1e1)
xlims_nfe_ = (1e1, 1e4)
ylims_nfe_ = (1e-13, 1e1)
function ordre(method, f, x0, tspan, sol, hspan, nfespan;
xlims=xlims_, ylims=ylims_, xlims_nfe=xlims_nfe_, ylims_nfe=ylims_nfe_)
plts = [] # Liste des graphiques
#
plt1 = plot(xaxis=:log, yaxis=:log); push!(plts, plt1)
plt2 = plot(xaxis=:log, yaxis=:log); push!(plts, plt2)
#
plts = ordre(plts, method, f, x0, tspan, sol, hspan, nfespan,
xlims=xlims, ylims=ylims, xlims_nfe=xlims_nfe, ylims_nfe=ylims_nfe)
# Mise en forme des graphiques
plot!(plts[1], titlefont = font(12, "Calibri"), legend=:topleft,
xlabel="h", ylabel="Error", left_margin=15px, bottom_margin=15px, top_margin=10px, right_margin=15px)
plot!(plts[2], titlefont = font(12, "Calibri"), legend=:topright,
xlabel="Calls to f", ylabel="Error", left_margin=15px, bottom_margin=15px, top_margin=10px, right_margin=15px)
return plts
end
function ordre(plts_in, method, f, x0, tspan, sol, hspan, nfespan;
xlims=xlims_, ylims=ylims_, xlims_nfe=xlims_nfe_, ylims_nfe=ylims_nfe_)
# Copie du graphique d'entrée
plts = deepcopy(plts_in)
# Récupération des informations sur la méthode
method_func, method_name, methode_stages = method_infos(method)
# Ecriture des choix de paramètres
println("Méthode d'intégration : ", method_name)
# Les différents nombre de pas de temps
Nspan = round.(Int, (tspan[2]-tspan[1]) ./ hspan)
# Calcul de l'erreur
err = []
for N ∈ Nspan
# Solution approchée
ts, xs = method_func(f, x0, tspan, N)
# On calcule l'erreur en norme infinie
push!(err, maximum(abs.(xs .- sol.(ts))))
end
# calcul par régression linéaire de la pente de la droite et de l'ordonnée à l'origine
reg = fit(log10.(hspan), log10.(err), 1)
K = 10^reg[0]
p = reg[1]
println("\nconstante du grand O : K = $(round(K, digits=5))")
println("ordre de convergence : p = $(round(p, digits=5))")
# Affichage de l'erreur en fonction du pas de temps: on enlève la constante K
plot!(plts[1], hspan, err ./ K, xaxis=:log, yaxis=:log, label="$method_name", marker=:circle)
# Affichage de la droite de régression
#plot!(plt, hspan, hspan .^ p, label="$method_name Regression", linestyle=:dash)
# Mise en forme du graphique
plot!(plts[1], xlims=xlims, ylims=ylims)
# Calcul de l'erreur en fonction du nombre d'appels à la fonction f
err = []
for Nfe ∈ Nfespan
#
N = round(Int, Nfe / methode_stages)
# Solution approchée
ts, xs = method_func(f, x0, tspan, N)
# On calcule l'erreur en norme infinie
push!(err, maximum(abs.(xs .- sol.(ts))))
end
# Affchage de l'erreur en fonction du nombre d'appels à la fonction f
# Nfespan = methode_stages .* Nspan
plot!(plts[2], Nfespan, err ./ K, xaxis=:log, yaxis=:log, label="$method_name", marker=:circle)
# Mise en forme du graphique
plot!(plts[2], xlims=xlims_nfe, ylims=ylims_nfe)
return plts
end;
hspan_ = 10 .^ range(-4, stop=0, length=20)
Nfespan_ = 10 .^ range(0, stop=4, length=20);
```
%% Cell type:markdown id: tags:
## La méthode d'Euler
On rappelle que la méthode d'Euler est donnée par :
$$
\left\{
\begin{array}{l}
x_{n+1} = x_n + h f(t_n, x_n) \\
x_0 = x(t_0)
\end{array}
\right.
$$
### Exercice 1
1. Calculer la solution analytique et la coder dans la fonction `sol` ci-dessous.
2. Implémenter la méthode d'Euler dans la fonction `euler` ci-dessous.
3. Vérifier la convergence de la méthode d'Euler sur l'exemple donné dans la cellule d'après.
%% Cell type:code id: tags:
``` julia
# Solution analytique
function sol(t)
return 0 # TO UPDATE
end;
```
%% Cell type:code id: tags:
``` julia
# Implémentation de la méthode d'Euler
function euler(f, x0, tspan, N)
t0, tf = tspan
h = (tf - t0) / N
t = t0
x = x0
ts = [t0]
xs = [x0]
for i in 1:N
x = x # TO UPDATE
t = t # TO UPDATE
push!(ts, t)
push!(xs, x)
end
return ts, xs
end;
```
%% Cell type:code id: tags:
``` julia
# Définition du problème de Cauchy
f(t, x) = (1-2t)x # Second membre f(t, x)
x0 = 1.0 # Condition initiale
tspan = (0.0, 3.0); # Intervalle de temps
```
%% Cell type:code id: tags:
``` julia
# Convergence de la méthode sur l'exemple
method = :euler
Nspan = [5, 10, 20]
# Calcul des graphiques
plts = convergence(method, f, x0, tspan, sol, Nspan)
# Affichage des graphiques
plot(plts..., layout=(1, length(Nspan)), size=(900, 500))
```
%% Cell type:markdown id: tags:
### Exercice 2
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 id: tags:
``` julia
# Ordre de convergence de la méthode
method = :euler
hspan = hspan_[1e-7 .≤ hspan_ .≤ 1]
Nfespan = Nfespan_[10 .≤ Nfespan_ .≤ 1e7]
# Calcul du graphique
plt_order_euler = ordre(method, f, x0, tspan, sol, hspan, Nfespan)
# Affichage du graphique
plot(plt_order_euler..., layout=(1, 2), size=(900, 500))
```
%% Cell type:markdown id: tags:
## La méthode de Runge
On rappelle que la méthode de Runge est donnée par :
$$
\left\{
\begin{array}{l}
\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) \\
x_0 = x(t_0)
\end{array}
\right.
$$
### Exercice 3
1. Implémenter la méthode de Runge dans la fonction `runge` ci-dessous.
2. Vérifier la convergence de la méthode de Runge sur l'exemple donné dans la cellule d'après.
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 id: tags:
``` julia
function runge(f, x0, tspan, N)
t0, tf = tspan
h = (tf - t0) / N
t = t0
x = x0
ts = [t0]
xs = [x0]
for i in 1:N
k1 = f(t, x)
k2 = 0 # TO UPDATE
x = x # TO UPDATE
t = t # TO UPDATE
push!(ts, t)
push!(xs, x)
end
return ts, xs
end;
```
%% Cell type:code id: tags:
``` julia
# Convergence de la méthode sur l'exemple
method = :runge
Nspan = [5, 10, 20]
# Calcul des graphiques
plts = convergence(method, f, x0, tspan, sol, Nspan)
# Affichage des graphiques
plot(plts..., layout=(1, length(Nspan)), size=(900, 500))
```
%% Cell type:code id: tags:
``` julia
# Ordre de convergence de la méthode
method = :runge
hspan = hspan_[1e-5 .≤ hspan_ .≤ 1]
Nfespan = Nfespan_[10 .≤ Nfespan_ .≤ 1e6]
# Calcul du graphique
plt_order_runge = ordre(plt_order_euler, method, f, x0, tspan, sol, hspan, Nfespan)
# Affichage du graphique
plot(plt_order_runge..., layout=(1, 2), size=(900, 500))
```
%% Cell type:markdown id: tags:
## La méthode de Heun
On rappelle que la méthode de Heun est donnée par :
$$
\left\{
\begin{array}{l}
\displaystyle k_1 = f(t_n, x_n) \\[1em]
\displaystyle k_2 = f(t_n + \frac{h}{3}, x_n + \frac{h}{3} k_1) \\[1em]
\displaystyle k_3 = f(t_n + \frac{2}{3}h, x_n + \frac{2}{3}h\, k_2) \\[1em]
\displaystyle x_{n+1} = x_n + h\left(\frac{1}{4} k_1 + \frac{3}{4} k_3\right), \quad
\displaystyle x_0 = x(t_0)
\end{array}
\right.
$$
### Exercice 4
1. Implémenter la méthode de Heun dans la fonction `heun` ci-dessous.
2. Vérifier la convergence de la méthode de Heun sur l'exemple donné dans la cellule d'après.
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 id: tags:
``` julia
function heun(f, x0, tspan, N)
t0, tf = tspan
h = (tf - t0) / N
t = t0
x = x0
ts = [t0]
xs = [x0]
for i in 1:N
k1 = f(t, x)
k2 = 0 # TO UPDATE
k3 = 0 # TO UPDATE
x = x # TO UPDATE
t = t # TO UPDATE
push!(ts, t)
push!(xs, x)
end
return ts, xs
end;
```
%% Cell type:code id: tags:
``` julia
# Convergence de la méthode sur l'exemple
method = :heun
Nspan = [5, 10, 20]
# Calcul des graphiques
plts = convergence(method, f, x0, tspan, sol, Nspan)
# Affichage des graphiques
plot(plts..., layout=(1, length(Nspan)), size=(900, 500))
```
%% Cell type:code id: tags:
``` julia
# Ordre de convergence de la méthode
method = :heun
hspan = hspan_[1e-4 .≤ hspan_ .≤ 1]
Nfespan = Nfespan_[10 .≤ Nfespan_ .≤ 1e5]
# Calcul du graphique
plt_order_heun = ordre(plt_order_runge, method, f, x0, tspan, sol, hspan, Nfespan)
# Affichage du graphique
plot(plt_order_heun..., layout=(1, 2), size=(900, 500))
```
%% Cell type:markdown id: tags:
## La méthode de Runge-Kutta d'ordre 4
On rappelle que la méthode de Runge-Kutta d'ordre 4 est donnée par :
$$
\left\{
\begin{array}{l}
\displaystyle k_1 = f(t_n, x_n) \\[1em]
\displaystyle k_2 = f(t_n + \frac{h}{2}, x_n + \frac{h}{2} k_1) \\[1em]
\displaystyle k_3 = f(t_n + \frac{h}{2}, x_n + \frac{h}{2} k_2) \\[1em]
\displaystyle k_4 = f(t_n + h, x_n + h\, k_3) \\[1em]
\displaystyle x_{n+1} = x_n + \frac{h}{6} (k_1 + 2 k_2 + 2 k_3 + k_4), \quad
\displaystyle x_0 = x(t_0)
\end{array}
\right.
$$
### Exercice 5
1. Implémenter la méthode de Runge-Kutta d'ordre 4 dans la fonction `rk4` ci-dessous.
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.
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 id: tags:
``` julia
function rk4(f, x0, tspan, N)
t0, tf = tspan
h = (tf - t0) / N
t = t0
x = x0
ts = [t0]
xs = [x0]
for i in 1:N
k1 = f(t, x)
k2 = 0 # TO UPDATE
k3 = 0 # TO UPDATE
k4 = 0 # TO UPDATE
x = x # TO UPDATE
t = t # TO UPDATE
push!(ts, t)
push!(xs, x)
end
return ts, xs
end;
```
%% Cell type:code id: tags:
``` julia
# Convergence de la méthode sur l'exemple
method = :rk4
Nspan = [5, 10, 20]
# Calcul des graphiques
plts = convergence(method, f, x0, tspan, sol, Nspan)
# Affichage des graphiques
plot(plts..., layout=(1, length(Nspan)), size=(900, 500))
```
%% Cell type:code id: tags:
``` julia
# Ordre de convergence de la méthode
method = :rk4
hspan = hspan_[1e-3 .≤ hspan_ .≤ 1]
Nfespan = Nfespan_[10 .≤ Nfespan_ .≤ 1e4]
# Calcul du graphique
plt_order_rk4 = ordre(plt_order_heun, method, f, x0, tspan, sol, hspan, Nfespan)
# Affichage du graphique
plot(plt_order_rk4..., layout=(1, 2), size=(900, 500))
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment