7. 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 codice_1 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 codice_2 ou bien par codice_3), alors on peut utiliser la fonction codice_4 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 codice_5 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
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 codice_6. 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 codice_7 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 codice_8 :
x = linspace(0, 4, 50);
y = f(x);
inttrap(x, y)
Pour une liste de valeurs.
La fonction codice_9 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 codice_8 permet de calculer une intégrale par la méthode des trapèzes.
I = inttrap(X, Y) // I = 22
La fonction codice_11 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) ("hi")1 ≤ "i" ≤ "m" et ("xi")1 ≤ "i" ≤ "n", avec "m" < "n". Le produit de convolution discret est défini par :
Il est obtenu par la commande codice_12 ou bien la commande codice_13. 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 codice_14 et codice_15 n'utilisent pas le même algorithme. La fonction codice_15 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 codice_17. Soit une fonction "y"("t") ; si l'équation différentielle est de premier ordre
alors ƒ ayant été définie (fonction externe), la syntaxe pour déterminer "y"("t") est
où
La fonction interne codice_17 admet des arguments permettant de résoudre des situations spécifiques : 
-->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
Équation différentielle ordinaire d'ordre plus élevé.
Une ÉDO d'ordre supérieur peut se réduire à une ÉDO d'ordre 1.
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 :
En posant
l'équation s'écrit
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 :
Prenons par exemple le cas d'une masse suspendue à un ressort. Le principe fondamental de la dynamique (PFD) s'écrit :
où :
soit
en posant la pulsation, exprimée en s–1 :
Nous avons donc ici :
Nous prenons arbitrairement une pulsation ω = , 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 formula_21
Choix du solveur.
Le solveur par défaut est LSODA, développé par le laboratoire Livermore. 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 : codice_23, codice_24, codice_25 ou codice_26
De manière générale, on distingue deux types de solveurs :
Options pour le solveur.
La commande codice_31 ouvre une fenêtre permettant de modifier de manière interactive les paramètres du solveur. On peut également modifier la variable codice_32 :
%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 :
Voir aussi.
"Régression" < ↑ > "Résolution d'équations"