Download zernew.eso

Back to the list

   1 : C ZERNEW    SOURCE    CHAT      05/01/13    04:22:41     5004           
   2 :       subroutine zernew(xr,xi,coef,m,xvr,xvi,kerre)
   3 :       IMPLICIT INTEGER(I-N)
   4 :       IMPLICIT REAL*8(A-H,O-Z)
   5 :       dimension coef(*)
   6 : * coef(1),coef(2),...coef(m+1)  sonrt les coef (reels) du polynome de 
   7 : * degré m .coef(1) est le terme constant
   8 : * xr et xi sont les parties reelle et imaginaire d'une racine approchée
   9 : * et on cherche la solution excate par un newton
  10 : * en sortie xvr et xvi sont les parties reelle et imaginaire de la
  11 : * solution
  12 :       eps=1.d-10
  13 :       xvr=xr
  14 :       xvi=xi
  15 :       icou=0
  16 :  1    continue
  17 :       icou=icou+1
  18 : * evaluation du polynome et de sa dérivée
  19 :       xpor=coef(1)
  20 :       xpoi=0.d0
  21 :       xder=coef(2)
  22 :       xdei=0.d0
  23 :       xpur=xvr
  24 :       xpui=xvi
  25 :       do 10 ll=1,m-1
  26 :       xpor=xpor+xpur*coef(1+ll)
  27 :       xpoi=xpoi+xpui*coef(1+ll)
  28 :       xder=xder+xpur*coef(2+ll)*(ll+1)
  29 :       xdei=xdei+xpui*coef(2+ll)*(ll+1)
  30 :       xpurr = xpur*xvr-xpui*xvi
  31 :       xpuii=xpur*xvi+xpui*xvr
  32 :       xpur=xpurr
  33 :       xpui=xpuii
  34 :    10 continue
  35 :       xpor=xpor+xpur*coef(1+m)
  36 :       xpoi=xpoi+xpui*coef(1+m)
  37 : *      write(6,*) 'xpor,xpoi,xder,xdei',xpor,xpoi,xder,xdei
  38 : * recherche de xdr et xdi
  39 : 
  40 :        denom=xder*xder+xdei*xdei
  41 :       if(denom.eq.0.d0) then
  42 : *        write(6,*) 'dérivée nulle'
  43 :          kerre=22
  44 :         return
  45 :       endif 
  46 :       xdr=-(xpor*xder +xpoi*xdei)/denom
  47 :       xdi= (xpor*xdei-xpoi*xder)/denom
  48 :       xvr=xvr+xdr
  49 :       xvi = xvi+xdi
  50 :        xcomp=max(abs(xvr),abs(xvi))*eps
  51 : *       write(6,*) 'denom,xdr,xdi,xvr,xvi',denom,xdr,xdi,xvr,xvi
  52 :       if(abs(xdr).le.xcomp.and.abs(xdi).le.xcomp.and.icou.ge.4)
  53 :      $      go to 2
  54 :       go to 1
  55 :     2 continue
  56 : *      write(6,*) ' icou',icou
  57 :       return
  58 :       end
  59 : 
  60 : 
  61 : 
  62 :  
  63 :  
  64 :  

© Cast3M 2003 - All rights reserved.
Disclaimer