Télécharger fzero.procedur

Retour à la liste

Numérotation des lignes :

  1. * FZERO PROCEDUR BP208322 19/02/22 21:15:01 10118
  2.  
  3. ************************************************************************
  4. * __Objet__ :
  5. * TROUVER LE ZERO D'UNE FONCTION f(x)
  6. *
  7. * __Syntaxe__ :
  8. * XSOL = FZERO FONCTION XA XB (XTOL);
  9. * avec :
  10. * * FONCTION : procedur definisant la fonction x -> y=f(x) par :
  11. * y = FONCTION x;
  12. * * [XA XB] : intervalle de recherche
  13. * * XTOL (facultatif) : précision sur la solution
  14. *
  15. * __Référence__ :
  16. * METHODE DE BRENT, 1973
  17. * De nombreuses mises en oeuvre numeriques existent : fortran, c++,
  18. * matlab, ... et gibiane !
  19. *
  20. * __Exemple__ :
  21. * debp f6 x;
  22. * y = (exp x) - 2. - (1. / ((10.*x)**2)) + (2./((100.*x)**3));
  23. * finp y;
  24. * x6 = FZERO f6 -4. 2.;
  25. * y6 = f6 x6;
  26. * * on obtient : x6=0.70320 et f6=-2.7E-16
  27. * * verification graphique :
  28. * x6_p = prog -4. PAS 0.01 2.;
  29. * f6_p = f6 x6_p;
  30. * ev6 = evol bleu manu 'x' x6_p 'f(x)' f6_p;
  31. * dess ev6 'YBOR' -5 5 'YGRA' 1 'AXES';
  32. *
  33. ************************************************************************
  34.  
  35. DEBP FZERO F*'PROCEDUR' XA*'FLOTTANT' XB*'FLOTTANT';
  36.  
  37. * option verbeux ?
  38. BLABLA = (VALE IMPI) >EG 1;
  39.  
  40. * recup des donnees d'entrees + valeurs par defaut
  41. FA = F XA;
  42. FB = F XB;
  43. FTOL = 1.E-14 * (MAXI 'ABS' (PROG FA FB 1.E-32));
  44. ARGU XTOL/'FLOTTANT';
  45. SI (NON (EXIS XTOL)); XTOL=1.E-14*(ABS (XB - XA)); FINSI;
  46. SI BLABLA; MESS 'tolerance fixee a ' XTOL; FINSI;
  47. * on utlise la convergnce de la dichotomie pour calculer NITER
  48. NITER = ENTIER SUPERIEUR (-1.*(LOG XTOL) / (LOG 2.));
  49.  
  50. * sauvegarde
  51. XASAVE=XA; XBSAVE=XB;
  52. FASAVE=FA; FBSAVE=FB;
  53.  
  54. * tests preliminaires
  55. SI ((ABS FA) < FTOL); XB=XA; FB=FA; QUIT FZERO; FINSI;
  56. FAFB = FA * FB ;
  57. SI (FAFB > 0.);
  58. mess 'Pas de changement de signe detecte!';
  59. erre 21; quit FZERO;
  60. FINSI;
  61.  
  62. FC=FB;
  63. ialgo=MOT 'initial';
  64. * ---------------- on itere ----------------
  65. REPE BITER NITER;
  66.  
  67. * test de convergence 1
  68. SI ( ((ABS (XB - XA)) &lt;EG XTOL ) OU ((ABS FB) < FTOL) );
  69. QUIT BITER;
  70. finsi;
  71.  
  72. * On s'assure que :
  73. * B=meilleur resultat, A=precedente valeur de B
  74. * C=resultat de signe oppose a B
  75. FBFC=FB*FC;
  76. SI (FBFC > 0.);
  77. XC = XA; FC = FA;
  78. XD = XB - XA; XE = XD;
  79. FINSI;
  80. SI ((ABS FC) < (ABS FB));
  81. XA = XB; XB = XC; XC = XA;
  82. FA = FB; FB = FC; FC = FA;
  83. FINSI;
  84.  
  85. * test de convergence 2
  86. * M=demi-intervalle
  87. XM = 0.5*(XC - XB);
  88. SI((ABS XB) > 1.); XTOLER = 2.*XTOL*(ABS XB);
  89. SINON; XTOLER = 2.*XTOL;
  90. FINSI;
  91. SI ( ((ABS XM) &lt;EG XTOLER) OU ((ABS FB) < FTOL) );
  92. QUIT BITER;
  93. FINSI;
  94. SI BLABLA; mess (&BITER + 1) ': x=' XB' -> F='FB ' ('ialgo')'; FINS;
  95.  
  96. * choix de la "meilleure" stratégie
  97. SI ( ((ABS XE) < XTOLER) OU ((ABS FA) &lt;EG (ABS FB)) );
  98. * methode de la bisection (dichotomie)
  99. XD = XM; XE = XM;
  100. ialgo = mot 'dichotomie';
  101. SINON;
  102. * Interpolation ...
  103. S = FB/FA;
  104. SI (XA EGA XC);
  105. * ... lineaire
  106. P = 2.0*XM*S;
  107. Q = 1.0 - S;
  108. ialgo = mot 'interp. lin.';
  109. SINON;
  110. * ...quadratique inverse
  111. Q = FA/FC;
  112. R = FB/FC;
  113. P = (S*(2.0*XM*Q*(Q - R)) - ((XB - XA)*(R - 1.0)));
  114. Q = (Q - 1.0)*(R - 1.0)*(S - 1.0);
  115. ialgo = mot 'interp. quad.';
  116. FINSI;;
  117. SI (P > 0.); Q = -1. * Q;
  118. SINON; P = -1. * P;
  119. FINSI;
  120. * point issu de l'interpolation acceptable ?
  121. SI ( ((2.0*P) < ((3.0*XM*Q) - (ABS (XTOLER*Q))))
  122. ET (P < (ABS (0.5*XE*Q))) );
  123. XE = XD; XD = P/Q;
  124. SINON;
  125. XD = XM; XE = XM;
  126. ialgo='dichotomie';
  127. FINSI;
  128. FINSI;
  129.  
  130. * point suivant
  131. XA = XB; FA = FB;
  132. SI ((ABS XD) > XTOLER);
  133. XB = XB + XD;
  134. SINON;
  135. SI(XB > XC); XB = XB - XTOLER;
  136. SINON; XB = XB + XTOLER;
  137. FINSI;
  138. FINSI;
  139. FB = F XB;
  140.  
  141. FIN BITER;
  142.  
  143. SI BLABLA; mess &BITER ': x=' XB' -> F='FB ' ('ialgo')'; FINSI;
  144.  
  145. FINP XB;
  146.  
  147.  
  148.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales