'**************************************************** '* Find all roots of a complex polynomial using * '* Newton's iterative formulation. * '* ------------------------------------------------ * '* SAMPLE RUN: * '* (Find roots of complex polynomial: * '* (2.5 + I) X3 + (7.5 - 12 I) X2 + (-3.75 + 0.4 I * '* ) X + (4.25 - I) ) * '* * '* Error code = 0 * '* * '* Complex roots are: * '* 0.25996740 -0.39966139 * '* -0.11184574 0.66800696 * '* -1.07915614 4.90406822 * '* * '* ------------------------------------------------ * '* Ref.: From Numath Library By Tuan Dang Trong in * '* Fortran 77 [BIBLI 18]. * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '**************************************************** DEFDBL A-H, O-Z DEFINT I-N N = 3 'order of polynomial DIM A(N + 1), B(N + 1), X(N), Y(N) ' A: real parts of complex coefficients ' B: imaginary parts of complex coefficients ' (stored in decreasing powers) ' X: real parts of roots ' Y: imaginary parts of roots */ DIM W1(N + 1), W2(N + 1), W3(N + 1), W4(N + 1) ' working zones 'define vectors A, B of coefficients A(1) = 2.5: B(1) = 1# A(2) = 7.5: B(2) = -12# A(3) = -3.75: B(3) = .4 A(4) = 4.25: B(4) = -1# ' Example #2 (real roots are -3, 1, 2) ' A(1)= 1. B(1)=0.0 ' A(2)= 0. B(2)=0.0 ' A(3)=-7. B(3)=0.0 ' A(4)= 6. B(4)=0.0 */ ' call main function GOSUB 1000 'call CNEWTON(N,A,B,RR,RI,IER,W1,W2,W3,W4) F\$ = "##.########" ' print results CLS PRINT PRINT " Error code = "; IER PRINT PRINT " Complex roots are:" FOR I = 1 TO N PRINT " "; PRINT USING F\$; X(I); PRINT " "; PRINT USING F\$; Y(I) NEXT I PRINT END 1000 'CNEWTON (N, A, B, X, Y, IER, W1, W2, W3, W4) '--------------------------------------------------------------------- ' CALCULATE ALL THE ROOTS OF A COMPLEX POLYNOMIAL USING ' NEWTON'S ITERATIVE FORMULATION ' ' INPUTS: ' N ORDER OF POLYNOMIAL ' A,B TABLES OF DIMENSION N+1, RESPECTIVELY REAL, ' IMAGINARY PART OF POLYNOMIAL COMPLEX COEFFICIENTS ' C = A+I*B, STORED IN DECREASING POWERS. ' OUTPUTS: ' X,Y TABLES OF DIMENSION N, RESPECTIVELY REAL, ' IMAGINARY PART OF COMPLEX ROOTS FOUND, Z = X+I*Y ' IER ERROR CODE = 0, CONVERGENCE OK ' = I, NO CONVERGENCE FOR I-TH ROOT ' WORKING ZONES: ' G,H,T,F TABLES OF DIMENSION N '---------------------------------------------------------------------- ' REFERENCE: ' E.DURAND -SOLUTIONS NUMERIQUES DES EQUATIONS ALGEBRIQUES, ' TOME 1, MASSON & CIE, 1960, PP.260-266 '---------------------------------------------------------------------- ' Labels: 5, 7, 10, 15, 30, 40, 45, 55, 60, 70, 75 ' int I,IT1,IT2,K,K1,K2,L,M ' double EPS,EPS1,P,Q,R,TOL,X0,X1,Y0,Y1 ' double CABS,DX,DY,GD,GX,GY,R1,R2,R3,X2,Y2 IT1 = 25: IT2 = 100: TOL = .01: EPS = 0# IER = 0 IF EPS <> 0! THEN GOTO 7 EPS = 1# 5 EPS = .5 * EPS EPS1 = EPS + 1# IF EPS1 > 1# THEN GOTO 5 7 L = N X0 = 1# Y0 = 1# W1(1) = A(1) W3(1) = A(1) W2(1) = B(1) W4(1) = B(1) 10 IF L = 1 THEN GOTO 55 M = L + 1 I = N - L + 1 K1 = 0 K2 = 0 X1 = X0 Y1 = Y0 15 FOR K = 2 TO M W1(K) = A(K) + X1 * W1(K - 1) - Y1 * W2(K - 1) W2(K) = B(K) + Y1 * W1(K - 1) + X1 * W2(K - 1) NEXT K FOR K = 2 TO L W3(K) = W1(K) + X1 * W3(K - 1) - Y1 * W4(K - 1) W4(K) = W2(K) + Y1 * W3(K - 1) + X1 * W4(K - 1) NEXT K GD = W3(L) * W3(L) + W4(L) * W4(L) GX = W1(M) * W3(L) + W2(M) * W4(L) GY = W2(M) * W3(L) - W1(M) * W4(L) DX = -GX / GD DY = -GY / GD X2 = X1 + DX Y2 = Y1 + DY R1 = (ABS(DX) + ABS(DY)) / (ABS(X2) + ABS(Y2)) K2 = K2 + 1 X1 = X2 Y1 = Y2 IF R1 < TOL THEN GOTO 30 IF K2 <= IT2 THEN GOTO 15 IER = I GOTO 60 30 K1 = K1 + 1 IF K1 <= IT1 THEN GOTO 15 X(I) = X2 Y(I) = Y2 IF ABS(X2) <= EPS THEN X2 = EPS R2 = ABS(Y2 / X2) IF R2 < .1 THEN GOTO 40 IF R2 * EPS > 1# THEN X(I) = 0# X0 = X2 Y0 = Y2 GOTO 45 40 IF R2 <= EPS THEN Y(I) = 0# X0 = X2 + Y2 Y0 = X2 - Y2 45 FOR K = 1 TO L A(K) = W1(K) B(K) = W2(K) NEXT K L = L - 1 GOTO 10 55 R3 = W1(1) * W1(1) + W2(1) * W2(1) X(N) = -(W1(1) * W1(2) + W2(1) * W2(2)) / R3 Y(N) = (W2(1) * W1(2) - W1(1) * W2(2)) / R3 IF ABS(X(N)) <= EPS THEN X(N) = EPS R2 = ABS(Y(N) / X(N)) IF R2 * EPS > 1# THEN X(N) = 0# IF R2 <= EPS THEN Y(N) = 0# 60 FOR L = 1 TO N P = X(L) Q = Y(L) R = P * P + Q * Q IF L = 1 THEN GOTO 70 FOR I = L TO 2 STEP -1 CABS = X(I - 1) * X(I - 1) + Y(I - 1) * Y(I - 1) IF R >= CABS THEN GOTO 75 X(I) = X(I - 1) Y(I) = Y(I - 1) NEXT I 70 I = 1 75 X(I) = P Y(I) = Q NEXT L RETURN 'CNEWTON ' end of file tnewton.bas