(71B) Aberth method
02-23-2022, 03:13 PM (This post was last modified: 02-23-2022 03:31 PM by robve.)
Post: #1
 robve Senior Member Posts: 379 Joined: Sep 2020
(71B) Aberth method
ABERTH POLYNOMIAL ROOT FINDER
Aberth method: Durand-Kerner with Newton's method
Estimate the complex roots of a polynomial with real coefficients
Note: Durand-Kerner is not guaranteed to converge

Optimized implementation:
- Finds complex roots z of polynomial p such that |p(z)|<T given tolerance T
- Prevents overflow/nan/inf by bounding the root shift adjustment magnitude during polishing
- Tested with a large sample of polynomials

For HP-71B + Math Pac

Degree: max polynomial degree limited by time and memory only
Levels: up to L=20 iterations (adjustable)

Example:

(Solve x^2+2x+3=0)
[RUN]
N=2 A(0)=1 A(1)=2 A(2)=3
T=1E-6 L=20
(-1,-1.414214)
(-1,1.414214)

(Solve x^3-2x^2-x+2=0)
[RUN]
N=3 A(0)=1 A(1)=-2 A(2)=-1 A(3)=2
T=1E-6 L=20
(1,0)
(2,0)
(-1,0)

(Solve x^3+3x^2+x+3=0)
[RUN]
N=3 A(0)=1 A(1)=3 A(2)=1 A(3)=3
T=1E-6 L=20
(0,1)
(0,-1)
(-3,0)

(Solve x^4+2999x^3-10003e3x^2-2399e7x+24e9=0)
[RUN]
N=4 A(0)=1 A(1)=2999 A(2)=-10003E3 A(3)=-2399E7 A(4)=24E9
T=1E-6 L=20
(1,0)
(-4000,0)
(-2000,0)
(3000,0)

(Solve 5x^6-45x^5+225x^4-425x^3+170x^2+370x-500=0)
[RUN]
N=6 A(0)=5 A(1)=-45 A(2)=225 A(3)=-425 A(4)=170 A(5)=370 A(6)=-500
T=1E-6 L=20
(1,1)
(3,4)
(-1,0)
(3,-4)
(1,-1)
(2,0)

Code:
' ABERTH POLYNOMIAL ROOT FINDER ' Estimate the complex roots of a polynomial with real coefficients ' For HP-71B + Math Pac by Robert van Engelen ' Optimized implementation: ' - Finds complex roots z of polynomial p such that |p(z)|<T given tolerance T ' - Prevents overflow/nan/inf by bounding the root stride modulus during polishing ' - Tested with a large sample of polynomials ' ' For details see: https://www.hpmuseum.org/forum/thread-18051.html ' See also: https://en.wikipedia.org/wiki/Durand–Kerner_method '           https://en.wikipedia.org/wiki/Aberth_method ' ' VARIABLES '  N         degree '  T         tolerance (1E-6) '  L         max iterations (20), updated to the iterations executed '  A(0..N)   N+1 real coefficients indexed from highest to lowest '  Z(0..N-1) N complex roots '  D         denominator '  E         estimated error '  H         root stride modulus '  P         p(z) and root stride p(z)/d '  Q         p'(z) '  R         root '  U         max root stride modulus '  I,J,K     loop iterators '  M         number of coefficients - 1 (M<=N-1) '  F         scratch ' ' DELAY 1 or smaller recommended 100 OPTION BASE 0 110 DESTROY ALL 120 INTEGER I,J,K,L,M,N @ REAL E,H,T,U 130 INPUT "N=";N @ IF N<2 THEN 130 140 REAL A(N) @ M=N-1 @ COMPLEX D,F,P,Q,R,Z(M) 150 MAT INPUT A 160 INPUT "T=","1E-6";T @ INPUT "L=","20";L ' remove leading zero coefficients (trailing A) 170 IF M>0 AND A(M+1)=0 THEN Z(M)=0 @ M=M-1 @ GOTO 170 180 IF M<0 OR A(0)=0 THEN DISP "ERROR" @ END ' generate initial roots 190 R=1 @ FOR I=0 TO M @ R=(.4,.9)*R @ Z(I)=R @ NEXT I ' max root stride modulus bounded to prevent overflow 200 U=MAXREAL^(1/(M+1)) @ IF U<1 THEN U=1 ' iterative root polishing until convergence or max iterations reached 210 FOR K=1 TO L @ E=0 220 FOR I=0 TO M ' compute p(z) and p'(z) using Horner's rule 230 P=A(0) @ Q=0 @ R=Z(I) @ D=0 240 FOR J=0 TO M @ Q=Q*R+P @ P=P*R+A(J+1) @ NEXT J ' optimized: if |p(z)|<tolerance then do not update root z 250 IF ABS(P)<T THEN 310 ' compute denominator 260 FOR J=0 TO M @ IF I=J THEN 280 270 F=R-Z(J) @ IF ABS(F)>=T THEN D=D+P/F 280 NEXT J ' compute p(z)/d if d is not too small to prevent overflow 290 D=Q-D @ IF ABS(D)>T THEN P=P/D ' polish root if its root stride modulus is within max bounds to prevent overflow 300 H=ABS(P) @ E=E+H @ IF H<U THEN Z(I)=R-P 310 NEXT I 320 DISP "K=";K;"E=";E 330 IF E<=T THEN J=L @ L=K @ K=J 340 NEXT K ' round results (optional) 350 U=1/T 360 FOR I=0 TO M @ Z(I)=(T*INT(U*REPT(Z(I))+.5),T*INT(U*IMPT(Z(I))+.5)) @ NEXT I ' estimate actual error 370 E=E+(E=0)*T ' display results 380 MAT DISP Z

"I count on old friends" -- HP 71B,Prime|Ti VOY200,Nspire CXII CAS|Casio fx-CG50...|Sharp PC-G850,E500,2500,1500,14xx,13xx,12xx...
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)