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. print(`This is the Nonlinear Finite-Difference Method.`): print(`Input the function F(X,Y,Z) in terms of x, y, z.`): print(`For example: (32+2*x^3-y*z)/8`): F := scanf(`%a`)[1]: print(`F(x,y,z) = `): print(F): 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 print(`Input left and right endpoints separated by blank.`): AA := scanf(`%f`)[1]: BB := scanf(`%f`)[1]: print(`Left endpoint =`): print(AA): print(`Right endpoint =`): print(BB): if AA >= BB then print(`Left endpoint must be less than right endpoint.`): else OK := TRUE: fi: od: print(`Input Y(a).`): ALPHA := scanf(`%f`)[1]: print(`y(a) =`): print(ALPHA): print(`Input Y(b).`): BETA := scanf(`%f`)[1]: print(`y(b) =`): print(BETA): OK := FALSE: while OK = FALSE do print(`Input an integer > 1 for the number of`): print(`subintervals. Note that h := (b-a)/(n+1)`): N := scanf(`%d`)[1]: print(`Number of subintervals =`): print(N): if N <= 1 then print(`Number must exceed 1.`): else OK := TRUE: fi: od: OK := FALSE: while OK = FALSE do print(`Input 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 maximum number of iterations.`): NN := scanf(`%d`)[1]: print(`Maximum number of iterations =`): print(NN): if NN <= 0 then print(`Must be positive integer.`): else OK := TRUE: fi: od: 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, `NONLINEAR FINITE-DIFFERENCE METHOD\134n\134n`): fprintf(OUP, ` I X(I) W(I)\134n`): # Step 1 N1 := N-1: H := (BB-AA)/(N+1): # Step 2 for I2 from 1 to N do W[I2-1] := ALPHA+I2*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: D2[0] := -(2*W[0]-W[1]-ALPHA+H*H*F(X,W[0],T)): # Step 6 for I2 from 2 to N1 do X := AA+I2*H: T := (W[I2]-W[I2-2])/(2*H): A[I2-1] := 2+H*H*FY(X,W[I2-1],T): B[I2-1] := -1+H*FYP(X,W[I2-1],T)/2: C[I2-1] := -1-H*FYP(X,W[I2-1],T)/2: D2[I2-1] := -(2*W[I2-1]-W[I2]-W[I2-2]+H*H*F(X,W[I2-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: D2[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] := D2[0]/L[0]: # Step 9 for I2 from 2 to N1 do L[I2-1] := A[I2-1]-C[I2-1]*U[I2-2]: U[I2-1] := B[I2-1]/L[I2-1]: Z[I2-1] := (D2[I2-1]-C[I2-1]*Z[I2-2])/L[I2-1]: od: # Step 10 L[N-1] := A[N-1]-C[N-1]*U[N-2]: Z[N-1] := (D2[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 I2 := N-J: V[I2-1] := Z[I2-1]-U[I2-1]*V[I2]: W[I2-1] := W[I2-1]+V[I2-1]: if abs(V[I2-1]) > VMAX then VMAX := abs(V[I2-1]): fi: od: # Step 13 # Test for accuracy if VMAX <= TOL then I2 := 0: fprintf(OUP, `%3d %13.8f %13.8f\134n`, I2, AA, ALPHA): # Step 14 for I2 from 1 to N do X := AA+I2*H: fprintf(OUP, `%3d %13.8f %13.8f\134n`, I2, X, W[I2-1]): od: I2 := N+1: # Step 15 fprintf(OUP, `%3d %13.8f %13.8f\134n`, I2, BB, BETA): fprintf(OUP, `Convergence in %d iterations\134n`, K): OK := FALSE: else # Step 16 K := K+1: fi: od: # Step 17 if K > NN then fprintf(OUP, `No convergence in %d iterations\134n`, NN): fi: if OUP <> default then fclose(OUP): print(`Output file `,NAME,` created successfully`): fi:

SVBUaGlzfmlzfnRoZX5Ob25saW5lYXJ+RmluaXRlLURpZmZlcmVuY2V+TWV0aG9kLkc2Ig==

SVFJbnB1dH50aGV+ZnVuY3Rpb25+RihYLFksWil+aW5+dGVybXN+b2Z+eCx+eSx+ei5HNiI=

SUBGb3J+ZXhhbXBsZTp+fn4oMzIrMip4XjMteSp6KS84RzYi

SSxGKHgseSx6KX49fkc2Ig==

LCgiIiUiIiIqJEkieEc2IiIiJCNGJEYjKiZJInlHRidGJEkiekdGJ0YkIyEiIiIiKQ==

SVNJbnB1dH5sZWZ0fmFuZH5yaWdodH5lbmRwb2ludHN+c2VwYXJhdGVkfmJ5fmJsYW5rLkc2Ig==

STBMZWZ0fmVuZHBvaW50fj1HNiI=

JCIiIiIiIQ==

STFSaWdodH5lbmRwb2ludH49RzYi

JCIiJCIiIQ==

SSxJbnB1dH5ZKGEpLkc2Ig==

SSd5KGEpfj1HNiI=

JCIkcSIhIiI=

SSxJbnB1dH5ZKGIpLkc2Ig==

SSd5KGIpfj1HNiI=

JCIrTExMTDkhIik=

SUdJbnB1dH5hbn5pbnRlZ2Vyfj5+MX5mb3J+dGhlfm51bWJlcn5vZkc2Ig==

SUpzdWJpbnRlcnZhbHMufn5Ob3RlfnRoYXR+aH46PX4oYi1hKS8obisxKUc2Ig==

STlOdW1iZXJ+b2Z+c3ViaW50ZXJ2YWxzfj1HNiI=

IiM+

STFJbnB1dH5Ub2xlcmFuY2UuRzYi

SSxUb2xlcmFuY2V+PUc2Ig==

JCIiIiEiJg==

SURJbnB1dH5tYXhpbXVtfm51bWJlcn5vZn5pdGVyYXRpb25zLkc2Ig==

ST9NYXhpbXVtfm51bWJlcn5vZn5pdGVyYXRpb25zfj1HNiI=

IiM6

STlDaG9pY2V+b2Z+b3V0cHV0fm1ldGhvZDpHNiI=

STQxLn5PdXRwdXR+dG9+c2NyZWVuRzYi

STcyLn5PdXRwdXR+dG9+dGV4dH5GaWxlRzYi

STVQbGVhc2V+ZW50ZXJ+MX5vcn4yLkc2Ig==

SSlJbnB1dH5pc0c2Ig==

IiIi

NONLINEAR FINITE-DIFFERENCE METHOD I X(I) W(I) 0 1.00000000 17.00000000 1 1.10000000 15.75450253 2 1.20000000 14.77173965 3 1.30000000 13.99567744 4 1.40000000 13.38629656 5 1.50000000 12.91425241 6 1.60000000 12.55753823 7 1.70000000 12.29932628 8 1.80000000 12.12652886 9 1.90000000 12.02881381 10 2.00000000 11.99791542 11 2.10000000 12.02714237 12 2.20000000 12.11101980 13 2.30000000 12.24502486 14 2.40000000 12.42538836 15 2.50000000 12.64894403 16 2.60000000 12.91301262 17 2.70000000 13.21531175 18 2.80000000 13.55388508 19 2.90000000 13.92704612 20 3.00000000 14.33333333 Convergence in 4 iterations

JSFH