restart:

# HEAT EQUATION BACKWARD-DIFFERENCE ALGORITHM 12.2

# To approximate the solution to the parabolic partial-differential

# equation subject to the boundary conditions

# u(0,t) = u(l,t) = 0, 0 < t < T = max t,

# and the initial conditions

# u(x,0) = F(x), 0 <= x <= l:

# INPUT: endpoint l: maximum time T: constant ALPHA: integers m, N.

# OUTPUT: approximations W(I,J) to u(x(I),t(J)) for each

# I = 1, ..., m-1 and J = 1, ..., N.

print(`This is the Backward-Difference Method for Heat Equation.`):

print(`Input the function F(X) in terms of x.`):

print(`For example: sin(3.141592654*x)`):

F := scanf(`%a`)[1]: print(`F(x) = `):print(F):

F := unapply(F,x):

print(`The lefthand endpoint on the X-axis is 0.`):

OK := FALSE:

while OK = FALSE do

print(`Input the righthand endpoint on the X-axis.`):

FX := scanf(`%f`)[1]: print(`Righthand endpoint = `):print(FX):

if FX <= 0 then

print(`Must be positive number.`):

else

OK := TRUE:

fi:

od:

OK := FALSE:

while OK = FALSE do

print(`Input the maximum value of the time variable T.`):

FT := scanf(`%f`)[1]: print(`Maximum time value = `):print(FT):

if FT <= 0 then

print(`Must be positive number.`):

else

OK := TRUE:

fi:

od:

print(`Input the constant alpha.`):

ALPHA := scanf(`%f`)[1]: print(`alpha = `):print(ALPHA):

OK := FALSE:

while OK = FALSE do

print(`Input integer m = number of intervals on X-axis`):

print(`and N = number of time intervals - separated by a blank.`):

print(`Note that m must be 3 or larger.`):

M := scanf(`%d`)[1]:

N := scanf(`%d`)[1]:

print(`Number of intervals on x-axis = `):print(M):

print(`Number of time intervals = `):print(N):

if M <= 2 or N <= 0 then

print(`Numbers are not within correct range.`):

else

OK := TRUE:

fi:

od:

if OK = TRUE then

M1 := M-1:

M2 := M-2:

N1 := N-1:

# Step 1

H := FX/M:

K := FT/N:

VV := ALPHA*ALPHA*K/(H*H):

# Step 2

for I2 from 1 to M1 do

W[I2-1] := F(I2*H):

od:

# Step 3

# Steps 3 - 11 solve a tridiagonal linear system using Algorithm 6.7

L[0] := 1+2*VV:

U[0] := -VV/L[0]:

# Step 4

for I2 from 2 to M2 do

L[I2-1] := 1+2*VV+VV*U[I2-2]:

U[I2-1] := -VV/L[I2-1]:

od:

# Step 5

L[M1-1] := 1+2*VV+VV*U[M2-1]:

# Step 6

for J from 1 to N do

# Step 7

# Current t

T := J*K:

Z[0] := W[0]/L[0]:

# Step 8

for I2 from 2 to M1 do

Z[I2-1] := (W[I2-1]+VV*Z[I2-2])/L[I2-1]:

od:

# Step 9

W[M1-1] := evalf(Z[M1-1]):

# Step 10

for I1 from 1 to M2 do

I2 := M2-I1+1:

W[I2-1] := evalf(Z[I2-1]-U[I2-1]*W[I2]):

od:

# Step 11

print(`Choice of output method:`):

print(`1. Output to screen`):

print(`2. Output to text file`):

print(`Please enter 1 or 2.`):

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, `BACKWARD-DIFFERENCE METHOD\134n\134n`):

fprintf(OUP, ` I X(I) W(X(I),%12.6e)\134n`, FT):

for I2 from 1 to M1 do

X := I2*H:

fprintf(OUP, `%3d %11.8f %14.8f\134n`, I2, X, W[I2-1]):

od:

fi:

if OUP <> default then

fclose(OUP):

print(`Output file `,NAME,` created successfully`):

fi:

SVpUaGlzfmlzfnRoZX5CYWNrd2FyZC1EaWZmZXJlbmNlfk1ldGhvZH5mb3J+SGVhdH5FcXVhdGlvbi5HNiI=

SUdJbnB1dH50aGV+ZnVuY3Rpb25+RihYKX5pbn50ZXJtc35vZn54Lkc2Ig==

SUFGb3J+ZXhhbXBsZTp+fnNpbigzLjE0MTU5MjY1NCp4KUc2Ig==

SShGKHgpfj1+RzYi

LUkkc2luRzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliRzYiNiMsJEkieEdGJyQiK2FFZlRKISIq

SUpUaGV+bGVmdGhhbmR+ZW5kcG9pbnR+b25+dGhlflgtYXhpc35pc34wLkc2Ig==

SUxJbnB1dH50aGV+cmlnaHRoYW5kfmVuZHBvaW50fm9ufnRoZX5YLWF4aXMuRzYi

STZSaWdodGhhbmR+ZW5kcG9pbnR+PX5HNiI=

JCIiIiIiIQ==

SVBJbnB1dH50aGV+bWF4aW11bX52YWx1ZX5vZn50aGV+dGltZX52YXJpYWJsZX5ULkc2Ig==

STZNYXhpbXVtfnRpbWV+dmFsdWV+PX5HNiI=

JCIiIiIiIQ==

STpJbnB1dH50aGV+Y29uc3RhbnR+YWxwaGEuRzYi

SSlhbHBoYX49fkc2Ig==

JCIiIiIiIQ==

SVBJbnB1dH5pbnRlZ2Vyfm1+PX5udW1iZXJ+b2Z+aW50ZXJ2YWxzfm9uflgtYXhpc0c2Ig==

SVlhbmR+Tn49fm51bWJlcn5vZn50aW1lfmludGVydmFsc34tfnNlcGFyYXRlZH5ieX5hfmJsYW5rLkc2Ig==

SUFOb3RlfnRoYXR+bX5tdXN0fmJlfjN+b3J+bGFyZ2VyLkc2Ig==

SUFOdW1iZXJ+b2Z+aW50ZXJ2YWxzfm9ufngtYXhpc349fkc2Ig==

IiM1

STxOdW1iZXJ+b2Z+dGltZX5pbnRlcnZhbHN+PX5HNiI=

IiNd

STlDaG9pY2V+b2Z+b3V0cHV0fm1ldGhvZDpHNiI=

STQxLn5PdXRwdXR+dG9+c2NyZWVuRzYi

STcyLn5PdXRwdXR+dG9+dGV4dH5maWxlRzYi

STVQbGVhc2V+ZW50ZXJ+MX5vcn4yLkc2Ig==

SSlJbnB1dH5pc0c2Ig==

IiIi

BACKWARD-DIFFERENCE METHOD I X(I) W(X(I),1.000000e+00) 1 0.10000000 0.00004051 2 0.20000000 0.00007705 3 0.30000000 0.00010605 4 0.40000000 0.00012467 5 0.50000000 0.00013108 6 0.60000000 0.00012467 7 0.70000000 0.00010605 8 0.80000000 0.00007705 9 0.90000000 0.00004051

JSFH