tel/TPmatlab/equation_nonlineaire/TP4e/sqrtm_iteration_newton.m

   1 function [er,X]=sqrtm_iteration_newton(A,X0,epsilon)\r
   2 \r
   3 % \r
   4 %  sqrtm_iteration_newton : 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_newton(A,X0,espilon)\r
  10 %      calcule les différents itérés de\r
  11 %      la méthode de Newton 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 if (norm(A*X0-X0*A,inf)~=0)\r
  41    disp('attention, la matrice initiale et A ne commutent pas !!');\r
  42 end\r
  43 \r
  44    \r
  45 \r
  46 % Corps d'algorithme\r
  47 er=[];\r
  48 x=X0;\r
  49 X=x;\r
  50 k=0;\r
  51 test=(norm(X0^2-A,inf)~=0);\r
  52 while(test)\r
  53    y=(x+A/x)/2;\r
  54    disp(['itération numéro ',int2str(k)]);\r
  55    y2=y^2;\r
  56    erloc=norm(y2-A);\r
  57    er=[er erloc];\r
  58    X=[X y];\r
  59    test=~((norm(y2-A,inf)==0)|(k>=nmax)|(norm(x-y,inf)<=epsilon));\r
  60    disp('itéré : ');\r
  61    disp(y);\r
  62    disp(['erreur : ',num2str(erloc)]);\r
  63    disp('appuyez sur une touche pour continuer');\r
  64    pause;\r
  65    k=k+1;\r
  66    x=y;\r
  67 end\r
  68 if (n>=nmax)\r
  69    disp('attention, n>nmax');\r
  70 end\r