alg028.txt

Numerical Anaysis 8th Edition Burden and Faires (Maple Source)

> 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.
> alg028 := proc() local P, OK, TOL, M, X, FLAG, NAME, OUP, F, H, DEL1, DEL, I, B, D, E, J;
> printf(`This is Mullers Method.\n`);
> printf(`Input the Polynomial P(x) in terms of x\n`);
> printf(`For example: x^2-2*x+2\n`);
> P := scanf(`%a`)[1];
> P := unapply(P,x);
> OK := TRUE;
> if degree(P(x)) = 1 then
> pp := -tcoeff(P(x))/coeff(P(x),x);
> printf(`Polynomial is linear: root is %11.8f\n`, pp);
> OK := FALSE;
> fi;
> if degree(P(x)) < 1 then
> printf(`Polynomial is a constant function.\n`);
> OK := FALSE;
> fi;
> if OK = TRUE then
> OK := FALSE;
> while OK = FALSE do
> printf(`Input tolerance\n`);
> TOL := scanf(`%a`)[1];
> if TOL <= 0 then
> printf(`Tolerance must be positive\n`);
> else 
> OK := TRUE;
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input maximum number of iterations - no decimal point\n`);
> M := scanf(`%d`)[1];
> if M <= 0 then
> printf(`Must be positive integer\n`);
> else 
> OK := TRUE;
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the first of three starting values\n`);
> X[0] := scanf(`%f`)[1];
> printf(`Input the second of three starting values\n`);
> X[1] := scanf(`%f`)[1];
> printf(`Input the third starting value\n`);
> X[2] := scanf(`%f`)[1];
> 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
>    printf(`The starting values must be distinct.\n`);
> else
> OK := 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 then
> printf(`Solution is the initial approximation %.10e \n`,X[J]);
> OK := FALSE;
> fi;
> if OK = TRUE then
> printf(`Select output destination\n`);
> printf(`1. Screen\n`);
> printf(`2. Text file\n`);
> printf(`Enter 1 or 2\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 2 then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`For example:   A:\\OUTPUT.DTA\n`);
> NAME := scanf(`%s`)[1];
> OUP := fopen(NAME,WRITE,TEXT);
> else
> OUP := default;
> fi;
> fprintf(OUP, `MULLER'S METHOD\n`);
> fprintf(OUP, `The output is i, approximation x(i), f(x(i))\n\n`);
> fprintf(OUP,`Three columns means the results are real numbers,\n`);
> fprintf(OUP,`Five columns means the results are complex\n`);
> fprintf(OUP,`numbers with real and imaginary parts of x(i)\n`);
> fprintf(OUP,`followed by real and imaginary parts of f(x(i)).\n\n`);
> # Step 1
> H[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]);
> I := 3;
> # Step 2
> while I <= M and OK = TRUE do
> # Step 3
> B := DEL1[1]+H[1]*DEL;
> D := B*B-4*F[2]*DEL;
> # Test to see if straight line
> if abs(DEL) <= 1.0e-20 then
> # straight line - test if horizontal line
> if abs(DEL1[1]) <= 1.0e-20 then
> printf(`Horizontal Line\n`);
> OK := FALSE;
> else
> # Straight line, but not horizontal
> X[3] := (F[2]-DEL1[1]*X[2])/DEL1[1];
> H[2] := X[3]-X[2];
> fi;
> else
> # Not a straight line
> D := sqrt(D);
> # Step 4
> E := B+D;
> if abs(B-D) > abs(E) then
> E := B-D;
> fi;
> # Step 5
> H[2] := -2*F[2]/E;
> X[3] := X[2]+H[2];
> fi;
> if OK = TRUE then
> F[3] := P(X[3]);
> fprintf(OUP, `%d   %a   %a\n`, I, X[3], F[3]);
> fi;
> # Step 6
> if ( abs(H[2]) < TOL or F[3] = 0 ) then
> fprintf(OUP, `\nMethod is successfull.  The approximation\n`);
> fprintf(OUP, `%a   is within %.10e\n`,X[3], TOL);
> fprintf(OUP, `in %d iterations.\n`, I);
> OK := FALSE;
> else
> # Step 7
> for J from 1 to 2 do
> H[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 then
> if abs(DEL1[1]) <= 1.0e-20 then
> printf(`Horizontal Line\n`);
> OK := FALSE;
> fi;
> fi;
> fi;
> I := I+1;
> od;
> # Step 8
> if I > M and OK = TRUE then
> fprintf (OUP, `Final approximation not within tolerance %.10e\n`,TOL);
> fi;
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg028();