tel/TPmatlab/equation_nonlineaire/TP4m/sqrtm_iteration_ordre3.m

   1 function [er,X]=sqrtm_iteration_ordre3(A,X0,epsilon)\r
   2 \r
   3 % \r
   4 %  sqrtm_iteration_ordre3 : calcul d'une racine carrée de matrice par la méthode de Newton.\r
   5 %\r
   6 %\r
   7 % *********************************************************\r
   8 %\r
   9 %  [er,X]=sqrtm_iteration_ordre3(A,X0,epsilon)\r
  10 %      calcule les différents itérés de\r
  11 %      la méthode d'ordre 3 pour la recherche d'une racine de matrice.\r
  12 %      A chaque itération, elle affiche et stocke le résultat X_k ainsi que l'erreur\r
  13 %         er(k)=|||X_k^2-A|||.\r
  14 %       \r
  15 %       variables d'entrées :\r
  16 %    * A :matrice dont on veut une racine.\r
  17 %    * X0 : matrice initiale\r
  18 %    * epsilon : paramètre de précision\r
  19 %       variables de sortie :\r
  20 %    * er: suite des erreurs, définies par er(k)=norm(X_k^2-A).\r
  21 %    * X : suite des itérés (matriciels) calculés.\r
  22 % \r
  23 % ATTENTION : nombre maximum d'itération défini par nmax\r
  24 % \r
  25 % ************ Fonctions auxiliaires utilisées ************\r
  26 %\r
  27 %       aucune\r
  28 %\r
  29 % *********************************************************\r
  30 %\r
  31 \r
  32 % nombre maximum d'itération défini par nmax\r
  33 nmax=100;\r
  34 \r
  35 % Contrôles d'entrée\r
  36 [n,p]=size(A);\r
  37 if (n~=p) \r
  38    error('attention, A n''est pas carrée : arrêt du programme');\r
  39 end\r
  40 sqrtmA=sqrtm(A);\r
  41 if (norm(sqrtmA*X0-X0*sqrtmA,inf)~=0)\r
  42    disp('attention, la matrice initiale et la racine de A ne commutent pas !!');\r
  43 end\r
  44 \r
  45 \r
  46 % Corps d'algorithme\r
  47 er=[];\r
  48 X=X0;\r
  49 x=X0;\r
  50 k=0;\r
  51 A2=A^2;\r
  52 test=(norm(x^2-A,inf)~=0);\r
  53 while(test)\r
  54    z=inv(x);\r
  55    y=(3/8)*x+(3/4)*A*z-(1/8)*A2*(z^3);\r
  56    disp(['itération numéro ',int2str(k)]);\r
  57    y2=y^2;\r
  58    erloc=norm(y2-A);\r
  59    er=[er erloc];\r
  60    X=[X y];\r
  61    test=~((norm(y2-A,inf)==0)|(k>=nmax)|(norm(x-y,inf)<=epsilon));\r
  62    disp('itéré : ');\r
  63    disp(y);\r
  64    disp(['erreur : ',num2str(erloc)]);\r
  65    disp('appuyez sur une touche pour continuer');\r
  66    pause;\r
  67    k=k+1;\r
  68    x=y;\r
  69 end\r
  70 if (n>=nmax)\r
  71    disp('attention, n>nmax');\r
  72 end\r