Procédures | Cast3M

* FZERO     PROCEDUR  BP208322  19/02/22    21:15:01     10118          
 
************************************************************************
*   __Objet__ :
*           TROUVER LE ZERO D'UNE FONCTION f(x)
*
*   __Syntaxe__ : 
*           XSOL = FZERO FONCTION XA XB (XTOL);
*   avec :
*   * FONCTION : procedur definisant la fonction x -> y=f(x) par :
*                y = FONCTION x;
*   * [XA XB]  : intervalle de recherche
*   * XTOL (facultatif) : précision sur la solution
*
*   __Référence__ :
*                   METHODE DE BRENT, 1973
*   De nombreuses mises en oeuvre numeriques existent : fortran, c++,
*   matlab, ... et gibiane !
*
*   __Exemple__ :
*   debp f6 x;
*   y  = (exp  x) - 2. - (1. / ((10.*x)**2)) + (2./((100.*x)**3));
*   finp y;
*   x6 = FZERO f6 -4. 2.;
*   y6 = f6 x6;
*   * on obtient : x6=0.70320 et f6=-2.7E-16
*   * verification graphique :
*   x6_p = prog -4. PAS 0.01 2.;
*   f6_p = f6 x6_p;
*   ev6 = evol bleu manu 'x' x6_p 'f(x)' f6_p;
*   dess ev6 'YBOR' -5 5 'YGRA' 1 'AXES';
*
************************************************************************
 
DEBP FZERO F*'PROCEDUR' XA*'FLOTTANT' XB*'FLOTTANT';
 
* option verbeux ?
  BLABLA = (VALE IMPI) >EG 1;
 
* recup des donnees d'entrees + valeurs par defaut
  FA = F XA;
  FB = F XB;
  FTOL = 1.E-14 * (MAXI 'ABS' (PROG FA FB 1.E-32));
  ARGU XTOL/'FLOTTANT';
  SI (NON (EXIS XTOL)); XTOL=1.E-14*(ABS (XB - XA)); FINSI;
  SI BLABLA; MESS 'tolerance fixee a ' XTOL; FINSI;
* on utlise la convergnce de la dichotomie pour calculer NITER
  NITER = ENTIER SUPERIEUR (-1.*(LOG XTOL) / (LOG 2.));
 
* sauvegarde  
  XASAVE=XA; XBSAVE=XB;
  FASAVE=FA; FBSAVE=FB;
 
* tests preliminaires   
  SI ((ABS FA) < FTOL);  XB=XA; FB=FA; QUIT FZERO;  FINSI;  
  FAFB = FA * FB ;
  SI (FAFB > 0.); 
    mess 'Pas de changement de signe detecte!'; 
    erre 21; quit FZERO;
  FINSI;
 
  FC=FB;
  ialgo=MOT 'initial';
* ---------------- on itere ----------------
  REPE BITER NITER;
 
*   test de convergence 1
    SI ( ((ABS (XB - XA)) &lt;EG XTOL ) OU ((ABS FB) < FTOL) );
      QUIT BITER;
    finsi;
 
*   On s'assure que : 
*   B=meilleur resultat, A=precedente valeur de B
*   C=resultat de signe oppose a B
    FBFC=FB*FC;
    SI (FBFC > 0.);
        XC = XA;  FC = FA;
        XD = XB - XA;  XE = XD;
    FINSI;
    SI ((ABS FC) < (ABS FB));
        XA = XB;  XB = XC;  XC = XA;
        FA = FB;  FB = FC;  FC = FA;
    FINSI;
 
*   test de convergence 2
*   M=demi-intervalle
    XM = 0.5*(XC - XB);
    SI((ABS XB) > 1.); XTOLER = 2.*XTOL*(ABS XB);
    SINON;             XTOLER = 2.*XTOL;
    FINSI;
    SI ( ((ABS XM) &lt;EG XTOLER) OU ((ABS FB) < FTOL) );
      QUIT BITER;
    FINSI;     
    SI BLABLA; mess (&BITER + 1) ': x=' XB' -> F='FB ' ('ialgo')'; FINS;
 
*   choix de la "meilleure" stratégie
    SI ( ((ABS XE) < XTOLER) OU ((ABS FA) &lt;EG (ABS FB)) );
*     methode de la bisection (dichotomie)
      XD = XM;  XE = XM;
      ialgo = mot 'dichotomie';  
    SINON;
*       Interpolation ...
        S = FB/FA;
        SI (XA EGA XC);
*           ... lineaire
            P = 2.0*XM*S;
            Q = 1.0 - S;
            ialgo = mot 'interp. lin.';  
        SINON;
*           ...quadratique inverse
            Q = FA/FC;
            R = FB/FC;
            P = (S*(2.0*XM*Q*(Q - R)) - ((XB - XA)*(R - 1.0)));
            Q = (Q - 1.0)*(R - 1.0)*(S - 1.0);
            ialgo = mot 'interp. quad.';  
        FINSI;;
        SI (P > 0.);  Q = -1. * Q;
        SINON;        P = -1. * P;
        FINSI;
*       point issu de l'interpolation acceptable ?
        SI ( ((2.0*P) < ((3.0*XM*Q) - (ABS (XTOLER*Q)))) 
          ET (P < (ABS (0.5*XE*Q))) );
            XE = XD;  XD = P/Q;
        SINON;
            XD = XM;  XE = XM;
            ialgo='dichotomie';
        FINSI;
    FINSI; 
 
*   point suivant
    XA = XB; FA = FB;
    SI ((ABS XD) > XTOLER);
        XB = XB + XD;
    SINON; 
        SI(XB > XC); XB = XB - XTOLER; 
        SINON;       XB = XB + XTOLER;
        FINSI;        
    FINSI;
    FB = F XB;
 
  FIN BITER;
 
  SI BLABLA; mess &BITER ': x=' XB' -> F='FB ' ('ialgo')'; FINSI;
 
FINP XB;