Découvrir Scilab/Calcul différentiel et intégral

Dérivation

La dérivation peut s'aborder de différentes manières selon les données disponibles.

Pour une fonction définie

S'il s'agit d'un polynôme ou d'une fraction rationnelle, la commande derivat() effectue la dérivée formelle, par exemple

-->q = (1 + %s^2)/(1+%s)
 q  =
 
         2  
    1 + s   
    -----   
    1 + s   
 
-->derivat(q)
 ans  =
 
              2  
  - 1 + 2s + s   
    ----------   
              2  
    1 + 2s + s

Si l'on dispose d'une fonction extérieure (définie par function() … endfunction ou bien par deff()), alors on peut utiliser la fonction numderivative()^[1] pour obtenir l'approximation de la matrice jacobienne et hessienne :

-->deff("y = f(x)", "y = x.^2")

-->numderivative(f, 2) // f'(x) = 2*x ; J = f'(2) = 4
 ans  =
 
    4.

-->[J, H] = numderivative(f, 2) // f''(x) = 2 ; H = f''(2) = 2
 H  =
 
    2.  
 J  =
 
    4.  

-->x=(1:3)';

-->numderivative(f, x)
 ans  =
 
    2.    0.    0.   
    0.    4.    0.  
    0.    0.    6.   

-->[J, H] = numderivative(f, x)
 H  =
 
    2.    0.    0.    0.    0.    0.    0.    0.    0.  
    0.    0.    0.    0.    2.    0.    0.    0.    0.  
    0.    0.    0.    0.    0.    0.    0.    0.    2.  
 J  =
 
    2.    0.    0.  
    0.    4.    0.  
    0.    0.    6.

Pour une liste de valeurs

Si l'on dispose d'une liste de valeurs, typiquement un vecteur x et un vecteur y de même dimension, il faut alors distinguer plusieurs cas.

Valeurs peu bruitées

Si l'on suppose que le bruit est faible, donc que les valeurs sont « exactes », et que de plus la dérivée seconde est faible, alors on peut se contenter de déterminer les pentes entre les points, c'est-à-dire de calculer le taux de variation, faire une interpolation linéaire. Par exemple :

x = [0 2 5 20 40 60 80];   
y = [0.0956820 0.0480457 0.0277857 0.0036214 0.0002543 0.0002543 0.0001265];

yprime = diff(y)./diff(x);

subplot(2, 1, 1)
plot(x, y)

subplot(2, 1, 2)
plot(x(1:$-1), yprime)

Si la dérivée seconde n'est pas négligeable, on peut alors utiliser une interpolation par un polynôme de degré trois (spline cubique) :

yprime1 = splin(x, y);

subplot(2, 1, 2)
plot(x, yprime1, "r")

Valeurs bruitées : régression et lissage

Les méthodes ci-dessus sont très sensibles au bruit. Si les valeurs sont bruitées, alors il faut lisser la courbe d'une manière ou d'une autre.

Si l'on dispose d'un modèle analytique pour les données, on peut déterminer les paramètre de ce modèle par régression, voir la section ci-dessus. Sinon, l'utilisation d'une spline par la commande lsq_splin() présentée précédemment permet également d'obtenir la dérivée : il suffit d'évaluer la spline avec la syntaxe suivante :

[yi, yprime2] = interp(x, xs, ys, d); // valeurs de la spline aux abscisses x
// et valeurs de la dérivée

Détermination des dérivée et dérivée seconde par l'algorithme de Savitzky-Golay.

L'algorithme de Savitzky-Golay présenté précédemment permet également d'obtenir la dérivée : comme on détermine localement un polynôme, on a les dérivées successives de ce polynôme.

Intégration

Intégrale numérique

Pour une fonction définie

Scilab peut calculer l'intégrale numérique d'une fonction entre deux bornes :

integrate("expression", "x", x0, x1)

calcule l’intégrale de la fonction décrite par expression ; expression est une chaîne de caractères interprétable par Scilab, comme par exemple "sin(x)". Le paramètre x est la variable d’intégration, et les bornes sont x0 et x1.

Par exemple :

// fonction carré
function [y] = f(x)
    y = x.^2
endfunction

// intégrale entre 0 et 4
I = integrate("f(x)", "x", 0, 4) // I = 21.3333

// vérification
I1 = 4^3/3

On peut aussi utiliser la fonction intg(a, b, f) où a et b sont les bornes d'intégration. Avec l'exemple précédent :

-->Y = intg(0, 4, f)
 Y  =
 
    21.333333

Il peut également effectuer une intégration numérique à partir d'une suite de points (x(i), y(i)), par la méthode des trapèzes, avec la commande inttrap() :

x = linspace(0, 4, 50);
y = f(x);
inttrap(x, y)

Pour une liste de valeurs

La fonction sum() permet de cumuler les valeurs. Si les valeurs de y correspondent à des valeurs de x équidistantes d'une largeur h, alors une première approche de l'intégrale (méthode du point milieu) est :

X = 0:4;
Y = X.*X;

h = 1;

I = sum(Y)*h // I = 30

La fonction inttrap() permet de calculer une intégrale par la méthode des trapèzes.

I = inttrap(X, Y) // I = 22

La fonction intsplin() permet de calculer une intégrale en interpolant les points par une cerce (spline).

I = intsplin(X, Y) // I = 21.333333

Produit de convolution

Considérons deux suite de valeurs (vecteurs ligne) (h_i)_{1 ≤ i ≤ m} et (x_i)_{1 ≤ i ≤ n}, avec m < n. Le produit de convolution discret est défini par :

y_{j}=\sum _{i=1}^{m}h_{j}\times x_{k+1-j}

Il est obtenu par la commande conv(h, x) ou bien la commande convol(h, x). Par exemple, pour convoluer une fonction sinus avec une fonction porte de largeur 5 points :

t = linspace(0, %pi, 50);
x = sin(t);
h = 1/5*ones(1, 5);
y = conv(h, x);
plot(x)
plot(y, "r")

Les commandes conv et convol n'utilisent pas le même algorithme. La fonction convol utilise la transformée de Fourier rapide, et est mieux adaptée lorsqu'il y a beaucoup de valeurs.

É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 ode. Soit une fonction y(t) ; si l'équation différentielle est de premier ordre

dy/dt = ƒ(t,y),

alors ƒ ayant été définie (fonction externe), la syntaxe pour déterminer y(t) est

y = ode(y0,t0,t,f)

où

y0 et t0 sont les valeurs initiales du système, y(t0) = y0,
t est un vecteur de valeurs pour lesquelles on calcule les solutions, et
y est le vecteur de solutions (plot(t,y) permet de tracer le graphique des solutions).

La fonction interne ode admet des arguments permettant de résoudre des situations spécifiques :

"roots" pour rajouter une équation g(t,y) = 0,
"discrete" pour calculer de manière récursive y(k+1) = ƒ(k,y(k)) à partir d'un état initial y(k0).

Exemple

On se place sur un segment [0,5;1] ; on veut déterminer les valeurs numériques de la fonction y vérifiant

{\frac {\mathrm {d} y}{\mathrm {d} t}}=t

et

y(0)=0

on a donc en fait ƒ(t,y) = t , t0 = y0 = 0

-->t=[0.5:0.1:1]; // points de [0,5;1] espacés de 0,1

-->deff("[u] = f(t,y)", "u = t") // définition de la fonction f

-->ode(0,0,t,f)
 ans  =
    0.125    0.18    0.245    0.32    0.405    0.5

qui sont bien les valeurs de 1/2×t²

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¹ 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.

Pour plus de détails voir : 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₁ = y

x₂ = y ’

l'équation s'écrit

x₂’ + ax₁’ + bx₁ + c = 0.

En posant le vecteur

\mathrm {Y} ={\begin{pmatrix}x_{1}\\x_{2}\end{pmatrix}}

on a

\mathrm {Y} '={\begin{pmatrix}x_{1}'\\x_{2}'\end{pmatrix}}={\begin{pmatrix}x_{2}\\-ax_{1}'-bx_{1}-c\end{pmatrix}}.

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 :

{\ddot {s}}+2\lambda {\dot {s}}+\omega _{0}^{2}s+c=0

où λ est un terme d'amortissement et ω₀ la pulsation périodique ou pseudo-périodique. On a donc :

{\dot {\mathrm {Y} }}={\begin{pmatrix}x_{2}\\-2\lambda {\dot {x}}_{1}-\omega _{0}^{2}x_{1}-c\end{pmatrix}}.

On définit également le facteur de qualité Q :

\mathrm {Q} ={\frac {\omega _{0}}{2\lambda }}.

Pour plus de détails voir : w:fr:Oscillateur harmonique, w:fr:Équation différentielle linéaire d'ordre deux et 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 :

m{\ddot {x}}+kx+p=0

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 ≃ 9,81 m/s² ;
$x$ est l'altitude de l'objet (le ressort est détendu lorsque x = 0), ${\dot {x}}$ est la vitesse (en mètre par seconde, m/s) et ${\ddot {x}}$ l'accélération (en mètre par seconde au carré, m/s²) ;

soit

{\ddot {x}}+{\frac {k}{m}}x+{\frac {p}{m}}=0\Longleftrightarrow {\ddot {x}}+\omega _{0}^{2}x+g=0

en posant la pulsation, exprimée en s^–1 :

\omega _{0}={\sqrt {\frac {k}{m}}}.

Nous avons donc ici :

$x_{1}=x$ , la position ;
$x_{2}={\dot {x}}$ , la vitesse ;
${\dot {\mathrm {Y} }}={\begin{pmatrix}x_{2}\\-\omega _{0}^{2}{\dot {x}}_{1}-g\end{pmatrix}}.$

Nous prenons arbitrairement une pulsation ω = 1 s^–1, et comme conditions initiales x = 0 et x ’ = 0 à t = 0.

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")

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

\mathrm {E} =\mathrm {L} {\frac {\mathrm {d} i}{\mathrm {d} t}}+\mathrm {R} i+{\frac {q}{\mathrm {C} }}

soit

{\ddot {q}}+{\frac {\mathrm {R} }{\mathrm {L} }}{\dot {q}}+{\frac {q}{\mathrm {L} \mathrm {C} }}={\frac {\mathrm {E} }{\mathrm {L} }}

donc

\lambda ={\frac {\mathrm {R} }{2\mathrm {L} }}\ ;\ \omega _{0}^{2}={\frac {1}{\mathrm {L} \mathrm {C} }}.

// [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);

On retrouve bien une courbe de la forme $q(t)=\mathrm {A} e^{-\lambda t}\cos \left(\omega _{0}{\sqrt {1-{\frac {\lambda ^{2}}{\omega _{0}^{2}}}}}+\varphi \right).$

Pour aller plus loin: Voir G. Sallet, « Ordinary Differential Equations with Scilab – WATS Lectures, Provisional notes (Université de Saint-Louis) », sur Université de Metz, 2004 (consulté le 23 février 2018)

Choix du solveur

Le solveur par défaut est LSODA, développé par le laboratoire Livermore^[2]. Selon la raideur du problème, il utilise la méthode Adams ou bien BDF. On peut utiliser un autre solveur en indiquant le type comme premier paramètre : adams, stiff, rk ou rkf

De manière générale, on distingue deux types de solveurs :

les solveurs explicites, adaptés aux équations peu raides : adams (méthodes d'Adams-Bashforth), rk (méthode de Runge-Kutta d'ordre 4), rkf (méthode de Runge-Kutta-Fehlberg pair d'ordre 4 et 5) ;
les solveurs implicites, adaptés aux équations raides : stiff (backward differentiation formula, BDF).

Pour plus de détails voir : w:fr:Équation différentielle raide, w:fr:Formulation implicite ou explicite d'un problème de dynamique, w:fr:Méthodes de Runge-Kutta et w:fr:Méthodes d'Adams-Bashforth.

Options pour le solveur

La commande odeoptions() ouvre une fenêtre permettant de modifier de manière interactive les paramètres du solveur. On peut également modifier la variable %ODEOPTIONS :

%ODEOPTIONS = [itask, tcrit, h0, hmax, hmin, jactyp, mxstep, maxordn, maxords, ixpr, ml, mu]
// par défaut :
// %ODEOPTIONS = [1.   0.   0.   Inf   0.   2.   500.   12.   5.   0.  -1.  -1]

avec :

h0 : premier pas ;
hmax : pas maximum ;
hmin : pas minimum.

Notes

↑ les fonctions derivative() et numdiff() sont déclarées obsolètes et ont été supprimées à partir de la version 6
↑ (en) « Description and Use of LSODE, the Livermore Solver for Ordinary Differential Equations », sur Lawrence Livermore National Laboratory (consulté le 21 mai 2019).

Voir aussi

Dans Wikipédia

Équation différentielle ordinaire

Régression < ↑ > Résolution d'équations

[1] s fonctions derivative() et numdiff() sont déclarées obsolètes et ont été supprimées à partir de la version 6

[2] (en) « Description and Use of LSODE, the Livermore Solver for Ordinary Differential Equations », sur Lawrence Livermore National Laboratory (consulté le 21 mai 2019).

[1]

[2]