%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function TraceInterpolee(f,a,b) // // calcule l'interpolée de f // aux 3 points a, b, et (a+b)/2 // et trace les deux courbes sur [a,b] h=b-a;c=(a+b)/2;b d=2*(f(c)-f(a))/h;e=2*(f(a)+f(b)-2*f(c))/h/h; x=linspace(a,b)'; // 100 points de l'intervalle [a,b] z=f(a)+d*(x-a)+e*(x-a).*(x-c) // interpolee de f, calculée sur [a,b] y=f(x); plot2d([x,x],[y,z]) // trace y et z en fonction de x %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [y,z]=projection(A,x) // z=A*x y=x-z; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function c=strassen(a,b) // // attention: les matrices sont carrées // de taille n x n avec n=2^k n=size(a,1); if n==1 then c=a*b; else m=n/2; // on decoupe la matrice a en 4 blocs a11=a(1:m,1:m);a12=a(1:m,m+1:n);a21=a(m+1:n,1:m);a22=a(m+1:n,m+1:n); // idem pour b b11=b(1:m,1:m);b12=b(1:m,m+1:n);b21=b(m+1:n,1:m);b22=b(m+1:n,m+1:n); // regle de calcul de strassen m1=strassen(a12-a22,b21+b22); m2=strassen(a11+a22,b11+b22); m3=strassen(a11-a21,b11+b12); m4=strassen(a11+a12,b22 ); m5=strassen(a11 ,b12-b22); m6=strassen(a22 ,b21-b11); m7=strassen(a21+a22,b11 ); // on definit c=a*b par blocs c=[m1+m2-m4+m6,m4+m5;m6+m7,m2-m3+m5-m7]; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [x,y]=CercleUnite(A,n) // // ajouter un test vérifiant que la // matrice est s.d.p. x=[];y=[];h=2*%pi/n; for alpha=0:h:2*%pi p=tan(alpha); ap=[1 p]*A*[1;p]; if cos(alpha)>0 then x=[x 1/sqrt(ap)]; y=[y p/sqrt(ap)]; else x=[x -1/sqrt(ap)]; y=[y -p/sqrt(ap)]; end; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function Tn=tcheby(n,x) // for i=1:length(x) t=x(i); if abs(t)>1. then y=sqrt(t*t-1); y=(t+y)^n+(t-y)^n;y=y/2; else theta=acos(t); y=cos(n*theta); end; Tn(i)=y; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function B = GS(A) // pres=1.e-12; [m,n]=size(A); B=zeros(m,n); for i=1:n s=zeros(m,1) for k=1:i-1 s=s+(A(:,i)'*B(:,k))*B(:,k); end; s=A(:,i)-s; if norm(s) > pres then B(:,i) =s/norm(s); end; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function B = GS1(A) // pres=1.e-12; colnn=0; // indiquera la colonne non nulle courante [m,n]=size(A); B=zeros(m,n); for i=1:n s=zeros(m,1) for k=1:i-1 // on peut ecrire k=1:colnn s=s+(A(:,i)'*B(:,k))*B(:,k); end; s=A(:,i)-s; if norm(s) > pres then colnn=colnn+1; B(:,colnn) =s/norm(s); end; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function B = GSM(A) // pres=1.e-12; [m,n]=size(A); B=zeros(m,n); for i=1:n s=A(:,i); if norm(s) > pres then B(:,i) =s/norm(s); for k=i+1:n A(:,k)=A(:,k)-(A(:,k)'*B(:,i))*B(:,i); end; else error('Vecteurs lies') end; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function y=DansLimage(A,b) // kerAt=kernel(A'); if (size(kerAt,1)==size(b,1)) then if norm(kerAt'*b) >1.e-6 then y='non'; else y='oui'; end; else error('Probleme de dimension') end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function cercles(x,y,r) // // trace des cercles // centre x_i,y_i, rayons r_i for i=1;length(x) xarc(x(i)-r(i),y(i)+r(i),2*r(i),2*r(i),0,64*360) end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function plotGersch(a) // // trace les cercles de Gerschgorin-Hadamard [m,n]=size(a); if m~=n then, error('la matrice n''est pas carree'), end; d=diag(a); rayons1=sum(abs(a),'c'); rayons1=rayons1-abs(d); // définir un grand rectangle contenant tous les cercles // détermine les coins ``inférieur-gauche'' for i=1:n; coinx1(i)=real(a(i,i))-rayons1(i); coiny1(i)=imag(a(i,i))-rayons1(i); end; mx=min(coinx1);my=min(coiny1); mx=round(mx-1);my=round(my-1); // détermine les coins ``supérieur-droit'' for i=1:n; coinx1(i)=real(a(i,i))+rayons1(i); coiny1(i)=imag(a(i,i))+rayons1(i); end; Mx=max(coinx1);My=max(coiny1); Mx=round(Mx+1);My=round(My+1); // on indique les valeurs propres par un + plot2d(real(spec(a)),imag(spec(a)),-1,"031"," ",[mx,my,Mx,My]) xgrid() square(mx,my,Mx,My) // détermine les coins superieur-gauche // des carrés qui delimitent les cercles pour a for i=1:n; coinx1(i)=real(a(i,i))-rayons1(i); coiny1(i)=imag(a(i,i))+rayons1(i); end; // on trace les cercles for i=1:n; diam=2*rayons1(i); xarc(coinx1(i),coiny1(i),diam,diam,0,64*360); end; xgrid() square(mx,my,Mx,My)