« Découvrir Scilab/Calcul différentiel et intégral » : différence entre les versions
Contenu supprimé Contenu ajouté
→Produit de convolution : EDO (import de Résolution d'équations) |
|||
Ligne 169 :
== Équation différentielle ordinaire ==
=== Équation différentielle ordinaire de premier ordre ===
Une d'équation différentielle ordinaire ''(ordinary differential equation)'', en abrégé ÉDO ''(ODE)'', peut être résolue de manière numérique avec la fonction <font id="ode1"><code>ode</code></span>. Soit une fonction ''y''(''t'') ; si l'équation différentielle est de premier ordre
Ligne 199 ⟶ 201 :
; Note
: La fonction ƒ doit être continue et différentiable en tout point. Si la fonction présente des singularité, alors il faut la résoudre sur la première partie D<sup>1</sup> du domaine, puis sur la seconde partie en prenant comme conditions limites les résultats du premier calcul.
=== Équation différentielle ordinaire d'ordre plus élevé ===
Une ÉDO d'ordre supérieur peut se réduire à une ÉDO d'ordre 1.
{{loupe|w:fr:Équation différentielle ordinaire#Réduction à 1 de l'ordre d'une équation}}
== Équations différentielles ordinaires d'ordres plus élevés ==
Une ÉDO d'ordre supérieur peut se réduire à une ÉDO d'ordre 1.
{{loupe|w:fr:Équation différentielle ordinaire#Réduction à 1 de l'ordre d'une équation}}
Nous pouvons donc utiliser la même fonction grâce à un changement de variables. Considérons par exemple l'ÉDO linéaire à coefficients constants d'ordre 2 suivante :
: ''y'' ’’ + ''ay'' ’ + ''by'' + c = 0.
En posant
: ''x''<sub>1</sub> = ''y''
: ''x''<sub>2</sub> = ''y'' ’
l'équation s'écrit
: ''x''<sub>2</sub>’ + ''ax''<sub>1</sub>’ + ''bx''<sub>1</sub> + ''c'' = 0.
En posant le vecteur
: <math>\mathrm{Y} = \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}</math>
on a
: <math>\mathrm{Y}' = \begin{pmatrix} x_1' \\ x_2' \end{pmatrix} = \begin{pmatrix} x_2 \\ -a x_1' - bx_1 - c \end{pmatrix}.</math>
En physique, on s'intéresse souvent à des fonctions du temps et l'on parle d'oscillations harmoniques ; on écrit souvent cette équation sous la forme :
: <math>\ddot{s} + 2 \lambda \dot{s} + \omega_0^2 s + c = 0</math>
où λ est un terme d'amortissement et ω<sub>0</sub> la pulsation périodique ou pseudo-périodique. On a donc :
: <math>\dot{\mathrm{Y}} = \begin{pmatrix} x_2 \\ -2 \lambda \dot{x}_1 - \omega_0^2 x_1 - c \end{pmatrix}.</math>
On définit également le facteur de qualité Q :
: <math>\mathrm{Q} = \frac{\omega_0}{2\lambda}.</math>
{{loupe|w:fr:Oscillateur harmonique|w:fr:Équation différentielle linéaire d'ordre deux|v:fr:Équation différentielle/Équation différentielle linéaire du deuxième ordre à coefficients constants}}
Prenons par exemple le cas d'une masse suspendue à un ressort. Le principe fondamental de la dynamique (PFD) s'écrit :
: <math>m\ddot{x} + kx + p = 0</math>
où :
* ''m'' est la masse en kilogrammes (kg) ;
* ''k'' est la raideur du ressort en newton par mètre (N/m) ;
* ''p'' est le poids de l'objet, ''p'' = ''mg'', l'accélération de la gravité valant ''g'' ≃ {{unité|9.81|m/s<sup>2</sup>}} ;
* <math>x</math> est l'altitude de l'objet (le ressort est détendu lorsque ''x'' = 0), <math>\dot{x}</math> est la vitesse (en mètre par seconde, m/s) et <math>\ddot{x}</math> l'accélération (en mètre par seconde au carré, m/s<sup>2</sup>) ;
soit
: <math>\ddot{x} + \frac{k}{m}x + \frac{p}{m} = 0 \Longleftrightarrow \ddot{x} + \omega_0^2 x + g = 0</math>
en posant la pulsation, exprimée en s<sup>–1</sup> :
: <math>\omega_0 = \sqrt{\frac{k}{m}}.</math>
Nous avons donc ici :
* <math>x_1 = x</math>, la position ;
* <math>x_2 = \dot{x}</math>, la vitesse ;
* <math>\dot{\mathrm{Y}} = \begin{pmatrix} x_2 \\ -\omega_0^2 \dot{x}_1 - g \end{pmatrix}.</math>
Nous prenons arbitrairement une pulsation ω = {{unité|1|s<sup>–1</sup>}}, et comme conditions initiales ''x'' = 0 et ''x'' ’ = 0 à ''t'' = 0.
<source lang="scilab">
k = 1; // raideur du ressort en N/m
m = 1; // masse du corps en kg
lmbd = 0;
omega02 = sqrt(k/m);
g = 9.81;
function [Res] = f(t, Y)
Res = [Y(2) ; -2*lmbd*Y(2) - omega02*Y(1) - g]
endfunction
Y0 = [0 ; 0]
t0 = 0
n = 10
t = linspace(0, 2*%pi, n)
pos = []
vit = []
solu = ode(Y0, t0, t, f)
plot(t, solu(1, :), "b")
plot(t, solu(2, :), "r")
legend("position", "vitesse")
</source>
On retrouve bien des courbes sinusoïdales.
Considérons maintenant un circuit électrique RLC en série soumis à un échelon de tension E à ''t'' = 0. La loi des mailles donne
: <math>\mathrm{E} = \mathrm{L}\frac{\mathrm{d}i}{\mathrm{d}t} + \mathrm{R}i + \frac{q}{\mathrm{C}}</math>
soit
: <math>\ddot{q} + \frac{\mathrm{R}}{\mathrm{L}}\dot{q} + \frac{q}{\mathrm{L}\mathrm{C}} = \frac{\mathrm{E}}{\mathrm{L}}</math>
donc
: <math>\lambda = \frac{\mathrm{R}}{2\mathrm{L}}\ ;\ \omega_0^2 = \frac{1}{\mathrm{L}\mathrm{C}}.</math>
<source lang="python">
// [x1' = x2,
// x2' = -(R/L)x2 - 1/(LC)*x1 + E/L]
// donc en posant y = [x1, x2], on a
// y' = [x2, -(R/L)x2 - 1/(LC)*x1 + E/(LC)]
// pour optimiser les calculs, on pose
// R/L = R*C/(L*C) ;
L = 0.5; // inductance en H
R = 2; // résistance en Ω
C = 1e-3; // capacité en F
E = 1; // échelon de tension en V
foo = 1/(L*C);
lmbd = 0.5*R*C*foo;
omega2 = foo;
function [Res] = equation(t, y)
Res = [y(2) ; -2*lmbd*y(2) - omega2*y(1) + E/L];
endfunction
t0 = 0;
T = 2*%pi/sqrt(omega2);
tmax = 6/lmbd;
t = 0:(T/7):tmax;
n = size(t, "*");
y0 = [0 ; 0];
resultat = ode(y0, t0, t, equation);
a1 = newaxes();
a1.tight_limits = "on";
plot(t, resultat(1, :), "b");
xtitle("Circui RLC", "t (s)", "q (C)")
legend("charge", 2);
a2 = newaxes();
a2.filled = "off";
a2.axes_visible(1) = "off";
a2.y_location = "right";
a2.tight_limits = "on";
plot(t, resultat(2, :), "r");
xtitle("", "", "i (A)")
legend("intensité", 1);
</source>
On retrouve bien une courbe de la forme <math>q(t) = \mathrm{A}e^{-\lambda t} \cos \left ( \omega_0 \sqrt{1 - \frac{\lambda^2}{\omega_0^2}} + \varphi \right ).</math>
; Pour aller plus loin
| |||