- Durée approximative : 1 séance (à finir à la maison)
%% Cell type:markdown id: tags:
## Introduction
Pour la documentation du package `Julia` sur l'intégration numérique, 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
usingForwardDiff
usingLinearAlgebra
usingPlots
include("utils.jl");# plot_traj!, plot_flow!
```
%% 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.
"""
Second membre de l'IVP
x : vecteur d'état
λ : vecteur de paramètres
t : variable de temps. Ici le temps intervient explicitement, le système est non autonome.
"""
function exemple1(x,λ,t)
xpoint=similar(x)
# A COMPLETER/MODIFIER
xpoint[:]=zeros(size(x))
returnxpoint
end
λ=[]# 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,λ)# 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 : $\mathrm{sc}_i = \mathrm{abstol}_i + \mathrm{reltol}_i \times \max(|x_{0i}|,|x_{1i}|)$. On calcule alors le pas de façon à avoir
Nous allons ici résoudre le problème à valeur initiale $\dot{x}(t) = 1 + x^p(t)$ avec $p \ge 1$ et $x(0)=0$.
- Coder la fonction $f$ qui permet d'écrire l'équation différentielle sous la forme $\dot{x}(t)=f(t,x(t))$.
- Déterminer à partir de quelle valeur de $p\ge 1$ la solution explose en temps fini (à $0.1$ près) inférieur à la valeur 3.
%% Cell type:code id: tags:
``` julia
"""
Second membre de l'IVP
x : vecteur d'état
p : scalaire
t : variable de temps. Ici le temps n'intervient pas explicitement, le système est autonome.
"""
function f(x,p,t)
return0.0
end
#
x0=0.0
t0=0.0
tf=3
tspan=(t0,tf)# instants initial et terminal
#
xspan=0:0.01:2
#
pf=plot()
px=plot()
#
p=1
prob=ODEProblem(f,x0,tspan,p,reltol=1.e-8,abstol=1.e-8)# définition du problème en Julia
sol=solve(prob)# intégration numérique
pf=plot!(pf,xspan,x->(f(x,p,t0)-1),label="p=$p")# affichage du second membre
px=plot!(px,sol,label="p=$p")# affichage de la solution
#
####################################
p=1######### A MODIFIER #########
####################################
prob=ODEProblem(f,x0,tspan,p,reltol=1.e-8,abstol=1.e-8)# définition du problème en Julia
sol=solve(prob)# intégration numérique
pf=plot!(pf,xspan,x->(f(x,p,t0)-1),label="p=$p")# affichage du second membre
px=plot!(px,sol,label="p=$p")# affichage de la solution
#
p=2
prob=ODEProblem(f,x0,tspan,p,reltol=1.e-8,abstol=1.e-8)# définition du problème en Julia
sol=solve(prob)# intégration numérique
pf=plot!(pf,xspan,x->(f(x,p,t0)-1),label="p=$p")# affichage du second membre
px=plot!(px,sol,label="p=$p")# affichage de la solution
#
plot!(pf,xlims=(0,2),ylims=(0,4))
plot!(px,xlims=(t0,tf),ylims=(0,100))
plot(pf,px,size=(900,400))
```
%% Cell type:markdown id: tags:
## Pendule simple
On s'intéresse ici au [pendule simple](https://fr.wikipedia.org/wiki/Pendule_simple). Les principes physiques de la mécanique classique donnent comme équation qui régit l'évolution du mouvement
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 ?
- Remplacer la vitesse angulaire initiale $\dot{\alpha}(0)$ par $2$. Commentaires.
%% Cell type:code id: tags:
``` julia
"""
Second membre de l'IVP
x : vecteur d'état
λ : vecteur de paramètres
t : variable de temps. Ici le temps n'intervient pas explicitement, le système est autonome.
"""
function pendule(x,λ,t)
xpoint=similar(x)
g,l,k,m=λ
# A COMPLETER/MODIFIER
xpoint[:]=zeros(size(x))
returnxpoint
end
#
g=9.81
l=10
k=0
m=1
λ=[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)# instants initial et terminal
prob=ODEProblem(pendule,x0,tspan,λ,reltol=1.e-8,abstol=1.e-8)# définition du problème en Julia
sol=solve(prob)# 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 ci-dessous et interprétez le graphique: où sont les oscillations, les rotations et les points d'équilibre stables et instables ?
- 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
λ=[g,l,k,m]# paramètres constants
plt=plot(xlabel="x₁",ylabel="x₂",legend=false)# initialisation du plot
fortheta0in0:(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
On calcule ci-dessous les valeurs propres de la matrice jacobienne du second membre de l'IVP au point d'équilibre $x_e = (0, 0)$. Commentaires (rappelez-vous votre cours d'[automatique](https://gitlab.irit.fr/toc/etu-n7/controle-optimal/-/raw/master/ressources/notes-autom-2021.pdf)).
%% Cell type:code id: tags:
``` julia
# Jacobienne de la fonction vdp
dvdp(x)=ForwardDiff.jacobian(x->vdp(x,[],0),x)
# Point d'équilibre et la matrice Jacobienne en ce point
xe=[0,0]
A=dvdp(xe)
# Valeurs propres de la matrice Jacobienne
eigvals(A)
```
%% Cell type:markdown id: tags:
## Modèle proie-prédateur (modèle de Lotka-Volterra)
L'équation différentielle considérée est donnée par les [équations de prédation de Lotka-Volterra](https://fr.wikipedia.org/wiki/Équations_de_prédation_de_Lotka-Volterra)
- $x(t)$ est l'effectif des proies en fonction du temps ;
- $y(t)$ est l'effectif des prédateurs en fonction du temps ;
Les paramètres suivants caractérisent les interactions entre les deux espèces :
- $\alpha$, taux de reproduction intrinsèque des proies (constant, indépendant du nombre de prédateurs) ;
- $\beta$, taux de mortalité des proies dû aux prédateurs rencontrés ;
- $\delta$, taux de reproduction des prédateurs en fonction des proies rencontrées et mangées ;
- $\gamma$, taux de mortalité intrinsèque des prédateurs (constant, indépendant du nombre de proies) ;
**Questions.**
- Le point $(0, 0)$ est clairement un point d'équilibre stable. Il existe un autre point d'équilibre, lequel ?
- Afficher l'autre point d'équilibre ainsi que quelques trajectoires du système dans le plan de phase, qui rentrent dans les limites du plot ci-dessous. Qu'observez-vous sur les trajectoires ?
%% Cell type:code id: tags:
``` julia
# ------------------------------
# DEBUT A COMPLETER/MODIFIER LE SECOND MEMBRE DE L'IVP
# ...
# FIN A COMPLETER/MODIFIER LE SECOND MEMBRE DE L'IVP
# ------------------------------
# paramètres
α=2/3
β=4/3
γ=1
δ=1
λ=[α,β,γ,δ]
t0=0
tf=20
tspan=(t0,tf)
plt=plot(xlabel="x₁",ylabel="x₂",legend=false)
# ------------------------------
# DEBUT A COMPLETER/MODIFIER L'AFFICHAGE
# ...
# FIN A COMPLETER/MODIFIER L'AFFICHAGE
# ------------------------------
plot!(xlims=(0,4),ylims=(0,2.5),size=(850,550))
```
%% Cell type:markdown id: tags:
## Phénomène de résonance
Soit $\omega$ et $\omega_0$ deux réels strictement positifs. On considère l'équation différentielle
- En prenant comme variable d'état $x=(x_1,x_2)=(y, \dot{y})$, é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 : les solutions sont-elles bornées ?
- Remplacer la vitesse angulaire pour avoir $\omega = \omega_0$. Commentaires.
Remarque : voir la page Wikipedia sur le [phénomène de résonance](https://fr.wikipedia.org/wiki/Résonance#:~:text=Le%20phénomène%20de%20résonance%20n,en%20phase%20»%20avec%20cette%20dernière.) pour plus d'informations.
%% Cell type:code id: tags:
``` julia
"""
Second membre de l'IVP
x : vecteur d'état
λ : vecteur de paramètres
t : variable de temps. Ici le temps intervient explicitement, le système est non autonome.
