tel/TPmatlab/equation_nonlineaire/TP4a/iteration_newton.m

   1 function [n,X]=iteration_newton(x0,nmax,epsilon,fcn)\r
   2 \r
   3 \r
   4 %       iteration_newton : calcul d'une racine la méthode de Newton\r
   5 %\r
   6 % *********************************************************\r
   7 %\r
   8 %       [n,X]=iteration_newton(x0,nmax,epsilon,fcn) renvoie\r
   9 %  les itérés de la méthode de Newton\r
  10 %  (la dérivée f' est calculée symboliquement)\r
  11 %\r
  12 %       variables d'entrées : \r
  13 %               * nmax : nombre maximal d'itérations\r
  14 %               * x0 : le  premier terme\r
  15 %     * epsilon : précision souhaitée\r
  16 %     * fcn est une chaîne de caractère représentant la fonction \r
  17 %      (de type inline, builtin ou par fichier M-file);\r
  18 %\r
  19 %       variables de sortie :\r
  20 %     * n est l'indice correspondant à x_n avec f(x_n)=0 ou |x_n-x_{n-1}|<epsilon \r
  21 %                                         ou n>=nmax\r
  22 %     * X est une vecteur qui contient les valeurs de x_0 à x_n\r
  23 %\r
  24 % \r
  25 % \r
  26 % ************ Fonctions auxiliaires utilisées ************\r
  27 %\r
  28 %       aucune\r
  29 %\r
  30 % *********************************************************\r
  31 %\r
  32 \r
  33 \r
  34 \r
  35 % Contrôles d'entrée\r
  36 qs=exist(fcn);\r
  37 if (qs~=2) & (qs~=5) & (qs~=1)\r
  38    error('fcn doit être le nom d''une fonction (built-in, M-file ou inline)');\r
  39 end\r
  40 \r
  41 % corps d'algorithme\r
  42 syms vx\r
  43 g=feval(fcn,vx);\r
  44 fcnpsymb=diff(g,vx);\r
  45 n=0;\r
  46 test=(feval(fcn,x0)~=0);\r
  47 x=x0;\r
  48 fx=feval(fcn,x);\r
  49 fxp=subs(fcnpsymb,vx,x);\r
  50 if ~(isnumeric(fxp));\r
  51    fxp=eval(fxp);\r
  52 end\r
  53 X=x0;\r
  54 while(test)\r
  55    xnp1=x-fx/fxp;\r
  56    n=n+1;\r
  57    fxnp1=feval(fcn,xnp1);\r
  58    test=~((fxnp1==0)|(n>=nmax)|(abs(x-xnp1)<=epsilon));\r
  59    x=xnp1;\r
  60    fx=feval(fcn,x);\r
  61    fxp=subs(fcnpsymb,vx,x);\r
  62    if ~(isnumeric(fxp));\r
  63       fxp=eval(fxp);\r
  64    end\r
  65    X=[X x];\r
  66 end   \r
  67 if (n>=nmax)\r
  68    disp('attention, n>nmax');\r
  69 end\r
  70 \r
  71 \r