(*  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
*
*  OUPUT:  coefficients C(1), ..., C(n) of the basis functions
*)
simpson[fn_,a_,b_] := Module[{h,i,Y,z,s},
   h = (b-a)/4;
   For[i = 0,i <= 4, i++,
      Y = a+i*h;
      If[fn == 1,
	 z[i] = N[(4-i)*i*h^2*QQ[Y]];
      ];
      If[fn == 2,
	 z[i] = N[(i*h)^2*QQ[Y]];
      ];
      If[fn == 3,
	 z[i] = N[(h*(4-i))^2*QQ[Y]];
      ];
      If[fn == 4,
	 z[i] = N[P[Y]];
      ];
      If[fn == 5,
	 z[i] = N[i*h*F[Y]];
      ];
      If[fn == 6,
	 z[i] = N[(4-i)*h*F[Y]];
      ];
   ];
   s = (z[0]+z[4]+2*z[2]+4*(z[1]+z[3]))*h/3;
   s
];
TEMP = Input["This is the Piecewise Linear Rayleigh-Ritz Method.\n
   Input the function F(x) in terms of x for example,\n
   2*3.141592654^2*Sin[3.141592654*x]\n"];
F[x_] := Evaluate[TEMP];
TEMP = Input["Input the function Q(x) in terms of x \n
    for example, 3.141592654^2\n"];
QQ[x_] := Evaluate[TEMP];
TEMP = Input["Input the function P(x) in terms of x\n
    for example, 1\n"];
P[x_] := Evaluate[TEMP];
AA = InputString["X(0), ..., X(n+1) are to be supplied.\n
   Are the preparations complete?\n
   Answer 'yes' or 'no'\n"];
OK = 0;
If[AA == "y" || AA == "yes" || AA =="Y",
   OK = 0;
   While[OK == 0,
      n = Input["Input the integer n where X(0) = 0, X(n+1) = 1.\n"];
      If[n < 1,
	 Input["N must be greater than one.\n
	 \n
	 Press 1 [enter] to continue\n"],
	 OK = 1;
      ];
   ];
   X[0] = 0;
   X[n+1] = 1;
   FLAG = Input["Choice of method to input X(1), ..., X(n):\n
   1. Input from keyboard at the prompt\n
   2. Equally spaced nodes to be calculated\n
   3. Input from text file\n
   Please enter 1, 2, or 3.\n"];
   If[FLAG == 2,
      HC = 1/(n+1);
      For[J = 1 , J<= n, J++,
	 X[J] = N[J*HC];
	 H[J-1] = N[HC];
      ];
      H[n] = HC,
      If[FLAG == 3,
	 AA = InputString["Has the input file been created?
	 Enter 'yes' or 'no'\n"];
	 If[AA == "yes" || AA == "y" || AA == "Y",
	    NAME = InputString["Enter the input file name using the format\n
	    - drive:\ name.ext,  for example: A:\data.dta\n"];
	    INP = OpenRead[NAME];
	    For[J = 1, J <= n, J++,
	       X[J] = Evaluate[Read[INP,Number]];
	    ];
	    For[J = 0, J <= n, J++,
	       H[J] = X[J+1]-X[J];
	    ];
	    Close[INP],
	    Input["This program will end so that the input\n
	    file can be created.\n"];
	    OK = 0;
	 ],
	 For[J = 1, J <= n, J++,
	    X[J] = Input["Input X("<>ToString[J]<>")\n"];
	    H[J-1] = N[X[J]-X[J-1]];
	 ];
	 H[n] = X[n+1]-X[n];
      ];
   ],
   Input["This program will end so that the preparations\n
   can be completed. \n
   \n
   Press 1 [enter] to continue\n"];
   OK = 0;
];
(* Step 1 - Done within the input procedure *)
If[OK == 1,
   N1 = n-1;
   (* Step 3 *)
   For[J = 1, J <= N1, J++,
      Q[0,J-1] = simpson[1,X[J],X[J+1]]/H[J]^2;
      Q[1,J-1] = simpson[2,X[J-1],X[J]]/H[J-1]^2;
      Q[2,J-1] = simpson[3,X[J],X[J+1]]/H[J]^2;
      Q[3,J-1] = simpson[4,X[J-1],X[J]]/H[J-1]^2;
      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];
   ];
   Q[1,n-1] = simpson[2,X[n-1],X[n]]/((H[n-1])^2);
   Q[2,n-1] = simpson[3,X[n],X[n+1]]/((H[n])^2);
   Q[3,n-1] = simpson[4,X[n-1],X[n]]/((H[n-1])^2);
   Q[3,n] = simpson[4,X[n],X[n+1]]/((H[n])^2);
   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 4 *)
   For[J = 1, J <= n, J++,
      ALPHA[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];
   ];
   (* Step 5 *)
   ALPHA[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 
      Algorrithm 6.7 *)
   (* Step 6 *)
   A[0] = ALPHA[0];
   ZETA[0] = BETA[0]/ALPHA[0];
   Z[0] = B[0]/A[0];
   (* Step 7 *)
   For[J = 2, J <= N1, J++,
      A[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];
   ];
   (* Step 8 *)
   If[n > 1,
   A[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 9 *)
   c[n-1] = Z[n-1];
   For[J = 1, J <= N1, J++,
      J1 = n-J;
      c[J1-1] = Z[J1-1]-ZETA[J1-1]*c[J1];
   ];
   FLAG = Input["Select output destination\n
      1. Screen\n
      2. Text file\n
    Enter 1 or 2\n"];
   If[FLAG == 2,
      NAME = InputString["Input the file name\n
	 For example:   output.dta\n"];
      OUP = OpenWrite[NAME,FormatType->OutputForm],
      OUP = "stdout";
   ];
   Write[OUP," PIECEWISE LINEAR RAYLEIGH-RITZ METHOD\n"];
   Write[OUP,"\n"];
   Write[OUP,"i   X(i-1)      X(i)       X(i+1)     C(i)\n"];
   Write[OUP,"\n"];
   For[J = 1, J <= n, J++,
      Write[OUP,J,"   ",SetPrecision[X[J-1],10],"     ",
            SetPrecision[X[J],10],"   ",SetPrecision[X[J+1],10],
            "   ",SetPrecision[c[J-1],10]];
   ];
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP]
   ];
];
Quit[];