" plot!(plt, [V[i][j][1] for j ∈ 1:N], [V[i][j][2] for j ∈ 1:N], legend = false, \n",
" color = :green, linewidth = 2)\n",
" end\n",
" return nothing\n",
"end;"
]
...
...
@@ -490,7 +497,7 @@
"plot_traj!(plt, F, z0, tf)\n",
"\n",
"# plot some endpoints of the flow\n",
"z0 = [[-1, p0] for p0 ∈ range(ymin, ymax, length=200)]\n",
"z0 = [[-1, p0] for p0 ∈ range(ymin, ymax, length=50)]\n",
"plot_flow!(plt, F, z0, tf)\n",
"\n",
"# additional plots\n",
...
...
@@ -528,7 +535,7 @@
" α = 1.5\n",
" zpoint = similar(z)\n",
" zpoint[1] = -z[1] + α*z[1]^2 + z[2]\n",
" zpoint[2] = (1. - 2*α*z[1])*z[2]\n",
" zpoint[2] = (1 - 2*α*z[1])*z[2]\n",
" return zpoint\n",
"end\n",
"\n",
...
...
@@ -545,7 +552,7 @@
"plot_traj!(plt, F, z0, tf)\n",
"\n",
"# plot some endpoints of the flow\n",
"z0 = [[-1, p0] for p0 ∈ range(ymin, ymax, length=200)]\n",
"z0 = [[-1, p0] for p0 ∈ range(ymin, ymax, length=100)]\n",
"plot_flow!(plt, F, z0, tf)\n",
"\n",
"# additional plots\n",
...
...
%% Cell type:markdown id: tags:
# 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
