restart:# PIECEWISE LINEAR RAYLEIGH-RITZ ALGORITHM 11.5## To approximate the solution of the boundary-value problem## -D(P(X)Y')/DX + Q(X)Y = F(X), 0 <= X <= 1,# Y(0) = Y(1) = 0,## with a piecewise linear function:## INPUT: integer N: mesh points X(0) = 0 < X(1) < ...# < X(N) < X(N+1) = 1## OUTPUT: coefficients C(1),...,C(N) of the basis functionsSIMPSON := proc(FN,A,B) local H, I1, Y, Z, simpson:H := (B-A)/4:for I1 from 0 to 4 doY := A+I1*H:if FN = 1 thenZ[I1] := (4-I1)*I1*H*H*QQ(Y):fi:if FN = 2 thenZ[I1] := (I1*H)*(I1*H)*QQ(Y):fi:if FN = 3 thenZ[I1] := (H*(4-I1))*(H*(4-I1))*QQ(Y):fi:if FN = 4 thenZ[I1] := P(Y):fi:if FN = 5 thenZ[I1] := I1*H*F(Y):fi:if FN = 6 thenZ[I1] := (4-I1)*H*F(Y):fi:od:simpson := (Z[0]+Z[4]+2*Z[2]+4*(Z[1]+Z[3]))*H/3:RETURN(simpson):end: print(`This is the Piecewise Linear Rayleigh-Ritz Method.`):print(`Input F(X), Q(X), and P(X) in terms of x separated by a space.`):print(`For example:`): print(`2*3.141592654^2*sin(3.141592654*x) 3.141592654^2 1`):F := scanf(`%a`)[1]:QQ := scanf(`%a`)[1]:P := scanf(`%a`)[1]: print(`F(x) = `):print(F):print(`Q(x) = `):print(QQ):print(`P(x) = `):print(P):F := unapply(F,x):QQ := unapply(QQ,y):P := unapply(P,z):print(` `):print(`X(0), ..., X(N+1) are to be supplied.`):print(` `):print(`Are the preparations complete? Answer 1 for yes or 2 for no.`):print(` `):AA := scanf(`%d`)[1]:OK := FALSE:if AA = 1 thenOK := FALSE:while OK = FALSE doprint(`Input integer N where X(0) = 0, X(N+1) = 1.`):N := scanf(`%d`)[1]:print(`N = `):print(N):if N <= 1 thenprint(`N must be greater than one.`):elseOK := TRUE:fi:od:X[0] := 0:X[N+1] := 1:print(`Choice of method to input X(1), ..., X(N):`):print(`1. Input from keyboard at the prompt`):print(`2. Equally spaced nodes to be calculated`):print(`3. Input from text file`):print(`Please enter 1, 2, or 3.`):FLAG := scanf(`%d`)[1]: print(`Input is `):print(FLAG):if FLAG = 2 thenHC := 1/(N+1):for J from 1 to N doX[J] := evalf(J*HC):H[J-1] := HC:od:H[N] := HC:elseif FLAG = 3 thenprint(`Has the input file been created? `):print(`Enter 1 for yes or 2 for no.`):AA := scanf(`%d`)[1]: print(`Input is `):print(AA):if AA = 1 thenprint(`Enter the input file name using the format`):print(` - drive:\134\134name.ext,`):print(`for example: A:\134\134DATA.DTA`):NAME := scanf(`%s`)[1]: print(`Input file is `):print(NAME):INP := fopen(NAME,READ,TEXT):for J from 1 to N doX[J] := evalf(fscanf(INP, `%f`)[1]):od:for J from 0 to N doH[J] := X[J+1]-X[J]:od:fclose(INP):elseprint(`The program will end so that the input `):print(`file can be created.`):OK := FALSE:fi:elsefor J from 1 to N doprint(`Input X(j) for j =`):print(J):X[J] := evalf(scanf(`%f`)[1]): print(`Input is `):print(X[J]):H[J-1] := X[J]-X[J-1]:od:H[N] := X[N+1]-X[N]:fi:fi:elseprint(`The program will end so that the preparations`):print(`can be completed.`):OK := FALSE:fi:# Step 1 is done in the input procedure.if OK = TRUE thenN1 := N-1:# Step 3for J from 1 to N1 doQ[0,J-1] := SIMPSON(1,X[J],X[J+1])/((H[J])*(H[J])):Q[1,J-1] := SIMPSON(2,X[J-1],X[J])/((H[J-1])*(H[J-1])):Q[2,J-1] := SIMPSON(3,X[J],X[J+1])/((H[J])*(H[J])):Q[3,J-1] := SIMPSON(4,X[J-1],X[J])/((H[J-1])*(H[J-1])):Q[4,J-1] := SIMPSON(5,X[J-1],X[J])/H[J-1] :Q[5,J-1] := SIMPSON(6,X[J],X[J+1])/H[J]:od:Q[1,N-1] := SIMPSON(2,X[N-1],X[N])/((H[N-1])*(H[N-1])):Q[2,N-1] := SIMPSON(3,X[N],X[N+1])/((H[N])*(H[N])):Q[3,N-1] := SIMPSON(4,X[N-1],X[N])/((H[N-1])*(H[N-1])):Q[3,N] := SIMPSON(4,X[N],X[N+1])/((H[N])*(H[N])):Q[4,N-1] := SIMPSON(5,X[N-1],X[N])/H[N-1]:Q[5,N-1] := SIMPSON(6,X[N],X[N+1])/H[N]:# Step 4for J from 1 to N1 doALPHA[J-1] := Q[3,J-1]+Q[3,J]+Q[1,J-1]+Q[2,J-1]:BETA[J-1] := Q[0,J-1]-Q[3,J]:B[J-1] := Q[4,J-1]+Q[5,J-1]:od:# Step 5ALPHA[N-1] := Q[3,N-1]+Q[3,N]+Q[1,N-1]+Q[2,N-1]:B[N-1] := Q[4,N-1]+Q[5,N-1]:# Steps 6 - 10 solve a tridiagonal linear system using Algorithm 6.7# Step 6A[0] := ALPHA[0]:ZETA[0] := BETA[0]/ALPHA[0]:Z[0] := B[0]/A[0]:# Step 7for J from 2 to N1 doA[J-1] := ALPHA[J-1]-BETA[J-2]*ZETA[J-2]:ZETA[J-1] := BETA[J-1]/A[J-1]:Z[J-1] := (B[J-1]-BETA[J-2]*Z[J-2])/A[J-1]:od:# Step 8A[N-1] := ALPHA[N-1] - BETA[N-2] * ZETA[N-2]:Z[N-1] := (B[N-1]-BETA[N-2]*Z[N-2])/A[N-1]:# Step 9C[N-1] := evalf(Z[N-1]):for J from 1 to N1 doJ1 := N - J:C[J1-1] := evalf(Z[J1-1]-ZETA[J1-1]*C[J1]):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 thenprint(`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):elseOUP := default:fi:fprintf(OUP, `PIECEWISE LINEAR RAYLEIGH-RITZ METHOD\134n\134n`):fprintf(OUP, ` I X(I-1) X(I) X(I+1) C(I)\134n\134n`):for J from 1 to N dofprintf(OUP, `%3d %11.8f %11.8f %11.8f %13.8f\134n`, J, X[J-1],X[J], X[J+1], C[J-1]):od:if OUP <> default thenfclose(OUP):print(`Output file `,NAME,` created successfully`):fi:fi:JSFH