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

Un livre de Wikilivres.
Sauter à la navigation Sauter à la recherche

Table des matièresIndex



7. Calcul différentiel et intégral


Dérivation[modifier | modifier le wikicode]

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

Pour une fonction définie[modifier | modifier le wikicode]

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[modifier | modifier le wikicode]

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[modifier | modifier le wikicode]

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[modifier | modifier le wikicode]

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[modifier | modifier le wikicode]

Intégrale numérique[modifier | modifier le wikicode]

Pour une fonction définie[modifier | modifier le wikicode]

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)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[modifier | modifier le wikicode]

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[modifier | modifier le wikicode]

Considérons deux suite de valeurs (vecteurs ligne) (hi)1 ≤ im et (xi)1 ≤ in, avec m < n. Le produit de convolution discret est défini par :

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[modifier | modifier le wikicode]

Équation différentielle ordinaire de premier ordre[modifier | modifier le wikicode]

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)

  • 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
et
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 D1 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é[modifier | modifier le wikicode]

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

x1 = y
x2 = y ’

l'équation s'écrit

x2’ + ax1’ + bx1 + c = 0.

En posant le vecteur

on a

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 :

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

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


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 :

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/s2 ;
  • est l'altitude de l'objet (le ressort est détendu lorsque x = 0), est la vitesse (en mètre par seconde, m/s) et l'accélération (en mètre par seconde au carré, m/s2) ;

soit

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

Nous avons donc ici :

  • , la position ;
  • , la vitesse ;

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

soit

donc

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


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 2018-02-23

Choix du solveur[modifier | modifier le wikicode]

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[modifier | modifier le wikicode]

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[modifier | modifier le wikicode]

  1. les 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 2019-05-21.

Voir aussi[modifier | modifier le wikicode]

Dans Wikipédia

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