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.\134n`):

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:

# 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:\134n`):

print(`1. Output to screen\134n`):

print(`2. Output to text file\134n`):

print(`Please enter 1 or 2.\134n`):

FLAG := scanf(`%d`)[1]: print(`Input is `):print(FLAG):

if FLAG = 2 then

print(`Input the file name in the form - drive:\134\134name.ext`):

print(`for example: A:\134\134OUTPUT.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\134n\134n`):

fprintf(OUP, ` I J X(I) Y(J) W(I,J)\134n\134n`):

for I2 from 1 to N1 do

for J from 1 to M1 do

fprintf(OUP, `%3d %2d %11.8f %11.8f %13.8f\134n`, I2, J, X[I2], Y[J], W[I2,J]):

od:

fprintf(OUP, `Convergence occurred on iteration number: %d\134n`, 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:

SWduVGhpc35pc350aGV+RmluaXRlLURpZmZlcmVuY2V+TWV0aG9kfmZvcn5FbGxpcHRpY35FcXVhdGlvbnMuRzYi

SWpvSW5wdXR+dGhlfmZ1bmN0aW9uc35GKFgsWSl+YW5kfkcoWCxZKX5pbn50ZXJtc35vZn54fmFuZH55fnwrc2VwYXJhdGVkfmJ5fmF+c3BhY2VHNiI=

SUBGb3J+ZXhhbXBsZTp+fngqZXhwKHkpfngqZXhwKHkpRzYi

SSpGKHgseSl+PX5HNiI=

KiZJInhHNiIiIiItSSRleHBHNiQlKnByb3RlY3RlZEdJKF9zeXNsaWJHRiQ2I0kieUdGJEYl

SSpHKHgseSl+PX5HNiI=

KiZJInhHNiIiIiItSSRleHBHNiQlKnByb3RlY3RlZEdJKF9zeXNsaWJHRiQ2I0kieUdGJEYl

SUxJbnB1dH5lbmRwb2ludHN+b2Z+aW50ZXJ2YWx+W0EsQl1+b25+WC1heGlzRzYi

STZzZXBhcmF0ZWR+Ynl+YX5ibGFuay5HNiI=

SSVBfj1+RzYi

JCIiIUYj

SSVCfj1+RzYi

JCIiIyIiIQ==

SUxJbnB1dH5lbmRwb2ludHN+b2Z+aW50ZXJ2YWx+W0MsRF1+b25+WS1heGlzRzYi

STZzZXBhcmF0ZWR+Ynl+YX5ibGFuay5HNiI=

SSVDfj1+RzYi

JCIiIUYj

SSVEfj1+RzYi

JCIiIiIiIQ==

SVBJbnB1dH5udW1iZXJ+b2Z+aW50ZXJ2YWxzfm5+b25+dGhlflgtYXhpc35hbmR+bUc2Ig==

SUNvbn50aGV+WS1heGlzfnNlcGFyYXRlZH5ieX5hfmJsYW5rRzYi

SVBOb3RlfnRoYXR+Ym90aH5ufmFuZH5tfnNob3VsZH5iZX5sYXJnZXJ+dGhhbn4yLkc2Ig==

STZOdW1iZXJ+b25+eC1heGlzfk5+PX5HNiI=

IiIn

STZOdW1iZXJ+b25+eS1heGlzfk1+PX5HNiI=

IiIm

STVJbnB1dH50aGV+VG9sZXJhbmNlLkc2Ig==

SS1Ub2xlcmFuY2V+PX5HNiI=

JCIiIiEjNQ==

SUhJbnB1dH50aGV+bWF4aW11bX5udW1iZXJ+b2Z+aXRlcmF0aW9ucy5HNiI=

SUBNYXhpbXVtfm51bWJlcn5vZn5pdGVyYXRpb25zfj1+RzYi

IiQrIg==

STpDaG9pY2V+b2Z+b3V0cHV0fm1ldGhvZDp8K0c2Ig==

STUxLn5PdXRwdXR+dG9+c2NyZWVufCtHNiI=

STgyLn5PdXRwdXR+dG9+dGV4dH5maWxlfCtHNiI=

STZQbGVhc2V+ZW50ZXJ+MX5vcn4yLnwrRzYi

SSpJbnB1dH5pc35HNiI=

IiIi

POISSON EQUATION FINITE-DIFFERENCE METHOD I J X(I) Y(J) W(I,J) 1 1 0.33333333 0.20000000 0.40726462 1 2 0.33333333 0.40000000 0.49748325 1 3 0.33333333 0.60000000 0.60759608 1 4 0.33333333 0.80000000 0.74200707 2 1 0.66666667 0.20000000 0.81452375 2 2 0.66666667 0.40000000 0.99495758 2 3 0.66666667 0.60000000 1.21518318 2 4 0.66666667 0.80000000 1.48400854 3 1 1.00000000 0.20000000 1.22176621 3 2 1.00000000 0.40000000 1.49240473 3 3 1.00000000 0.60000000 1.82274271 3 4 1.00000000 0.80000000 2.22599270 4 1 1.33333333 0.20000000 1.62896369 4 2 1.33333333 0.40000000 1.98977830 4 3 1.33333333 0.60000000 2.43022656 4 4 1.33333333 0.80000000 2.96792825 5 1 1.66666667 0.20000000 2.03604232 5 2 1.66666667 0.40000000 2.48695839 5 3 1.66666667 0.60000000 3.03750550 5 4 1.66666667 0.80000000 3.70972407 Convergence occurred on iteration number: 61

JSFH