Calcul numérique

Algèbre linéaire

Tout traitement mathématique est vu par Scilab comme une manipulation de matrices. Pour l'algèbre linéaire, les fonctionnalités sont donc particulièrement développées. Dans une utilisation un peu avancée, on aura intérêt à tirer parti de la structure de matrice creuse (sparse) que nous n'aborderons pas.

Il est important de garder à l'esprit que tous les calculs numériques dépendent de la représentation des nombres réels en machine, et que donc les résultats ``exacts'' n'existent pas. La détection de singularités dépend d'une précision prédéfinie et n'est donc pas infaillible.

rank([1,1;1,1])
rank([1.00000000000001,1;1,1])
rank([1.000000000000001,1;1,1])
inv([1,1;1,1])
A=[1,2,3;4,5,6;7,8,9]
det(A)
det(A')
rank(A)
inv(A)     // le calcul est effectue
inv(A')    // message d'erreur
A*inv(A)   // ce n'est pas l'identite

Si A et B sont deux matrices, alors A\B retourne une matrice X telle que A*X=B, et B/A une matrice X telle que X*A=B pourvu que les dimensions de A et B soient compatibles. Si A n'est pas une matrice carrée inversible, ces commandes retourneront un résultat, qu'il y ait une infinité de solutions ou qu'il y en ait aucune.

A=[1,2,3;4,5,6;7,8,9]
b=[1;1;1]
x=A\b       // le systeme a une infinite de solutions
A*x
b=[1,1,1]
x=b/A       // le systeme a une infinite de solutions  
x*A        
b=[1;1;2]
x=A\b       // le systeme n'a pas de solution
A*x
b=[1,1,2]
x=b/A       // le systeme n'a pas de solution 
x*A

Si A est carrée et régulière x=A\b est unique et équivalent mathématiquement à x=inv(A)*b mais plus rapide.

b=[1:500]';
A=diag(b);
timer(); A\b ; timer()
timer(); inv(A)*b ; timer()

Pour résoudre un système linéaire quand on craint une singularité, il vaut mieux utiliser linsolve, qui détecte les impossibilités et peut retourner une base du noyau de la matrice en cas de solutions multiples. Attention, linsolve(A,b) résout A*x=-b.

A=[1,2,3;4,5,6;7,8,9]
b=[1;1;1]
[x,k]=linsolve(A,b)
A*x
A*k
b=[1;1;2]
[x,k]=linsolve(A,b)  // erreur

La commande [D,U]=bdiag(A) retourne (en théorie) une matrice diagonale D et une matrice inversible U telles que U*D*inv(U)=A.

A=diag([1,2,3]); P=rand(3,3); A=P*A*inv(P)
spec(A)
[D,U]=bdiag(A)       // matrice diagonalisable
inv(U)*A*U
U*D*inv(U)-A
A=[1,0,0;0,1,1;0,0,1]; P=rand(3,3); A=P*A*inv(P)
spec(A)
[D,U]=bdiag(A)       // matrice non diagonalisable
inv(U)*A*U
U*D*inv(U)-A
[X,diagevals]=spec(A)

La décomposition LU d'une matrice, c-à-d la décomposition en produit d'une matrice triangulaire inférieure L et d'une matrice triangulaire supérieure U, à des permutations près, s'obtient grâce à la commande lu.

A= [ 1, 2, 3; 4, 5, 6 ; 7, 8, 9  ]
[L,U,E]=lu(A)
E*A-L*U
det(U)
det(A)

Opérations matricielles
A'	transposée de A
rank	rang
lu	décomposition LU
inv	inverse
expm	exponentielle matricielle
det	déterminant
trace	trace
poly(A,"x")	polynôme caractéristique de A
spec	valeurs propres de A
bdiag	diagonalisation
svd	décomposition en valeurs singulières
A\b	solution de A*x = b
b/A	solution de x*A = b
linsolve(A,b)	solution de A*x = -b

Intégration

Des fonctions d'intégration numérique sont disponibles pour les fonctions réelles et complexes, à une, deux et trois variables. Les fonctions integrate, intg, int2d, int3d, intc et intl prennent en entrée une fonction externe, ou définie par une chaîne de caractères. Les fonctions inttrap et intsplin prennent en entrée des vecteurs d'abscisses et d'ordonnées.

Calculs d'intégrales
integrate	fonction définie par une chaîne de caractères
intg	fonction externe
inttrap	méthode des trapèzes
intsplin	approximation par splines
int2d	fonction de deux variables
int3d	fonction de trois variables
intc	fonction complexe le long d'un segment
intl	fonction complexe le long d'un arc de cercle

Dans l'exemple ci-dessous, on calcule avec les différentes fonctions d'intégration disponibles, la valeur de :

$$\int_{-10}^0 e^x dx\;=\;1-e^{-10}\;\simeq\;0.9999546\;.$$

x = [-10:0.1:0];
y=exp(x);
1-1/%e^10
inttrap(x,y)
intsplin(x,y)
integrate("exp(x)","x",-10,0)
deff("y=f(x)","y=exp(x)")
intg(-10,0,f)

Les algorithmes de calcul numérique des transformées classiques sont disponibles. Le tableau ci-dessous en donne quelques-unes, voir apropos transform pour les autres.

Transformées
dft	transformée de Fourier discrète
fft	transformée de Fourier rapide
convol	produit de convolution

Résolutions et optimisation

La fonction fsolve résout numériquement un système d'équations, mis sous la forme $f(x)=0$ . Comme toutes les résolutions numériques, celle-ci part d'une valeur initiale $x_0$ , et itère une suite censée converger vers une solution. Le résultat dépend évidemment de la valeur initiale.

deff("y=f(x)","y=sin(%pi*x)")
fsolve(0.2,f)
fsolve(0.4,f)
fsolve(0.45,f)
fsolve(0.5,f)
fsolve([0.45:0.01:0.5],f)
help fsolve
[x,v,info]=fsolve(0.5,f)
[x,v,info]=fsolve([0.45:0.01:0.5],f)

Résolutions
fsolve	systèmes d'équations
roots	racines d'un polynôme
factors	facteurs irréductibles réels d'un polynôme
linsolve	systèmes linéaires

La fonction optim recherche un minimum (local) d'une fonction dont on connaît le gradient. La définition de la fonction en entrée est un peu particulière. Le choix de trois algorithmes différents est offert en option.

help optim
deff("[y,yprim,ind]=f(x,ind)","y=sin(%pi*x),yprim=%pi*cos(%pi*x)")
[x,v]=optim(f,0.2)
[v,xopt]=optim(f,0.50000000000000001)
[v,xopt]=optim(f,0.5000000000000001)

Optimisation
optim	optimisation
linpro	programmation linéaire
quapro	programmation quadratique

Equations différentielles

La fonction ode est en fait un environnement qui donne accès à la plupart des méthodes numériques classiques pour la résolution des équations différentielles ordinaires (voir help ode).

A titre d'exemple nous commençons par le problème de Cauchy en dimension 1 suivant :

$$\left\{ \begin{array}{lcl} y'(t) &=& y(t) \cos(t)\;,\ y(0) &=& 1\;. \end{array} \right.$$

La solution explicite est $y(t)=\exp(\sin(t))$ . La résolution numérique par ode est très stable.

deff("yprim=f(t,y)","yprim=y*cos(t)")
t0=0; y0=1; t=[0:0.01:10];
sol=ode(y0,t0,t,f);
max(abs(sol-exp(sin(t))))     // l'erreur est tres faible
t0=0; y0=1; t=[0:0.1:100];
sol=ode(y0,t0,t,f);
max(abs(sol-exp(sin(t))))     // l'erreur reste tres faible

La situation n'est pas toujours aussi favorable. Considérons par exemple $y(t)=\exp(t)$ , solution du problème de Cauchy suivant.

$$\left\{ \begin{array}{lcl} y'(t) &=& y(t)\;,\ y(0) &=& 1\;. \end{array} \right.$$

deff("yprim=f(t,y)","yprim=y")
t0=0; y0=1; t=[0:0.01:10];
sol=ode(y0,t0,t,f);
max(abs(sol-exp(t)))          // l'erreur est raisonnable
t0=0; y0=1; t=[0:0.1:100];
sol=ode(y0,t0,t,f);
max(abs(sol-exp(t)))          // l'erreur est enorme

Voici en dimension 2 la résolution d'un système différentiel de Lotka-Volterra, avec superposition de la trajectoire calculée et du champ de vecteurs défini par le système (figure 8).

//
// definition du systeme 
//
deff("yprim=f(t,y)",..
    ["yprim1=y(1)-y(1)*y(2)";..
     "yprim2=-2*y(2)+2*y(1)*y(2)";..
     "yprim=[yprim1;yprim2]"])
//
// champ de vecteurs
//
xmin=0; xmax=3; ymin=0; ymax=3;
fx=xmin:0.3:xmax; fy=ymin:0.3:ymax;
clf()
fchamp(f,1,fx,fy)
//
//resolution du systeme
//
t0=0; tmax=5; pas=0.1
t=t0:pas:tmax;
y0=[2;2];
sol=ode(y0,t0,t,f);
//
// trace de la trajectoire
//
plot2d(sol(1,:),sol(2,:),5,"000")

Figure 8: Résolution d'un système de Lotka-Volterra par ode, représentation de plusieurs solutions et du champ de vecteurs.