'******************************************************** '* Equations of degree 2 to 4 * '* ---------------------------------------------------- * '* This program calculates the real or complex roots of * '* algebraic equations of degree 2, 3 and 4. * '* ---------------------------------------------------- * '* Reference: Mathematiques et statitiques - Programmes * '* en BASIC. Editions du P.S.I., 1981 * '* [BIBLI 13]. * '* ---------------------------------------------------- * '* SAMPLE RUN: * '* * '* ROOTS OF EQUATIONS OF DEGREE 2, 3 OU 4 * '* Example: X^4+4X^2-3X+3 = 0 * '* * '* INSTRUCTIONS: * '* * '* 1. Input degree of equation * '* 2. Input coefficients beginning with highest exponent* '* 3. Solutions are given under the form: * '* (real part) + I (imaginary part) * '* * '* Degree of equation (2 to 4): 4 * '* * '* Coefficient of X^4 = 1 * '* Coefficient of X^3 = 0 * '* Coefficient of X^2 = 4 * '* Coefficient of X^1 = -3 * '* Coefficient of X^0 = 3 * '* * '* Root 1 = 0.449702650 + I 0.731070280 * '* Root 2 = 0.449702650 - I 0.731070280 * '* Root 1 = -0.449702650 + I 1.967231870 * '* Root 2 = -0.449702650 - I 1.967231870 * '* * '******************************************************** CLS DIM a(4, 4), r(4) F\$="##.#########" PRINT " ROOTS OF EQUATIONS OF DEGREE 2, 3 OU 4" PRINT " Example: X^4+4X^2-3X+3 = 0" PRINT PRINT " INSTRUCTIONS:" PRINT PRINT " 1. Input degree of equation" PRINT " 2. Input coefficients beginning by highest exponent" PRINT " 3. Solutions are given under the form:" PRINT " (real part) + I (imaginary part)" PRINT INPUT " Degree of equation (2 to 4): ", n IF n < 2 THEN n = 2 END IF IF n > 4 THEN n = 4 END IF PRINT FOR i = n TO 0 STEP -1 PRINT " Coefficient of X^"; i; : INPUT " = ", a(n, i) NEXT i PRINT : PRINT ' calling root search subroutine GOSUB 1000 ' printing results ON n - 1 GOTO 400, 500, 600 ' case n=2 400 IF im = 0 THEN PRINT " Root 1 = "; : print using F\$; r(1) PRINT " Root 2 = "; : print using F\$; r(2) ELSE PRINT " Root 1 = "; : print using F\$; r(1); PRINT " +I "; : print using F\$; im PRINT " Root 2 = "; : print using F\$; r(2); PRINT " -I "; : print using F\$; im END IF END ' case n=3 500 PRINT " Root 1 = "; : print using F\$; r(3) IF im = 0 THEN PRINT " Root 2 = "; : print using F\$; r(1) PRINT " Root 3 = "; : print using F\$; r(2) ELSE PRINT " Root 2 = "; : print using F\$; r(1); PRINT " +I "; : print using F\$; im PRINT " Root 3 = "; : print using F\$; r(2); PRINT " -I "; : print using F\$; im END IF END ' cas n=4 600 IF im = 0 THEN PRINT " Root 1 = "; : print using F\$; r(1) PRINT " Root 2 = "; : print using F\$; r(2) ELSE PRINT " Root 1 = "; : print using F\$; r(1); PRINT " +I "; : print using F\$; im PRINT " Root 2 = "; : print using F\$; r(2); PRINT " -I "; : print using F\$; im END iF IF ii = 0 THEN PRINT " Root 3 = "; : print using F\$; r(3) PRINT " Root 4 = "; : print using F\$; r(4) ELSE PRINT " Root 3 = "; : print using F\$; r(3); PRINT " +I "; : print using F\$; ii PRINT " Root 4 = "; : print using F\$; r(4); PRINT " -I "; : print using F\$; ii END IF END 'Roots search subroutine 1000 'Normalization of coefficients FOR i = 0 TO n - 1 a(n, i) = a(n, i) / a(n, n) NEXT i ' Branch according to n ON n - 1 GOTO 1400, 1500, 1600 ' cas n=2 1400 GOSUB 3000 RETURN ' cas n=3 1500 GOSUB 2000 RETURN ' cas n=4 1600 a = a(4, 3) b = a(4, 2) c = a(4, 1) d = a(4, 0) q = b - (3 * a * a / 8) r = c - (a * b / 2) + (a * a * a / 8) s = d - (a * c / 4) + (a * a * b / 16) - (3 * a * a * a * a / 256) ' Definition of coefficients of cubic equation a(3, 2) = q / 2 a(3, 1) = (q * q - 4 * s) / 16 a(3, 0) = -(r * r / 64) ' Calculates a real root of this equation IF (r <> 0) OR (a(3, 1) >= 0) THEN GOTO 1650 ' Case when this equation is reduced to 2nd degree a(2, 1) = a(3, 2): a(2, 0) = a(3, 1) GOSUB 3000 ' One takes positive root r(3) = r(1) GOTO 1700 ' Calling subroutine 2000 with sw=-1 to calculate only one root 1650 sw = -1 GOSUB 2000 ' real root of above equation 1700 k = SQR(r(3)) ' Calculates L and M when k=0 IF k = 0 THEN r = SQR(q * q - 4 * s): GOTO 1750 END IF q = q + (4 * r(3)) r = r / (2 * k) 1750 l = (q - r) / 2: m = (q + r) / 2 ' Solving two 2nd degree equations a(2, 1) = 2 * k: a(2, 0) = l ' 1st equation GOSUB 3000 ' Solutions transferred in r(3), r(4), ii r(3) = r(1) - (a(4, 3) / 4): r(4) = r(2) - (a(4, 3) / 4): ii = im a(2, 1) = -2 * k: a(2, 0) = m ' 2nd equation GOSUB 3000 r(2) = r(2) - (a(4, 3) / 4): r(1) = r(1) - (a(4, 3) / 4) RETURN '******************************************** '* This subroutine calculates the roots of * '* X^3 + A(3,2)*X^2 + A(3,1)*X + A(3,0) = 0 * '******************************************** ' Case when a root equals zero 2000 IF a(3, 0) = 0 THEN xc = 0: GOTO 2500 END IF ' Searching highest coefficient in absolute value am = ABS(a(3, 0)) FOR i = 1 TO 3 tt = ABS(a(3, i)) IF (am < tt) THEN am = tt NEXT i ' Define interval where a real root exists ' according to sign of A(3,0) IF a(3, 0) > 0 THEN a = -am - 1 b = 0 GOTO 2100 END IF a = 0: b = am + 1 2100 ' Searching for xc = real root in interval (a,b) ' (by Bisection method) DEF fnf (x) = x * x * x + a(3, 2) * x * x + a(3, 1) * x + a(3, 0) ' Define segment (xa,xb) including the root xa = a: xb = b: y1 = fnf(a): y2 = fnf(b) ' Desired precision er = .000001 IF y1 = 0 THEN xc = a: GOTO 2500 END IF IF y2 = 0 THEN xc = b: GOTO 2500 END IF ' xc = middle of segment 2200 xc = (xa + xb) / 2: te = fnf(xc) IF te = 0 THEN GOTO 2500 IF (xb - xa) < er THEN GOTO 2500 IF (fnf(xa) * te) > 0 THEN xa = xc: GOTO 2200 END IF xb = xc: GOTO 2200 ' r(3) is a real root 2500 r(3) = xc IF sw = -1 THEN RETURN ' Calculates the roots of remaining quadratic equation ' Define the equation coefficients a(2, 1) = a(3, 2) + xc a(2, 0) = a(3, 1) + a(3, 2) * xc + xc * xc GOSUB 3000 RETURN 'End subroutine 2000 '******************************************* '* This Subroutine calculates the roots of * '* X^2 + A(2,1)*X + A(2,0) = 0 * '* W = Determinant of equation * '******************************************* 3000 w = a(2, 1) * a(2, 1) - 4 * a(2, 0) q1 = -a(2, 1) / 2 ON (SGN(w) + 2) GOTO 3200, 3300, 3400 ' Cas W < 0 3200 q2 = SQR(-w) / 2 im = q2 r(1) = q1: r(2) = q1 RETURN ' Cas W = 0 3300 r(1) = q1: r(2) = q1: im = 0 ' Cas W > 0 3400 q2 = SQR(w) / 2: im = 0 r(1) = q1 + q2: r(2) = q1 - q2 RETURN 'End of file root4.bas 