alg041.txt

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

> restart;
> # SIMPSON'S COMPOSITE ALGORITHM 4.1
> #
> # To approximate I = integral ( ( f(x) dx ) ) from a to b:
> #
> # INPUT:   endpoints a, b; even positive integer n.
> #
> # OUTPUT:  approximation XI to I.
> alg041 := proc() local F, OK, A, B, N, H, XI0, XI1, XI2, NN, I, X, XI;
> printf(`This is Simpsons Method.\n\n`);
> printf(`Input the function F(x) in terms of x\n`);
> printf(`For example: cos(x)\n`);
> F := scanf(`%a`)[1];
> F := unapply(F,x);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input lower limit of integration and `);
> printf(`upper limit of integration\n`);
> printf(`separated by a blank\n`);
> A := scanf(`%f`)[1];
> B := scanf(`%f`)[1];
> if A > B then
> printf(`Lower limit must be less than upper limit\n`);
> else
> OK := TRUE;
> fi;
> od; 
> OK := FALSE;
> while OK = FALSE do
> printf(`Input an even positive integer N.\n`);
> N := scanf(`%d`)[1];
> if N > 0 and N mod 2 = 0 then
> OK := TRUE;
> else
> printf(`Input must be even and positive\n`);
> fi;
> od;
> if OK = TRUE then
> # Step 1
> H := (B-A)/N;
> # Step 2
> XI0 := F(A) + F(B);
> # Summation of f(x(2*I-1))
> XI1 := 0.0;
> # Summation of f(x(2*I))
> XI2 := 0.0;
> # Step 3
> NN := N - 1;
> for I from 1 to NN do
> # Step 4
> X := A + I * H;
> # Step 5
> if I mod 2 = 0 then 
> XI2 := XI2 + F(X);
> else
> XI1 := XI1 + F(X);      
> fi;
> od;
> # Step 6
> XI := (XI0 + 2.0 * XI2 + 4.0 * XI1) * H / 3.0;
> # Step 7
> printf(`\nThe integral of F from %12.8f to %12.8f is\n`, A, B);
> printf(`%12.8f\n`, XI);
> fi;
> RETURN(0);
> end;
> alg041();