'***************************************************** '* Program to demonstrate the Lin's method * '* ------------------------------------------------- * '* Reference: BASIC Scientific Subroutines, Vol. II * '* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. * '* ------------------------------------------------- * '* Example: Find two complex roots of polynomial: * '* f(x) = x^5+5x^4+10x^3+10x^2+5x+1 * '* * '* SAMPLE RUN: * '* * '* Input order of polynomial: 5 * '* * '* Input the polynomial coefficients: * '* * '* A(0) = ? 1 * '* A(1) = ? 5 * '* A(2) = ? 10 * '* A(3) = ? 10 * '* A(4) = ? 5 * '* A(5) = ? 1 * '* * '* Convergence factor: 0 * '* Maximum number of iterations: 1000 * '* * '* The roots found are: * '* * '* X1 = -.93763285 * '* Y1 = 2.5730332E-02 * '* * '* X2 = -.93763285 * '* Y2 = -2.2730332E-02 * '* * '* The number of iterations was: 1000 * '* * '***************************************************** defint i-n defdbl a-h,o-z cls dim A(10),B(10),C(10) print input " Order of polynomial: ", m print print " Input the polynomial coefficients:" print for i=0 to m print " A(";i;") = "; : input A(i) next print input " Convergence factor: ", e input " Maximum number of iterations: ",n aa=3.1415926535# : bb=SQR(2#) gosub 1000 'Call LIN subroutine print print " The roots found are:" print print " X1 = "; x1 print " Y1 = "; y1 print print " X2 = "; x2 print " Y2 = "; y2 print print " The number of iterations was: "; k print end '************************************************* '* Polynomial complex roots subroutine * '* --------------------------------------------- * '* This routine uses the Lin's method * '* --------------------------------------------- * '* Reference: A Practical Guide to Computer * '* Method for Engineers by Shoup * '* --------------------------------------------- * '* INPUTS: * '* Polynomial coefficients : A(0) to A(m) * '* Order of polynomial : m * '* Initial guess : a and b * '* Convergence factor : e * '* Maximum number of iterations : n * '* OUTPUTS: * '* Two complex roots : x,y1 x2,y2 * '* Number of iterations : k * '************************************************* 1000 'Normalize the A(i) series for i=0 to m C(i)=A(i)/A(m) next 'Start iteration 'Set initial guess for the quadratic coefficients B(0)=0 : B(1)=0 1100 B(m-1)=C(m-1)-aa B(m-2)=C(m-2)-aa*B(m-1)-bb for j=3 to m B(m-j)=C(m-j)-aa*B(m+1-j)-bb*B(m+2-j) next 'Guard against divide by zero if B(2)<>0 then goto 1200 aa=aa+1e-7 : bb=bb-1e-7 1200 a1=(C(1)-bb*B(3))/B(2) b1=C(0)/B(2) k=k+1 'Test for the number of iterations if k>=n then goto 1300 'Test for convergence if ABS(aa-a1)+ABS(bb-b1) < e*e then goto 1300 aa=a1 : bb=b1 'Next iteration goto 1100 1300 aa=a1 : bb=b1 cc=aa*aa-4#*bb 'Is there an imaginary part? if cc>0 then goto 1400 y1=SQR(-cc) : y2=-y1 : x1=-aa : x2=x1 goto 1500 1400 y1=0 : y2=y1 : x1=-aa+SQR(cc) : x2=-aa-SQR(cc) 1500 x1=x1/2# x2=x2/2# : y1=y1/2# : y2=y2/2# return 'End of file Lin.bas 