subroutine nuraci (rfonc,rtab,ra,rb,rh,repsi,rsol,inbr,ier) c c*********************************************************************** c* c* "Copyright [c] CNES 98 - tous droits reserves" c* ********************************************** c* c*PRO MAGLIB c* c*VER 99.03.31 - V 1.0 c*VER 01.06.05 - V 2.0 c*VER 03.01.06 - V 2.1 c* c*AUT spec. CNES - BERGES - mars 1986 c*AUT (c) MSLIB 0029 24.09.87 v1.0 c*AUT port. CISI c* c*ROL Theme : Calculs mathematiques c*ROL Resolution equation non lineaire sur [a,b]. c*ROL Recherche d'une racine reelle de l'equation non lineaire c*ROL f(x) = h sur l'intervalle [a,b] par la methode de la c*ROL secante avec une acceleration de la convergence par c*ROL dichotomie. c* c*PAR rfonc (I) : nom de la fonction calculant f(x) c*PAR : (definie par appelant) c* c*PAR rtab (I) : tableau des parametres definissant f(x) c* c*PAR ra (I) : borne a de l'intervalle [a,b] c*PAR rb (I) : borne b de l'intervalle [a,b] c* c*PAR rh (O) : valeur du second membre de l'equation f(x) = h c* c*PAR repsi (O) : seuil de precision demande c* c*PAR rsol (O) : valeur de la racine c* c*PAR inbr (O) : nombre de racines sur [a,b] c* c*PAR ier (O) : code de retour c* c*NOT ier : sans objet c* c*NOT inbr : 0 : 0 ou 2 * n racines sur [a,b] c*NOT inbr : (par convention rsol=0) c*NOT inbr : 1 : une racine c*NOT inbr : 2 : a et b sont racines (par convention rsol = a) c* c*INF utilise : rfonc c* c*HST version 1.0 - 99.03.31 - creation de la maglib au CDPP c*HST version 2.0 - 01.06.05 - correction de commentaires de code c*HST version 2.1 - 03.01.06 - corrections en compilation avec g77 c* c*********************************************************************** c* implicit none c c --------------------------------- c*FON Declaration identificateur rcs_id c --------------------------------- c character rcs_id*100 c c -------------------------- c*FON Declaration des parametres c -------------------------- c double precision rtab(*),ra,rb,rh,repsi,rsol integer inbr integer ier c c ---------------------------------- c*FON Declaration des fonctions externes c ---------------------------------- c external rfonc double precision rfonc c c --------------------------------- c*FON Declaration des variables locales c --------------------------------- c integer i1,i2 c*LOC i1,i2 : variables de controle c double precision y,y1,y2,s,s1,s2 c*LOC y,y1,y2,s,s1,s2 : valeurs de la fonction et signes c double precision x,x1,x2 c*LOC x,x1,x2 : valeurs pour la recherche de la racine c double precision epsdp c*LOC epsdp : epsilon sun double precision c double precision rmin c*LOC rmin : plus petit reel positif c double precision eps100 c*LOC eps100 : epsilon sun double precision * 100 c double precision r c*LOC r : valeur a comparer a epsilon c SAVE c c --------------------------------- c*FON Affectation identificateur rcs_id c --------------------------------- c data rcs_id /" >$Id$"/ c c ****************** c Debut de programme c ****************** c c ------------------------------------------------- c*FON Initialisation precision pour le calcul (epsilon) c ------------------------------------------------- c epsdp = 0.2220446d-15 eps100 = 100.d0 * epsdp rmin = eps100 eps100 = max(repsi,eps100) ier = 0 c c --------------------------------------------------------- c*FON Initialisations c*FON x1 = borne a de l'intervalle de recherche et x2 = borne b c --------------------------------------------------------- c i1 = 0 i2 = 0 x1 = ra x2 = rb c c ------------------------------------------ c*FON y1 : valeur de la fonction f(x) - rh en x1 c*FON y2 : valeur de la fonction f(x) - rh en x2 c ------------------------------------------ c y1 = rfonc(x1,rtab,ier) - rh y2 = rfonc(x2,rtab,ier) - rh c c ------------------------------------------- c*FON Si y1 = 0 alors x1 est racine de l'equation c ------------------------------------------- c if (abs(y1) .lt. rmin) then inbr = 1 if (abs(y2) .lt. rmin) inbr = 2 rsol = x1 go to 30 c c ------------------------------------------------- c*FON Sinon si y2 = 0 alors x2 est racine de l'equation c ------------------------------------------------- c elseif (abs(y2) .lt. rmin) then inbr = 1 rsol = x2 go to 30 endif c c c -------------------------------------------- c*FON Si y1 et y2 de meme signe 0 ou 2 * n racines c --------------------------------------------- c s1 = sign(1.d0,y1) s2 = sign(1.d0,y2) c if (s1 * s2 .gt. 0.d0) then rsol = 0.d0 inbr = 0 go to 30 endif c c -------------------------------------------------- c*FON Recherche de la racine par iteration : c*FON on procede par dichotomie sur l'intervalle [x1,x2] c -------------------------------------------------- c c evaluation de la valeur x c 10 continue x = (x1 * y2 - x2 * y1) / (y2 - y1) c 20 continue c c*FON calcul de la fonction en x (valeur comprise entre x1 et x2) c y = rfonc(x,rtab,ier) - rh c c ------------------------------------------ c*FON Si rapport < epsilon on a trouve la racine c ------------------------------------------ c r = abs(x2 - x1) / (1.d0 + abs(x1)) if (r .lt. eps100) then inbr = 1 rsol = x go to 30 endif c c ------------------------------------------------------------ c*FON Si y1 et y de signe oppose recherche sur l'intervalle [x1,x] c ------------------------------------------------------------ c s1 = sign(1.d0,y1) s = sign(1.d0,y) if (s1 * s .lt. 0.d0) then x2 = x y2 = y i1 = i1 + 1 i2 = 0 c c acceleration par dichotomie if (i1 .ne. 2) then if (i1 .gt. 2) i1 = 0 go to 10 else x = (x1 + x2) / 2.d0 go to 20 endif c c --------------------------------- c*FON Sinon si y = 0 alors x est racine c --------------------------------- c else if (abs(y) .lt. rmin) then inbr = 1 rsol = x go to 30 c else c c* ------------------------------------------------------------- c*FON Sinon y1 et y de meme signe recherche sur l'intervalle [x,x2] c* ------------------------------------------------------------- c x1 = x y1 = y i2 = i2 + 1 i1 = 0 c c acceleration par dichotomie c if (i2 .ne. 2) then if (i2 .gt. 2) i2 = 0 go to 10 else x = (x1 + x2) / 2.d0 go to 20 endif endif c c **************** c Fin de programme c **************** c 30 continue return end