%% MAIN is the main program for FEM1D_PMETHOD.
%
%  Discussion:
%
%    FEM1D_PMETHOD implements the P-version of the finite element method.
%
%    Program to solve the one dimensional problem:
%
%      - d/dX (P dU/dX) + Q U  =  F
%
%    by the finite-element method using a sequence of polynomials
%    which satisfy the boundary conditions and are orthogonal
%    with respect to the inner product:
%
%      (U,V)  =  Integral (-1 to 1) P U' V' + Q U V dx
%
%    Here U is an unknown scalar function of X defined on the
%    interval [-1,1], and P, Q and F are given functions of X.
%
%    The boundary values are U(-1) = U(1)=0.
%
%    Sample problem #1:
%
%      U=1-x**4,        P=1, Q=1, F=1.0+12.0*x**2-x**4
%
%    Sample problem #2:
%
%      U=cos(0.5*pi*x), P=1, Q=0, F=0.25*pi*pi*cos(0.5*pi*x)
%
%    The program should be able to get the exact solution for
%    the first problem, using NP = 2.
%
%  Modified:
%
%    03 November 2006
%
%  Author:
%
%    Max Gunzburger
%    Teresa Hodge
%
%    MATLAB translation by John Burkardt
%
%  Parameters:
%
%    real A(0:NP), the squares of the norms of the 
%    basis functions.
%
%    real ALPHA(NP).
%    ALPHA(I) contains one of the coefficients of a recurrence
%    relationship that defines the basis functions.
%
%    real BETA(NP).
%    BETA(I) contains one of the coefficients of a recurrence
%    relationship that defines the basis functions.
%
%    real F(1:NP+1).
%    F contains the basis function coefficients that form the
%    representation of the solution U.  That is,
%      U(X)  =  SUM (I=0 to NP) F(I+1) * BASIS(I)(X)
%    where "BASIS(I)(X)" means the I-th basis function
%    evaluated at the point X.
%
%    integer NP.
%    The highest degree polynomial to use.
%
%    integer NPRINT.
%    The number of points at which the computed solution
%    should be printed out at the end of the computation.
%
%    integer PROBLEM, indicates the problem being solved.
%    1, U=1-x**4, P=1, Q=1, F=1.0+12.0*x**2-x**4.
%    2, U=cos(0.5*pi*x), P=1, Q=0, F=0.25*pi*pi*cos(0.5*pi*x).
%
%    integer QUAD_NUM, the order of the quadrature rule.
%
%    real QUAD_W(QUAD_NUM), the quadrature weights.
%
%    real QUAD_X(QUAD_NUM), the quadrature abscissas.
%
  clear

  np = 2;
  quad_num = 10;
  nprint = 10;
  problem = 2;

  timestamp;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_PMETHOD\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve the two-point boundary value problem\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  - d/dX (P dU/dX) + Q U  =  F\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  on the interval [-1,1], with\n' );
  fprintf ( 1, '  U(-1) = U(1) = 0.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The P method is used, which represents U as\n' );
  fprintf ( 1, '  a weighted sum of orthogonal polynomials.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Highest degree polynomial to use is %d\n', np );
  fprintf ( 1, '  Number of points to be used for output = %d\n', nprint );

  if ( problem == 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Problem #1:\n' );
    fprintf ( 1, '  U=1-x**4,\n' );
    fprintf ( 1, '  P=1,\n' );
    fprintf ( 1, '  Q=1,\n' );
    fprintf ( 1, '  F=1 + 12 * x**2 - x**4\n' );
  elseif ( problem == 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Problem #2:\n' );
    fprintf ( 1, '  U=cos(0.5*pi*x),\n' );
    fprintf ( 1, '  P=1,\n' );
    fprintf ( 1, '  Q=0,\n' );
    fprintf ( 1, '  F=0.25*pi*pi*cos(0.5*pi*x)\n' );
  end
%
%  Get quadrature abscissas and weights for interval [-1,1].
%
  [ quad_w, quad_x ] = quad ( quad_num );
%
%  Compute the constants for the recurrence relationship
%  that defines the basis functions.
%
  [ a, alpha, beta ] = alpbet ( np, problem, quad_num, quad_w, quad_x );
%
%  Test the orthogonality of the basis functions.
%
  ortho ( a, alpha, beta, np, problem, quad_num, quad_w, quad_x );
%
%  Solve for the solution of the problem, in terms of coefficients
%  of the basis functions.
%
  f = sol ( a, alpha, beta, np, problem, quad_num, quad_w, quad_x );
%
%  Print out the solution, evaluated at each of the NPRINT points.
%
  out ( alpha, beta, f, np, nprint );
%
%  Compare the computed and exact solutions.
%
  exact ( alpha, beta, f, np, nprint, problem, quad_num, quad_w, quad_x );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PMETHOD\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp;