alg121.txt

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

> restart;
> Digits := 20;
> # POISSON EQUATION FINITE-DIFFERENCE ALGORITHM 12.1
> #
> # To approximate the solution to the Poisson equation
> #            DEL(u) = F(x,y), a <= x <= b, c <= y <= d,
> # SUBJECT TO BOUNDARY CONDITIONS:
> #                 u(x,y) = G(x,y),
> #     if x = a or x = b for c <= y <= d,
> #     if y = c or y = d for a <= x <= b
> #
> # INPUT:   endpoints a, b, c, d; integers m, n; tolerance TOL;
> #          maximum number of iterations M
> #
> # OUTPUT:  approximations W(I,J) to u(X(I),Y(J)) for each
> #          I = 1,..., n-1 and J=1,..., m-1 or a message that the
> #          maximum number of iterations was exceeded.
> alg121 := proc() local F, G, OK, A, B, C, D, N, M, TOL, NN, M1, M2, N1, N2, H, K, I, X, J, Y, W, V, VV, L, Z, E, LL, FLAG, NAME, OUP;
> printf(`This is the Finite-Difference Method for Elliptic Equations.\n`);
> printf(`Input the functions F(X,Y) and G(X,Y) in terms of x and y 
> separated by a space\n`);
> printf(`For example:  x*exp(y) x*exp(y)\n`);
> F := scanf(`%a`)[1];
> G := scanf(`%a`)[1];
> F := unapply(F,x,y);
> G := unapply(G,x,y);
> OK :=FALSE;
> while OK = FALSE do
> printf(`Input endpoints of interval [A,B] on X-axis\n`);
> printf(`separated by a blank.\n`);
> A := scanf(`%f`)[1];
> B := scanf(`%f`)[1];
> printf(`Input endpoints of interval [C,D] on Y-axis\n`);
> printf(`separated by a blank.\n`);
> C := scanf(`%f`)[1];
> D := scanf(`%f`)[1];
> if A >= B or C >= D then
> printf(`Left endpoint must be less than right endpoint.\n`);
> else
> OK := TRUE;
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input number of intervals n on the X-axis and m\n`);
> printf(`on the Y-axis separated by a blank\n`);
> printf(`Note that both n and m should be larger than 2.\n`);
> N := scanf(`%d`)[1];
> M := scanf(`%d`)[1];
> if M <= 2 or N <= 2
> then printf(`Numbers must exceed 2.\n`);
> else
> OK := TRUE;
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the Tolerance.\n`);
> TOL := scanf(`%f`)[1];
> if TOL <= 0 then
> printf(`Tolerance must be positive.\n`);
> else
> OK := TRUE;
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the maximum number of iterations.\n`);
> NN := scanf(`%d`)[1];
> if NN <= 0 then
> printf (`Number must be a positive integer.\n`);
> else
> OK := TRUE;
> fi;
> od;
> if OK = TRUE then
> M1 := M-1;
> M2 := M-2;
> N1 := N-1;
> N2 := N-2;
> # Step 1
> H := (B-A)/N;
> K := (D-C)/M;
> # Steps 2 and 3 construct mesh points
> # Step 2
> for I from 0 to N do
> X[I] := A+I*H;
> od;
> # Step 3
> for J from 0 to M do
> Y[J] := C+J*K;
> od;
> # Step 4
> for I from 1 to N1 do
> W[I,0] := G(X[I],Y[0]);
> W[I,M] := G(X[I],Y[M]);
> od;
> for J from 0 to M do
> W[0,J] := G(X[0],Y[J]);
> W[N,J] := G(X[N],Y[J]);
> od;
> for I from 1 to N1 do
> for J from 1 to M1 do
> W[I,J] := 0;
> od;
> od;
> # Step 5
> # Use V in place of lambda, VV in place of mu
> V := H*H/(K*K);
> VV := 2*(1+V);
> L := 1;
> OK := FALSE;
> # Z is a new value of W(I,J) to be used in computing the norm of the
> # error E
> # Step 6
> while L <= NN and OK = FALSE do
> # Steps 7 - 20 perform Gauss-Seidel iterations
> # Step 7
> Z := (-H*H*F(X[1],Y[M1])+G(A,Y[M1])+V*G(X[1],D)+V*W[1,M2]+W[2,M1])/VV;
> E := abs( W[1,M1]-Z);
> W[1,M1] := Z;
> # Step 8
> for I from 2 to N2 do
> Z := (-H*H*F(X[I],Y[M1])+V*G(X[I],D)+W[I-1,M1]+W[I+1,M1]+V*W[I,M2])/VV;
> if abs(W[I,M1]-Z) > E then
> E := abs( W[I,M1] - Z );
> fi;
> W[I,M1] := Z;
> od;
> # Step 9
> Z := 
> (-H*H*F(X[N1],Y[M1])+G(B,Y[M1])+V*G(X[N1],D)+W[N2,M1]+V*W[N1,M2])/VV;
> if abs( W[N1,M1]-Z) > E then
> E := abs( W[N1,M1]-Z);
> fi;
> W[N1,M1] := Z;
> # Step 10
> for LL from 2 to M2 do
> J := M2-LL+2;
> # Step 11
> Z := (-H*H*F(X[1],Y[J])+G(A,Y[J])+V*W[1,J+1]+V*W[1,J-1]+W[2,J])/VV;
> if abs(W[1,J]-Z) > E then
> E := abs(W[1,J]-Z);
> fi;
> W[1,J] := Z;
> # Step 12
> for I from 2 to N2 do
> Z := (-H*H*F(X[I],Y[J])+W[I-1,J]+V*W[I,J+1]+V*W[I,J-1]+W[I+1,J])/VV;
> if abs(W[I,J]-Z) > E then
> E := abs(W[I,J]-Z);
> fi;
> W[I,J] := Z;
> od;
> # Step 13
> Z := (-H*H*F(X[N1],Y[J])+G(B,Y[J])+W[N2,J]+V*W[N1,J+1]+V*W[N1,J-1])/VV;
> if abs(W[N1,J]-Z) > E then
> E := abs(W[N1,J]-Z);
> fi;
> W[N1,J] := Z;
> od;
> # Step 14
> Z := (-H*H*F(X[1],Y[1])+V*G(X[1], C)+G(A,Y[1])+V*W[1,2]+W[2,1])/VV;
> if abs(W[1,1]-Z) > E then
> E := abs(W[1,1]-Z);
> fi;
> W[1,1] := Z;
> # Step 15
> for I from 2 to N2 do
> Z := (-H*H*F(X[I],Y[1])+V*G(X[I],C)+W[I+1,1]+W[I-1,1]+V*W[I,2])/VV;
> if abs(W[I,1]-Z) > E then
> E := abs(W[I,1]-Z);
> fi;
> W[I,1] := Z;
> od;
> # Step 16
> Z := (-H*H*F(X[N1],Y[1])+V*G(X[N1],C)+G(B,Y[1])+W[N2,1]+V*W[N1,2])/VV;
> if abs(W[N1,1]-Z) > E then
> E := abs(W[N1,1]-Z);
> fi;
> W[N1,1] := Z;
> # Step 17
> if E <= TOL then
> # Step 18
> printf(`Choice of output method:\n`);
> printf(`1. Output to screen\n`);
> printf(`2. Output to text file\n`);
> printf(`Please 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, `POISSON EQUATION FINITE-DIFFERENCE METHOD\n\n`);
> fprintf(OUP,                                             `  I  J    X(I)        Y(J)         W(I,J)\n\n`);
> for I from 1 to N1 do
> for J from 1 to M1 do
> fprintf(OUP, `%3d %2d %11.8f %11.8f %13.8f\n`, I, J, X[I], Y[J], W[I,J]);
> od;
> od;
> fprintf(OUP, `Convergence occurred on iteration number: %d\n`, L); 
> # Step 19
> OK := TRUE;
> else
> # Step 20
> L := L+1;
> fi;
> od;
> # Step 21
> if L > NN then
> printf(`Method fails after iteration number %d\n`, NN)
> else
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> fi;
> RETURN(0);
> end;
> alg121();