> 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:
> 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:\\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, `CRANK-NICOLSON METHOD\n\n`):
> fprintf(OUP, `  I     X(I)         W(X(I),%12.6e)\n`, FT):
> for I2 from 1 to M1 do
> X := I2*H:
> fprintf(OUP, `%3d %11.8f %13.8f\n`, I2, X, V[I2-1]):
> od:
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi:
> fi: