"sol = solve(prob) # résolution du problème IVP\n",
"p1 = plot(sol, label = [\"x₁(t)\" \"x₂(t)\"], title = \"default options\") # affichage de la solution"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Le contrôle du pas\n",
"Les bons intégrateurs numériques ont une mise à jour automatique du pas. Pour les \n",
"[méthodes de Runge-Kutta](https://fr.wikipedia.org/wiki/Méthodes_de_Runge-Kutta) \n",
"par exemple, sur un pas $[t_0, t_{1}]$ on calcule par 2 schémas différents deux solutions approchées à l'instant $t_{1}$ : $x_{1}$ et $\\hat{x}_{1}$ dont l'erreur locale est respectivement en $O(h^p)$ et $O(h^{p+1})$ (les erreurs globales sont elles en $O(h^{p-1})$ et $O(h^p)$). Ainsi on peut estimer l'erreur locale du schéma le moins précis à l'aide de la différence $x_{1}-\\hat{x}_{1}$. On peut alors estimer le pas de façon à avoir la norme de cette différence inférieure à une tolérance donnée `Tol` ou encore à avoir \n",
"En pratique on considère des tolérences absolue et relative composante par composante et on définit : $sc_i = \\mathrm{Atol}_i + \\max(|x_{0i}|,|x_{1i}|)\\, \\mathrm{Rtol}_i$. On calcule alors le pas de façon à avoir \n",
"où $\\ddot{\\alpha}(t)$ désigne la dérivée seconde de l'angle $\\alpha$ par rapport au temps $t$. \n",
"\n",
"- En prenant comme variable d'état $x=(x_1,x_2)=(\\alpha, \\dot{\\alpha})$, écrire la fonction $f$ qui permet d'écrire l'équation différentielle sous la forme $\\dot{x}(t) = f(t,x(t))$\n",
"- Coder cette fonction ci-dessous et exécuter le code. D'après-vous, observez-vous une oscillation ou une rotation du pendule ?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"function pendule(x,p,t)\n",
" # Second membre de l'IVP\n",
" # p : vecteur de paramètres\n",
" # t : variable de temps. Ici le temps n'intervient pas explicitement.\n",
" xpoint = similar(x)\n",
" \n",
" # A COMPLETER/MODIFIER\n",
" xpoint = zeros(size(x))\n",
" \n",
" return xpoint \n",
"end\n",
"\n",
"#\n",
"g = 9.81; l = 10; k = 0; m = 1;\n",
"p = [g, l, k, m] # paramètres constants\n",
"\n",
"theta0 = pi/3\n",
"x0 = [theta0, 0] # état initial\n",
"\n",
"t0 = 0.\n",
"tf = 3*pi*sqrt(l/g)*(1 + theta0^2/16 + theta0^4/3072) # 2*approximation de la période\n",
"tspan = (t0, tf) # instant initial et terminal\n",
"\n",
"prob = ODEProblem(pendule,x0,tspan,p) # défini le problème en Julia\n",
# Intégration numérique avec Julia, une introduction
%% Cell type:markdown id: tags:
## Introduction
Pour la documentation voir
[documentation de DifferentialEquations.jl](https://diffeq.sciml.ai/stable/)
Il nous faut tout d'abord importer les bons packages de Julia.
%% Cell type:code id: tags:
``` julia
usingDifferentialEquations
usingPlots
```
%% Cell type:markdown id: tags:
## Exemple 1
Nous allons ici résoudre le problème à valeur initiale
$$(IVP)\left\{\begin{array}{l}
\dot{x}_1(t)=x_1(t)+x_2(t)+\sin t\\
\dot{x}_2(t)=-x_1(t)+3x_2(t)\\
x_1(0)=-9/25\\
x_2(0)=-4/25,
\end{array}\right.
$$
- Ecrire la fonction $f$ qui permet d'écrire l'équation différentielle sous la forme $\dot{x}(t)=f(t,x(t))$.
- Vérifier que la solution de $(IVP$ est
$$\begin{align*}
x_1(t)&= (-1/25)(13\sin t+9\cos t)\\
x_2(t)&= (-1/25)(3\sin t+4\cos t).
\end{align*}$$
- Coder la fonction $f$ et exécuter le code ci-dessous. Commentaires ?
%% Cell type:code id: tags:
``` julia
# Remarque :
# un problème IVP est une équation différentielle ordinaire
# (EDO en français, ODE en anglais)
# avec en plus une condition initiale.
function exemple1(x,p,t)
# Second membre de l'IVP
# p : vecteur de paramètres
# t : variable de temps
xpoint=similar(x)
# A COMPLETER/MODIFIER
xpoint=zeros(size(x))
returnxpoint
end
p=[]# pas de paramètres
t0=0
tf=10
tspan=(t0,tf)# intervalle d'intégration
x0=[-9/25,-4/25]# condition initiale
prob=ODEProblem(exemple1,x0,tspan,p)# définition du problème IVP
sol=solve(prob)# résolution du problème IVP
p1=plot(sol,label=["x₁(t)""x₂(t)"],title="default options")# affichage de la solution
```
%% Cell type:markdown id: tags:
## Le contrôle du pas
Les bons intégrateurs numériques ont une mise à jour automatique du pas. Pour les
[méthodes de Runge-Kutta](https://fr.wikipedia.org/wiki/Méthodes_de_Runge-Kutta)
par exemple, sur un pas $[t_0, t_{1}]$ on calcule par 2 schémas différents deux solutions approchées à l'instant $t_{1}$ : $x_{1}$ et $\hat{x}_{1}$ dont l'erreur locale est respectivement en $O(h^p)$ et $O(h^{p+1})$ (les erreurs globales sont elles en $O(h^{p-1})$ et $O(h^p)$). Ainsi on peut estimer l'erreur locale du schéma le moins précis à l'aide de la différence $x_{1}-\hat{x}_{1}$. On peut alors estimer le pas de façon à avoir la norme de cette différence inférieure à une tolérance donnée `Tol` ou encore à avoir
$$\frac{\|x_{1}-\hat{x}_{1}\|}{\mathrm{Tol}}<1.$$
En pratique on considère des tolérences absolue et relative composante par composante et on définit : $sc_i = \mathrm{Atol}_i + \max(|x_{0i}|,|x_{1i}|)\,\mathrm{Rtol}_i$. On calcule alors le pas de façon à avoir
où $\ddot{\alpha}(t)$ désigne la dérivée seconde de l'angle $\alpha$ par rapport au temps $t$.
- En prenant comme variable d'état $x=(x_1,x_2)=(\alpha, \dot{\alpha})$, écrire la fonction $f$ qui permet d'écrire l'équation différentielle sous la forme $\dot{x}(t) = f(t,x(t))$
- Coder cette fonction ci-dessous et exécuter le code. D'après-vous, observez-vous une oscillation ou une rotation du pendule ?
%% Cell type:code id: tags:
``` julia
```
%%Celltype:codeid:tags:
``` julia
function pendule(x,p,t)
# Second membre de l'IVP
# p : vecteur de paramètres
# t : variable de temps. Ici le temps n'intervient pas explicitement.
xpoint = similar(x)
# A COMPLETER/MODIFIER
xpoint = zeros(size(x))
return xpoint
end
#
g = 9.81; l = 10; k = 0; m = 1;
p = [g, l, k, m] # paramètres constants
theta0 = pi/3
x0 = [theta0, 0] # état initial
t0 = 0.
tf = 3*pi*sqrt(l/g)*(1 + theta0^2/16 + theta0^4/3072) # 2*approximation de la période
tspan = (t0, tf) # instant initial et terminal
prob = ODEProblem(pendule,x0,tspan,p) # défini le problème en Julia
sol = solve(prob) # réalise l'intégration numérique
plot(sol, label = ["x₁(t)" "x₂(t)"]) # affichage de la solution
```
%% Cell type:markdown id: tags:
## Diagramme de phases
- Exécutez le code et interprétez le graphique.
- On considère le cas où $k=0.15$ (on introduit un frottement). Que se passe-t-il ?
%% Cell type:code id: tags:
``` julia
#
g = 9.81; l = 1.5; k = 0; m = 1;
p = [g, l, k, m] # paramètres constants
plot()
for theta0 in 0:(2*pi)/10:2*pi
theta0_princ = theta0
tf = 3*pi*sqrt(l/g)*(1 + theta0_princ^2/16 + theta0_princ^4/3072) # 2*approximation of the period
