restart:# MULLER'S ALGORITHM 2.8## To find a solution to f(x) = 0 given three approximations x0, x1# and x2:## INPUT: x0,x1,x2: tolerance TOL: maximum number of iterations NO.## OUTPUT: approximate solution p or message of failure.## This implementation allows for a switch to complex arithmetic.# The coefficients are stored in the vector A, so the dimension# of A may have to be changed.print(`This is Mullers Method.`):print(`Input the Polynomial P(x) in terms of x `):print(`For example: x^2-2*x+2 `):P := scanf(`%a`)[1]:print(`P(x) = `):print(P):P := unapply(P,x):OK := TRUE:if degree(P(x)) = 1 thenpp := -tcoeff(P(x))/coeff(P(x),x):print(`Polynomial is linear: root is `):print(pp):OK := FALSE:fi:if degree(P(x)) < 1 thenprint(`Polynomial is a constant function. `):OK := FALSE:fi:if OK = TRUE thenOK := FALSE:while OK = FALSE doprint(`Input tolerance `):TOL := scanf(`%a`)[1]:print(`Tolerance = `):print(TOL):if TOL <= 0 thenprint(`Tolerance must be positive `):else OK := TRUE:fi:od:OK := FALSE:while OK = FALSE doprint(`Input maximum number of iterations - no decimal point `):M := scanf(`%d`)[1]:print(`Maximum number of iterations = `):print(M):if M <= 0 thenprint(`Must be positive integer `):else OK := TRUE:fi:od:OK := FALSE:while OK = FALSE doprint(`Input the first of three starting values `):X[0] := scanf(`%f`)[1]:print(`First starting value = `):print(X[0]):print(`Input the second of three starting values `):X[1] := scanf(`%f`)[1]:print(`Second starting value = `):print(X[1]):print(`Input the third starting value `):X[2] := scanf(`%f`)[1]:print(`Third starting value = `):print(X[2]):if ((abs(X[0]-X[1]) < 1.0e-20) or (abs(X[0]-X[2]) < 1.0e-20) or (abs(X[2]-X[1]) < 1.0e-20)) then print(`The starting values must be distinct. `):elseOK := TRUE:fi:od:fi:F[0] := P(X[0]):F[1] := P(X[1]):F[2] := P(X[2]):J := -1:if F[0] = 0 then J := 0: fi:if F[1] = 0 then J := 1: fi:if F[2] = 0 then J := 2: fi:if J > -1 thenprint(`Solution is the initial approximation `):print(X[J]):OK := FALSE:fi:if OK = TRUE thenprint(`Select output destination `):print(`1. Screen `):print(`2. Text file `):print(`Enter 1 or 2 `):FLAG := scanf(`%d`)[1]:print(`Input = `):print(FLAG):if FLAG = 2 thenprint(`Input the file name in the form - drive:\134\134name.ext `):print(`For example: A:\134\134OUTPUT.DTA `):NAME := scanf(`%s`)[1]:print(`Input = `):print(NAME):OUP := fopen(NAME,WRITE,TEXT):elseOUP := default:fi:fprintf(OUP, `MULLER'S METHOD \134n`):fprintf(OUP, `The output is i, approximation x(i), f(x(i)) \134n `):fprintf(OUP,`Three columns means the results are real numbers,\134n `):fprintf(OUP,`Five columns means the results are complex\134n`):fprintf(OUP,`numbers with real and imaginary parts of x(i)\134n`):fprintf(OUP,`followed by real and imaginary parts of f(x(i)).\134n\134n`):I1:=0:fprintf(OUP, `%d %a %a\134n`, I1, X[I1], F[I1]):I1:=1:fprintf(OUP, `%d %a %a\134n`, I1, X[I1], F[I1]):I1:=2;fprintf(OUP, `%d %a %a\134n`, I1, X[I1], F[I1]):# Step 1H[0] := X[1]-X[0]:H[1] := X[2]-X[1]:DEL1[0] := (F[1]-F[0])/H[0]:DEL1[1] := (F[2]-F[1])/H[1]:DEL := (DEL1[1]-DEL1[0])/(H[1]+H[0]):I1 := 3:# Step 2while I1 <= M and OK = TRUE do# Step 3B := DEL1[1]+H[1]*DEL:D1 := B*B-4*F[2]*DEL:# Test to see if straight lineif abs(DEL) <= 1.0e-20 then# straight line - test if horizontal lineif abs(DEL1[1]) <= 1.0e-20 thenprint(`Horizontal Line `):OK := FALSE:else# Straight line, but not horizontalX[3] := (F[2]-DEL1[1]*X[2])/DEL1[1]:H[2] := X[3]-X[2]:fi:else# Not a straight lineD1 := sqrt(D1):# Step 4E := B+D1:if abs(B-D1) > abs(E) thenE := B-D1:fi:# Step 5H[2] := -2*F[2]/E:X[3] := X[2]+H[2]:fi:if OK = TRUE thenF[3] := P(X[3]):fprintf(OUP, `%d %a %a\134n`, I1, X[3], F[3]):fi:# Step 6if ( abs(H[2]) < TOL or F[3] = 0 ) thenfprintf(OUP, `\134nMethod is successfull. The approximation\134n`):fprintf(OUP, `%a is within %.10e\134n`,X[3], TOL):fprintf(OUP, `in %d iterations.\134n`, I1):OK := FALSE:else# Step 7for J from 1 to 2 doH[J-1] := H[J]:X[J-1] := X[J]:F[J-1] := F[J]:od:X[2] := X[3]:F[2] := F[3]:DEL1[0] := DEL1[1]:DEL1[1] := (F[2]-F[1])/H[1]:DEL := (DEL1[1]-DEL1[0])/(H[1]+H[0]):if abs(DEL) <= 1.0e-20 thenif abs(DEL1[1]) <= 1.0e-20 thenprint(`Horizontal Line `):OK := FALSE:fi:fi:fi:I1 := I1+1:od:# Step 8if I1 > M and OK = TRUE thenfprintf (OUP, `Final approximation not within tolerance %.10e\134n`,TOL):fi:if OUP <> default thenfclose(OUP):print(`Output file `, NAME,` created successfully`):fi:fi:JSFH