> 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.
> print(`This is the Finite-Difference Method for Elliptic Equations.`):
> print(`Input the functions F(X,Y) and G(X,Y) in terms of x and y 
> separated by a space`):
> print(`For example:  x*exp(y) x*exp(y)`):
> F := scanf(`%a`)[1]:
> G := scanf(`%a`)[1]:print(`F(x,y) = `):print(F):print(`G(x,y) = `):print(G):
> F := unapply(F,x,y):
> G := unapply(G,x,y):
> OK :=FALSE:
> while OK = FALSE do
> print(`Input endpoints of interval [A,B] on X-axis`):
> print(`separated by a blank.`):
> A := scanf(`%f`)[1]:
> B := scanf(`%f`)[1]:print(`A = `):print(A):print(`B = `):print(B):
> print(`Input endpoints of interval [C,D] on Y-axis`):
> print(`separated by a blank.`):
> C := scanf(`%f`)[1]:
> D2 := scanf(`%f`)[1]:print(`C = `):print(C):print(`D = `):print(D2):
> if A >= B or C >= D2 then
> print(`Left endpoint must be less than right endpoint.\n`):
> else
> OK := TRUE:
> fi:
> od:
> OK := FALSE:
> while OK = FALSE do
> print(`Input number of intervals n on the X-axis and m`):
> print(`on the Y-axis separated by a blank`):
> print(`Note that both n and m should be larger than 2.`):
> N := scanf(`%d`)[1]:print(`Number on x-axis N = `):print(N):
> M := scanf(`%d`)[1]:print(`Number on y-axis M = `):print(M):
> if M <= 2 or N <= 2
> then print(`Numbers must exceed 2.`):
> else
> OK := TRUE:
> fi:
> od:
> OK := FALSE:
> while OK = FALSE do
> print(`Input the Tolerance.`):
> TOL := scanf(`%f`)[1]: print(`Tolerance = `):print(TOL):
> if TOL <= 0 then
> print(`Tolerance must be positive.`):
> else
> OK := TRUE:
> fi:
> od:
> OK := FALSE:
> while OK = FALSE do
> print(`Input the maximum number of iterations.`):
> NN := scanf(`%d`)[1]:print(`Maximum number of iterations = `):print(NN):
> if NN <= 0 then
> print (`Number must be a positive integer.`):
> 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 := (D2-C)/M:
> # Steps 2 and 3 construct mesh points
> # Step 2
> for I2 from 0 to N do
> X[I2] := A+I2*H:
> od:
> # Step 3
> for J from 0 to M do
> Y[J] := C+J*K:
> od:
> # Step 4
> for I2 from 1 to N1 do
> W[I2,0] := G(X[I2],Y[0]):
> W[I2,M] := G(X[I2],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 I2 from 1 to N1 do
> for J from 1 to M1 do
> W[I2,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],D2)+V*W[1,M2]+W[2,M1])/VV:
> E := abs( W[1,M1]-Z):
> W[1,M1] := Z:
> # Step 8
> for I2 from 2 to N2 do
> Z := (-H*H*F(X[I2],Y[M1])+V*G(X[I2],D2)+W[I2-1,M1]+W[I2+1,M1]+V*W[I2,M2])/VV:
> if abs(W[I2,M1]-Z) > E then
> E := abs( W[I2,M1] - Z ):
> fi:
> W[I2,M1] := Z:
> od:
> # Step 9
> Z := 
> (-H*H*F(X[N1],Y[M1])+G(B,Y[M1])+V*G(X[N1],D2)+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 I2 from 2 to N2 do
> Z := (-H*H*F(X[I2],Y[J])+W[I2-1,J]+V*W[I2,J+1]+V*W[I2,J-1]+W[I2+1,J])/VV:
> if abs(W[I2,J]-Z) > E then
> E := abs(W[I2,J]-Z):
> fi:
> W[I2,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 I2 from 2 to N2 do
> Z := (-H*H*F(X[I2],Y[1])+V*G(X[I2],C)+W[I2+1,1]+W[I2-1,1]+V*W[I2,2])/VV:
> if abs(W[I2,1]-Z) > E then
> E := abs(W[I2,1]-Z):
> fi:
> W[I2,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
> print(`Choice of output method:\n`):
> print(`1. Output to screen\n`):
> print(`2. Output to text file\n`):
> print(`Please enter 1 or 2.\n`):
> FLAG := scanf(`%d`)[1]: print(`Input is `):print(FLAG):
> if FLAG = 2 then
> print(`Input the file name in the form - drive:\\name.ext`):
> print(`for example:  A:\\OUTPUT.DTA`):
> NAME := scanf(`%s`)[1]: print(`Output file is `):print(NAME):
> 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 I2 from 1 to N1 do
> for J from 1 to M1 do
> fprintf(OUP, `%3d %2d %11.8f %11.8f %13.8f\n`, I2, J, X[I2], Y[J], W[I2,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
> print(`Method fails after iteration number `):print(NN):
> else
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi:
> fi:
> fi: