tel/TPmatlab/equation_nonlineaire/TP4f/iteration_ordre3_gen.m

   1 function [n,X]=iteration_ordre3_gen(x0,nmax,epsilon,fcn)\r
   2 \r
   3 \r
   4 %  iteration_ordre3_gen : calcul d'une racine la méthode d'ordre 3 générale\r
   5 %\r
   6 % *********************************************************\r
   7 %\r
   8 %  [n,X]=iteration_ordre3_gen(x0,nmax,epsilon,fcn) renvoie\r
   9 %  l'approximation de la racine d'une fonction par la méthode d'ordre trois\r
  10 %  générale (les dérivées f' et f'' sont calculées 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 %     * X est une vecteur qui contient les valeurs x_0, x_1, ...,x_n\r
  21 %     * n est l'indice correspondant à x_n avec f(x_n)=0 ou |x_n-x_{n-1}|<epsilon \r
  22 %      ou n>=nmax\r
  23 % \r
  24 % \r
  25 % ************ Fonctions auxiliaires utilisées ************\r
  26 %\r
  27 %       aucune\r
  28 %\r
  29 % *********************************************************\r
  30 %\r
  31 \r
  32 \r
  33 % Contrôles d'entrée\r
  34 qs=exist(fcn);\r
  35 if (qs~=2) & (qs~=5) & (qs~=1)\r
  36    error('fcn doit être le nom d''une fonction (built-in, M-file ou inline)');\r
  37 end\r
  38 \r
  39 % corps d'algorithme\r
  40 syms vx\r
  41 g=feval(fcn,vx);\r
  42 fcnpsymb=diff(g,vx);\r
  43 fcnppsymb=diff(fcnpsymb,vx);\r
  44 n=0;\r
  45 test=(feval(fcn,x0)~=0);\r
  46 x=x0;\r
  47 X=x;\r
  48 fx=feval(fcn,x);\r
  49 fxp=subs(fcnpsymb,vx,x);\r
  50 if ~(isnumeric(fxp));\r
  51    xpp=eval(fxp);\r
  52 end\r
  53 fxpp=subs(fcnppsymb,vx,x);\r
  54 if ~(isnumeric(fxpp));\r
  55    xpp=eval(fxpp);\r
  56 end\r
  57 while(test)\r
  58    A=fxpp/2;\r
  59    B=fxp;\r
  60    C=fx;\r
  61    if (C==0)\r
  62       d=0;    \r
  63    else   \r
  64       r=sqrt(B^2-4*A*C);\r
  65       Bp=B+r;\r
  66       Bm=B-r;\r
  67       Bpn=abs(Bp);\r
  68       Bmn=abs(Bm);\r
  69       if (Bpn<Bmn)\r
  70          d=-2*C/Bm;     \r
  71       else\r
  72          d=-2*C/Bp;     \r
  73       end\r
  74    end\r
  75    xnp1=x+d;\r
  76    fxnp1=feval(fcn,xnp1);\r
  77    test=~((fxnp1==0)|(n>=nmax)|(abs(x-xnp1)<=epsilon));\r
  78    x=xnp1;\r
  79    fx=feval(fcn,x);\r
  80    fxp=subs(fcnpsymb,vx,x);\r
  81    if ~(isnumeric(fxp));\r
  82       xpp=eval(fxp);\r
  83    end\r
  84    fxpp=subs(fcnppsymb,vx,x);\r
  85    if ~(isnumeric(fxpp));\r
  86       xpp=eval(fxpp);\r
  87     end\r
  88    n=n+1;\r
  89    X=[X x];\r
  90 end   \r
  91 if (n>=nmax)\r
  92    disp('attention, n>nmax');\r
  93 end\r