tel/TPmatlab/equation_nonlineaire/TP4h/muller_elem.m

   1 function [rac,nb_iter]=muller_elem(init,exp_f,arret)\r
   2 % Produit une racine d'équation (E) f(x)=0 par méthode de Muller.\r
   3 % Le calcul se termine par un forçage à une écriture symbolique de rac,\r
   4 % qui est affichée en plus de la numérique!\r
   5 % Ceci peut etre aisément modifié!\r
   6 %\r
   7 % variables d'entrée\r
   8 % init est un vecteur de trois valeurs initiales distinctes\r
   9 % pour la méthode de Muller  classique. On peut généraliser!\r
  10 % exp_f est l'expression de la fonction f sous la forme 'f(x)';\r
  11 % arret est un vecteur de trois réels [f_seuil,nul_seuil,n_max]\r
  12 % dont la signification est la suivante:\r
  13 % si |f(x)|<f_seuil on considère f(x) nul;\r
  14 % si |x|<nul_seuil on considère x nul;\r
  15 % n_max est un entier naturel qui désigne le nombre maximal\r
  16 % d'itérations autorisées.\r
  17 %\r
  18 % variables de sortie\r
  19 % rac est une valeur approchée de la racine obtenue à partir d'init\r
  20 % par algorithme de Muller classique; elle peut etre réelle ou complexe.\r
  21 % nb_iter désigne le nombre d'itérations exécutée pour sortir rac.\r
  22 %\r
  23 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
  24 %                       Fonctions connexes appelées\r
  25 %\r
  26 %\r
  27 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
  28 % spécification\r
  29 format long e;\r
  30 \r
  31 % Tests d'entrée\r
  32 \r
  33 if nargin~=3\r
  34     error('Nombre des arguments incorrect, ou écriture inattendue');\r
  35 end\r
  36 % Les tests de controle suivants peuvent etre plus ou moins développés,\r
  37 % selon la dangerosité pressentie de l'utilisateur... ordinaire!\r
  38 if (size(init,1)~=1)|(size(init,2)~=3)|((sum(diff(init)==0)>0))\r
  39     error('Champ init mal écrit. Voir controle');\r
  40 end\r
  41 if (size(arret,1)~=1)|(size(arret,2)~=3)\r
  42     error('Champ arret mal écrit. Voir controle');\r
  43 end\r
  44 % compléments éventuels\r
  45 \r
  46 \r
  47 % Algorithme proprement dit\r
  48 \r
  49 % Etape 1:initialisations\r
  50 X=init;Y=zeros(1,3);\r
  51 \r
  52 %  La suite serait avantageusement remplacée par une écriture vectorielle\r
  53 % mais l'utilisateur calculera faux sans le voir..; Il faudrait écrire\r
  54 % 'x.^2+1' au lieu de 'x^2+1'\r
  55 for k=1:3\r
  56     x=X(k);Y(k)=eval(exp_f);\r
  57 end\r
  58 \r
  59 d=diff(X);dif_div=diff(Y)./d;\r
  60 nb_iter=0;\r
  61 bool_arret=0;\r
  62 \r
  63 % Etape 2\r
  64 while bool_arret==0\r
  65     nb_iter=nb_iter+1;\r
  66     % détermination du polynome d'interpolation de f sur le support x\r
  67     % et résolution de l'équation approchée\r
  68     An=diff(dif_div)/sum(d);\r
  69     Bn=dif_div(2)+An*d(2);\r
  70     x=X(3);Cn=eval(exp_f);\r
  71     if abs(Cn)<arret(1)\r
  72         bool_arret=1;\r
  73     else\r
  74         % détermination du bon epsilon fournissant le dénominateur optimal\r
  75         den=Bn+sqrt(Bn^2-4*An*Cn).*[1 -1];[val,ind]=max(abs(den));\r
  76         den_opt=den(ind);\r
  77         ajout=-2*Cn/den_opt;\r
  78         \r
  79         % préparation de la boucle suivante\r
  80         \r
  81         X(1:2)=X(2:3);X(3)=X(2)+ajout;\r
  82         Y(1:2)=Y(2:3);\r
  83         x=X(3);Y(3)=eval(exp_f);\r
  84         d=diff(X);dif_div=diff(Y)./d;\r
  85         %mise à jour de bool_arret\r
  86         bool_arret=...\r
  87             (nb_iter>arret(3))|(abs(ajout)<arret(2))|(abs(Y(3))<arret(1));\r
  88        \r
  89     end\r
  90 end\r
  91 \r
  92 % toilettage de la racine:\r
  93 % on la force réelle si elle l'est presque (au sens de nul_seuil)\r
  94 % si elle est complexe, on renvoie aussi sa conjuguée dont on sait\r
  95 % qu'elle est aussi solution par théorie des corps si l'expression\r
  96 % de f est à coefficient réels.\r
  97 \r
  98 if abs(real(X(3)))<arret(2)\r
  99     X(3)=imag(X(3))*i;\r
 100 end\r
 101 if abs(imag(X(3)))<arret(2)\r
 102     rac=real(X(3));\r
 103 else\r
 104     rac=[X(3) conj(X(3))];\r
 105 end\r
 106 \r
 107 % sortie des résultats sous forme symbolique\r
 108 % Libérer du commentaire, si souhaité, les deux lignes suivantes.\r
 109 % Voir aussi la fonction affiche_racines.\r
 110 \r
 111 disp('Sortie des résultats sous forme symbolique proposée par matlab');\r
 112 symb=sym(rac,'r');disp('rac = ');disp (symb);\r
 113 \r
 114 % fin de fonction