alg114.txt

来自「Numerical Anaysis 8th Edition Burden and」· 文本 代码 · 共 177 行

> restart;
> # NONLINEAR FINITE-DIFFERENCE ALGORITHM 11.4
> #
> # To approximate the solution to the nonlinear boundary-value problem
> #
> #    Y'' = F(X,Y,Y'), A<=X<=B, Y(A) = ALPHA, Y(B) = BETA:
> #
> # INPUT:   Endpoints A,B; boundary conditions ALPHA, BETA;
> #          integer N; tolerance TOL; maximum number of iterations M.
> #
> # OUTPUT:  Approximations W(I) TO Y(X(I)) for each I=0,1,...,N+1
> #          or a message that the maximum number of iterations was
> #          exceeded.
> alg114 := proc() local F, FY, FYP, OK, AA, BB, ALPHA, BETA, N, TOL, NN, FLAG, NAME, OUP, N1, H, I, W, K, X, T, A, B, D, C, L, U, Z, V, VMAX, J;
> printf(`This is the Nonlinear Finite-Difference Method.\n`);
> printf(`Input the function F(X,Y,Z) in terms of x, y, z.\n`);
> printf(`For example:   (32+2*x^3-y*z)/8\n`);
> F := scanf(`%a`)[1];
> FY := diff(F,y);
> FYP := diff(F,z);
> F := unapply(F,x,y,z);
> FY := unapply(FY,x,y,z);
> FYP := unapply(FYP,x,y,z);
> OK := FALSE;
> while  OK = FALSE do
> printf(`Input left and right endpoints separated by blank.\n`);
> AA := scanf(`%f`)[1];
> BB := scanf(`%f`)[1];
> if AA >= BB then
> printf(`Left endpoint must be less than right endpoint.\n`);
> else
> OK := TRUE;
> fi;
> od;
> printf(`Input Y(  %.10e).\n`, AA);
> ALPHA := scanf(`%f`)[1];
> printf(`Input Y(  %.10e).\n`, BB);
> BETA := scanf(`%f`)[1];
> OK := FALSE;
> while OK = FALSE do
> printf(`Input an integer > 1 for the number of\n`);
> printf(`subintervals.  Note that h := (b-a)/(n+1)\n`);
> N := scanf(`%d`)[1];
> if N <= 1 then
> printf(`Number must exceed 1.\n`);
> else
> OK := TRUE;
> fi;
> od;
> OK := FALSE;
> while  OK = FALSE do
> printf(`Input Tolerance.\n`);
> TOL := scanf(`%f`)[1];
> if TOL <= 0 then
> printf(`Tolerance must be positive.\n`);
> else
> OK := TRUE;
> fi;
> od;
> OK := FALSE;
> while  OK = FALSE do
> printf(`Input maximum number of iterations.\n`);
> NN := scanf(`%d`)[1];
> if NN <= 0 then
> printf(`Must be positive integer.\n`);
> else
> OK := TRUE;
> fi;
> od;
> if OK = TRUE then
> printf(`Choice of output method:\n`);
> printf(`1. Output to screen\n`);
> printf(`2. Output to text File\n`);
> printf(`Please enter 1 or 2.\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 2 then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example   A:\\OUTPUT.DTA\n`);
> NAME := scanf(`%s`)[1];
> OUP := fopen(NAME,WRITE,TEXT);
> else
> OUP := default;
> fi;
> fprintf(OUP, `NONLINEAR FINITE-DIFFERENCE METHOD\n\n`);
> fprintf(OUP, `  I    X(I)         W(I)\n`);
> # Step 1
> N1 := N-1;
> H := (BB-AA)/(N+1);
> # Step 2
> for I from 1 to N do
> W[I-1] := ALPHA+I*H*(BETA-ALPHA)/(BB-AA);
> od;
> # Step 3
> K := 1;
> # Step 4
> while K <= NN and OK = TRUE do
> # Step 5
> X := AA+H;
> T := (W[1]-ALPHA)/(2*H);
> A[0] := 2+H*H*FY(X,W[0],T);
> B[0] := -1+H*FYP(X,W[0],T)/2;
> D[0] := -(2*W[0]-W[1]-ALPHA+H*H*F(X,W[0],T));
> # Step 6
> for I from 2 to N1 do
> X := AA+I*H;
> T := (W[I]-W[I-2])/(2*H);
> A[I-1] := 2+H*H*FY(X,W[I-1],T);
> B[I-1] := -1+H*FYP(X,W[I-1],T)/2;
> C[I-1] := -1-H*FYP(X,W[I-1],T)/2;
> D[I-1] := -(2*W[I-1]-W[I]-W[I-2]+H*H*F(X,W[I-1],T));
> od;
> # Step 7
> X := BB - H;
> T := (BETA-W[N-2])/(2*H);
> A[N-1] := 2+H*H*FY(X,W[N-1],T);
> C[N-1] := -1-H*FYP(X,W[N-1],T)/2;
> D[N-1] := -(2*W[N-1]-W[N-2]-BETA+H*H*F(X,W[N-1],T));
> # Step 8
> # Steps 8 - 12 solve a tridiagonal linear system using Algorithm 6.7
> L[0] := A[0];
> U[0] := B[0]/A[0];
> Z[0] := D[0]/L[0];
> # Step 9
> for I from 2 to N1 do
> L[I-1] := A[I-1]-C[I-1]*U[I-2];
> U[I-1] := B[I-1]/L[I-1];
> Z[I-1] := (D[I-1]-C[I-1]*Z[I-2])/L[I-1];
> od;
> # Step 10
> L[N-1] := A[N-1]-C[N-1]*U[N-2];
> Z[N-1] := (D[N-1]-C[N-1]*Z[N-2])/L[N-1];
> # Step 11
> V[N-1] := Z[N-1];
> VMAX := abs(V[N-1]);
> W[N-1] := W[N-1]+V[N-1];
> # Step 12
> for J from 1 to N1 do
> I := N-J;
> V[I-1] := Z[I-1]-U[I-1]*V[I];
> W[I-1] := W[I-1]+V[I-1];
> if abs(V[I-1]) > VMAX then
> VMAX := abs(V[I-1]);
> fi;
> od;
> # Step 13
> # Test for accuracy
> if VMAX <= TOL then
> I := 0;
> fprintf(OUP, `%3d %13.8f %13.8f\n`, I, AA, ALPHA);
> # Step 14
> for I from 1 to N do
> X := AA+I*H;
> fprintf(OUP, `%3d %13.8f %13.8f\n`, I, X, W[I-1]);
> od;
> I := N+1;
> # Step 15
> fprintf(OUP, `%3d %13.8f %13.8f\n`, I, BB, BETA);
> fprintf(OUP, `Convergence in %d iterations\n`, K);
> OK := FALSE;
> else
> # Step 16
> K := K+1;
> fi;
> od;
> # Step 17
> if K > NN then
> fprintf(OUP, `No convergence in %d iterations\n`, NN);
> fi;
> fi;
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> RETURN(0);
> end;
> alg114();