mode(7)
//Le script aide_TP4 étudie les erreurs et le temps mis pour
//résoudre un système linéaire par différentes méthodes.

disp('Norme matricielle :')
A=[1 2; 3 4]
norm(A,1), norm(A,2), norm(A,'inf'), norm(A,'fro')
sqrt(max(spec(A'*A)))
norm(A,2)==sqrt(max(spec(A'*A)))  //attentions aux arrondis !

disp('Systèmes linéaires et erreurs d''arrondis :')
A=[1 2 3 ; 4 5 6 ; 7 8 10]
condA=cond(A) //erreurs importantes quand cond(A) est grand
X0=[1 1 1]'; B=A*X0 //solution X0 de l'équation AX=B
X1=inv(A)*B;
norm_r1=norm(A*X1-B) //norme du "résidu" par inversion de  A
err1=norm(X1-X0)     //erreur par inversion de A
X2=A\B;
norm_r2=norm(A*X2-B)    //norme du "résidu" par "\"
err2=norm(X2-X0)        //erreur par "\"
Xgauss=lusolve(sparse(A),B); //on impose ici la méthode de Gauss
//variante : "rref" échelonne une matrice A|B, de gauche à droite
//A_aug=[A B]; A_aug_ech=rref(A_aug);
//Xgaussbis=A_aug_ech(:,size(A,'c')+1);
norm_r_gauss=norm(A*Xgauss-B)    //norme du "résidu" par Gauss
errgauss=norm(Xgauss-X0)        //erreur par Gauss

disp('Systèmes linéaires et moindres carrés :')
Abis=A; Abis(3,3)=9        //on modifie A en Abis
condAbis=cond(Abis), rgAbis=rank(Abis)    //Abis non inversible
Base_KerAbis=kernel(Abis) //base du noyau de Abis
X3=pinv(Abis)*B;   //résolution mod Ker(Abis) par moindres carrés
norm_r3=norm(Abis*X3-B) //norme du "résidu"
err3=norm(X3-X0) //erreur minimale
X4=Abis\B;         //cf. "help \", puis "help backslash"
norm_r4=norm(Abis*X4-B) //norme du "résidu" par Gauss "\"
err4=norm(X4-X0)         //erreur par Gauss "\"

disp('Temps mis pour résoudre un système linéaire en interne :')
n=200; An=rand(n,n); Bn=rand(n,1);
timer(); An\Bn; t=timer();         //résolution de AnX=Bn
timer(); inv(An); tbis=timer();     //inversion de An
timer(); An*An; tter=timer(); clear An Bn; //calcul de An^2
v=[rand(1,10000)+1;1 ./(rand(1,10000)+1)];
timer(); for i=1:100, s=sum(v); end; ts=timer()/2000000; //ts : +
timer(); for i=1:100, p=prod(v); end; tp=timer()/2000000; //tp : *
printf(' \nQuand n=%3d, Scilab met environ \n..
%.2f secondes pour résoudre un système lineaire nxn,\n..
%.2f secondes pour inverser une matrice nxn,\n..
%.2f secondes pour calculer le carré d''une matrice nxn,\n..
et environ %.2f secondes pour effectuer n^3 opérations itérées.\n\n'..
,n,t,tbis,tter,((ts+tp)/2)*n^3)
clear

disp('Stockage des matrices creuses :')
C=floor(10*sprand(10,10,0.20)); //matrice creuse à 80%
P=full(C)  //full(C) : matrice pleine associée à C (permet "norm(P)")
[noms,memoire]=who('local'); //données locales en mémoire
memoire(find(noms=='C')) //nombre de mémoires occupées par C
memoire(find(noms=='P')) //nombre de mémoires occupées par P
exec spy.sci;
spy(C)    //emplacements des termes non nuls de C
[i,j] = find(C); indices_termes_non_nuls=[i',j']

clear