i found some sources where algorithms are suggested to determine complex roots in general analytic functions (no polynomials). E.g. "A Fortran program for solving a nonlinear equation by Muller's method", written by I. Barrodale and K. B. Wilson. Unfortunately the code is only available in pdf format and the transfer to ASCII is quite cumbersome... And also Henning Bach, Volume 121, Number 12, page 675-678, December 1969. Its available in ASCII. The letter one works fine but is quite slow for my purpose. Is there a standard routine for computing complex roots in Fortran?
Thanks :)
For the ones who are interested. I took the one from I. Barrodale and K. B. Wilson. Since the article is free to read online, I wanna post the code here. It is by far the fastest and most reliable I found. To use it, define a subroutine calcf(x,f). f is the function value for x. x and f are complex.
for more Information -> http://www.sciencedirect.com/science/article/pii/0771050X78900414
This is the version for double precision:
SUBROUTINE ROOTS(KNOWN,NMORE,KM,REALRT,EP1,EP2,GUESS,MAXITS,RTS,NFOUND, FNV, IFAIL)
COMPLEX*16 RTS(KM), FNV(KM), X1, X2, X3, F1, F2, F3, &
X21, X32, F21, F32, F31, F321, B, XNEW, FNEW, DENOM, &
FSAVE, CSQRT, RADICL, RT, DIF, CSTEP, FSLAST
REAL*8 EP1, EP2, ONE, HUNDRD, TWO, FOUR, HALF, STEP, &
ZERO, EP1DEF, EP2DEF
INTEGER GUESS
LOGICAL REALRT
DATA IZERO,IONE,ITWO,ITHREE,ZERO,HALF,ONE,TWO,FOUR,HUNDRD,ITSDEF,EP1DEF,EP2DEF &
/0,1,2,3,0.0d0,0.5d0,1.0d0,2.0d0,4.0d0,100.0d0,100,0.5d-5,1.0d-6/
EP1 = AMAX1(EP1,EP1DEF)
EP2 = AMAX1(EP2,EP2DEF)
MAXITS = MIN0(MAXITS,ITSDEF)
IFAIL = IZERO
NFOUND = IZERO
IF (KNOWN.LT.IONE) GO TO 30
DO 10 I=1,KNOWN
II = I
CALL CALCF(RTS(I), FNV(I))
IF (ABS(FNV(II)).GE.EP2) GO TO 20
10 CONTINUE
20 IFAIL =IONE
GUESS = GUESS + KNOWN - II +IONE
NMORE = NMORE + KNOWN - II +IONE
KNOWN = II - IONE
30 CONTINUE
LOOP1 = KNOWN +IONE
LOOP2 = KNOWN + NMORE
IF (LOOP1.GT.LOOP2 .OR. LOOP1.LT.IONE) GO TO 180
IF (GUESS-NMORE) 40, 70, 60
40 ILO = GUESS + IONE
DO 50 I=ILO,LOOP2
RTS(I) = ZERO
50 CONTINUE
GO TO 70
60 CONTINUE
GUESS = NMORE
70 CONTINUE
STEP = HALF
DO 160 NEW=LOOP1,LOOP2
KOUNT = ITHREE
NEWM1 = NEW - IONE
RT = RTS(NEW)
X1 = RT - STEP
X2 = RT + STEP
X3 = RT
CALL CALCF(X1, F1)
CALL CALCF(X2, F2)
CALL CALCF(X3, F3)
FSLAST = F3
IF (NEW.GT.IONE) CALL TEST(X1, F1, FSAVE, RTS,NEWM1, EP2, KOUNT)
IF (NEW.GT.IONE) CALL TEST(X2, F2, FSAVE, RTS,NEWM1, EP2, KOUNT)
IF (NEW.GT.IONE) CALL TEST(X3, F3, FSLAST, RTS,NEWM1, EP2, KOUNT)
F21 = (F2-F1)/(X2-X1)
80 X32 = X3 - X2
F32 = (F3-F2)/X32
F321 = (F32-F21)/(X3-X1)
B = F32 + X32*F321
RADICL = B**ITWO - FOUR*F321*F3
IF (REALRT .AND. REAL(RADICL).LT.ZERO) RADICL = ZERO
RADICL = SQRT(RADICL)
IF (REAL(B)*REAL(RADICL)+AIMAG(B)*AIMAG(RADICL).LT.ZERO) RADICL = -RADICL
DENOM = B + RADICL
IF (ABS(DENOM).NE.ZERO) GO TO 100
IF (ABS(F3).GE.EP2) GO TO 90
XNEW = X3
GO TO 120
90 XNEW = X3 + X32
GO TO 120
100 CSTEP = TWO*F3/DENOM
IF (.NOT.REALRT .OR. ABS(F3).EQ.ZERO .OR. ABS(F32).EQ.ZERO) GO TO 110
CSTEP = F32/ABS(F32)*F3/ABS(F3)*ABS(CSTEP)
110 XNEW = X3 - CSTEP
120 CALL CALCF(XNEW, FNEW)
FSAVE = FNEW
IF (NEW.LE.IONE) GO TO 130
CALL TEST(XNEW, FNEW, FSAVE, RTS, NEWM1, EP2,KOUNT)
130 KOUNT = KOUNT +IONE
IF (KOUNT.GT.MAXITS) GO TO 170
DIF = XNEW - X3
IF (ABS(DIF).LT.EP1*ABS(XNEW)) GO TO 150
IF (ABS(FSAVE).LT.EP2) GO TO 150
IF (REALRT .OR. ABS(FSAVE).LT.HUNDRD*ABS(FSLAST)) GO TO 140
CSTEP = CSTEP*HALF
XNEW = XNEM + CSTEP
GO TO 120
140 X1 = X2
X2 = X3
X3 = XNEW
F1 = F2
F2 = F3
F3 = FNEW
FSLAST = FSAVE
F21 = F32
GO TO 80
150 CONTINUE
RTS(NEW) = XNEW
FNV(NEW) = FSAVE
NFOUND = NEW
160 CONTINUE
GO TO 190
170 CONTINUE
RTS(NEW) = XNEW
FNV(NEW) = FSAVE
IFAIL = ITHREE
GO TO 190
180 CONTINUE
IFAIL = ITWO
190 CONTINUE
RETURN
END
SUBROUTINE TEST(X, F, FSAVE, RTS, K, EPS, KOUNT)
COMPLEX*16 X, F, RTS(K), D, FSAVE
REAL*8 EPS!, CABS
DATA IONE, PERTB /1,0.01d0/
10 CONTINUE
DO 20 I=1,K
D = X - RTS(I)
IF (ABS(D).LE.EPS) GO TO 30
F = F/D
20 CONTINUE
GO TO 40
30 CONTINUE
X = X + PERTB
CALL CALCF(X, F)
FSAVE = F
KOUNT = KOUNT +IONE
GO TO 10
40 RETURN
END