11-29-2010, 06:29 PM

After all the discussions with my problems with the HP71B. Here is the program that solves the roots of a polynomial using the Abterth method. Edit the Data statements in lines 100, 120, and 170 to change the operational parameters and/or the polynomial coefficients.

10 REM POLYNOMIAL ROOTS USING ABERTH ALGORITHM BY NAMIR C SHAMMAS

20 REM NOVEMBER 30, 2010

30 OPTION BASE 0

40 DESTROY ALL

50 INTEGER M, N, I, J, K

60 COMPLEX X, P, D, S

70 DISP "REMINDER: DATA STMTS PROVIDE INPUT"

80 REM ===== START DATA INPUT SECTION ======

90 REM NEXT DATA IS MAX ITERS, TOLERANCE

100 DATA 100, 0.00001

110 REM NEXT DATA IS POLYNOMIAL ORDER

120 DATA 3

130 REM NEXT DATA IS SEQUENCE OF REAL PART, IMAGINARY PART OF

140 REM POLYNOMIAL COEFFICIENTS. START WITH COEFF. OF HIGHEST TERM

150 REM AND MOVE DOWN UNTIL YOU LIST CONSTANT TERM

160 REM DATA RLCOEF(N), IMCOEF(N), ... RLCOF(1), IMCOEF(1), RLCOEF(0), IMCOEF(0)

170 DATA 1, 0, 2, 0, 3, 0, 4, 0

180 REM ========= END INPUT SECTION ==========

190 READ M, T0, N

200 COMPLEX R(N),D1(N),C(N)

210 FOR I = N TO 0 STEP -1

220 READ X1,X2

230 C(I)=(X1,X2)

240 NEXT I

250 X = (0.4,0.9)

260 R(0)=(1,0)

270 FOR I = 1 TO N

280 R(I) = X*R(I-1)

290 NEXT I

300 FOR K = 1 TO M

310 FOR I = 1 TO N

320 CALL POLY(C(), N, R(I), P)

330 CALL DERIV(C(), N, R(I), D)

340 S = (0,0)

350 FOR J=1 TO N

360 IF I <> J THEN S = S + 1 / (R(I) - R(J))

370 NEXT J

380 D1(I) = 1 / (D / P - S)

390 R(I) = R(I) - D1(I)

400 NEXT I

410 REM TEST CONVERGENCE CRITERIA

420 C0 = 0

430 FOR I = 1 TO N

440 IF ABS(D1(I)) > T0 THEN C0 = 1

450 NEXT I

460 IF C0 = 0 THEN 'EXIT'

470 NEXT K

480 'EXIT':

490 FOR I = 1 TO N

500 DISP R(I)

510 PAUSE

520 NEXT I

530 END

540 SUB POLY(C(), N, X, S)

550 COMPLEX P

560 S = C(0)

570 P = X

580 FOR I = 1 TO N

590 S = S + C(I) * P

600 P = P * X

610 NEXT I

620 END SUB

630 SUB DERIV(C(), N, X, S)

640 COMPLEX P

650 S = C(1)

660 P = X

670 FOR I = 2 TO N

680 S = S + I * C(I) * P

690 P = P * X

700 NEXT I

710 END SUB

Enjoy!

Namir

*Edited: 30 Nov 2010, 2:45 p.m. after one or more responses were posted*