restart:

# CRANK-NICOLSON ALGORITHM 12.3

# To approximate the solution of 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 Crank-Nicolson Method.`):

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:

# Step 1

H := FX/M:

K := FT/N:

# VV is used in place of lambda

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

# Set V(M) to zero

V[M-1] := 0:

# Step 3

for I2 from 1 to M1 do

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

od:

# Step 3

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

L[0] := 1+VV:

U[0] := -VV/(2*L[0]):

# Step 4

for I2 from 2 to M2 do

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

U[I2-1] := -VV/(2*L[I2-1]):

od:

# Step 5

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

# Step 6

for J from 1 to N do

# Step 7

# Current t

T := J*K:

Z[0] := ((1-VV)*V[0]+VV*V[1]/2)/L[0]:

# Step 8

for I2 from 2 to M1 do

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

od:

# Step 9

V[M1-1] := Z[M1-1]:

# Step 10

for I1 from 1 to M2 do

I2 := M2-I1+1:

V[I2-1] := Z[I2-1]-U[I2-1]*V[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, `CRANK-NICOLSON 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 %13.8f\134n`, I2, X, V[I2-1]):

od:

if OUP <> default then

fclose(OUP):

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

fi:

SUNUaGlzfmlzfnRoZX5DcmFuay1OaWNvbHNvbn5NZXRob2QuRzYi

SUdJbnB1dH50aGV+ZnVuY3Rpb25+RihYKX5pbn50ZXJtc35vZn54Lkc2Ig==

SUJGb3J+ZXhhbXBsZTp+fn5zaW4oMy4xNDE1OTI2NTQqeClHNiI=

SShGKHgpfj1+RzYi

LUkkc2luRzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliRzYiNiMsJEkieEdGJyQiK2FFZlRKISIq

SUpUaGV+bGVmdGhhbmR+ZW5kcG9pbnR+b25+dGhlflgtYXhpc35pc34wLkc2Ig==

SUxJbnB1dH50aGV+cmlnaHRoYW5kfmVuZHBvaW50fm9ufnRoZX5YLWF4aXMuRzYi

STZSaWdodGhhbmR+ZW5kcG9pbnR+PX5HNiI=

JCIiIiIiIQ==

SVBJbnB1dH50aGV+bWF4aW11bX52YWx1ZX5vZn50aGV+dGltZX52YXJpYWJsZX5ULkc2Ig==

STZNYXhpbXVtfnRpbWV+dmFsdWV+PX5HNiI=

JCIiJiEiIg==

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==

SSpJTnB1dH5pc35HNiI=

IiIi

CRANK-NICOLSON METHOD I X(I) W(X(I),5.000000e-01) 1 0.10000000 0.00230512 2 0.20000000 0.00438461 3 0.30000000 0.00603489 4 0.40000000 0.00709444 5 0.50000000 0.00745954 6 0.60000000 0.00709444 7 0.70000000 0.00603489 8 0.80000000 0.00438461 9 0.90000000 0.00230512

JSFH