program pmethod
!
!*******************************************************************************
!
!! PMETHOD implements the P-version of the finite element method.
!
!
!  Programmed by Max Gunzburger and Teresa Hodge.
!
!
!  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.
!
!
!  Two sample problems:
!
!    U = 1-x**4,    P=1, Q=1, F=1.0+12.0*x**2-x**4
!
!    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.
!
!  Parameters:
!
!    REAL A(0:NP).
!    A(I) contains the square of the norm of the I-th basis
!    function.
!
!    REAL ALPHA(NP).
!    ALPHA(I) contains one of the coefficients of a recurrence
!    relationship that defines the basis functions.
!
!    REAL B(0:NP,0:NP).
!    B(I,J) contains the dot product of the I-th and J-th
!    basis functions.  The basis functions should be orthogonal,
!    so B(I,J) should be 0 unless I equals J.
!
!    REAL BETA(NP).
!    BETA(I) contains one of the coefficients of a recurrence
!    relationship that defines the basis functions.
!
!    REAL F(0:NP).
!    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) * 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 QUAD_NUM is the number of points used in the quadrature rule.
!
!    REAL QUAD_W(QUAD_NUM) are the quadrature weights.
!
!    REAL QUAD_X(QUAD_NUM) are the quadrature abscissas.
!
  implicit none
!
  integer, parameter :: np = 2
  integer, parameter :: nprint = 10
  integer, parameter :: quad_num = 10
!
  real a(0:np)
  real alpha(np)
  real b(0:np,0:np)
  real beta(np)
  real f(0:np)
  real quad_w(quad_num)
  real quad_x(quad_num)
!
  call timestamp ( )

  print *, ' '
  print *, 'PMETHOD'
  print *, ' '
  print *, 'Solve the two-point boundary value problem'
  print *, ' '
  print *, '  - d/dX (P dU/dX) + Q U  =  F'
  print *, ' '
  print *, 'on the interval [-1,1], with'
  print *, 'U(-1) = U(1)=0.'
  print *, ' '
  print *, 'The P method is used, which represents U as'
  print *, 'a weighted sum of orthogonal polynomials.'
  print *, ' '
  print *, ' '
  print *, 'Highest degree polynomial to use is ', np
  print *, 'Number of points to be used for output = ', nprint
!
!  Get quadrature abscissas and weights for interval [-1,1].
!
  call quad ( quad_num, quad_w, quad_x )
!
!  Compute the constants for the recurrence relationship
!  that defines the basis functions.
!
  call alpbet ( a, alpha, beta, np, quad_num, quad_w, quad_x )
!
!  Test the orthogonality of the basis functions.
!
  call ortho ( a, alpha, b, beta, np, quad_num, quad_w, quad_x )
!
!  Solve for the solution of the problem, in terms of coefficients
!  of the basis functions.
!
  call sol ( a, alpha, beta, f, np, quad_num, quad_w, quad_x )
!
!  Print out the solution, evaluated at each of the NPRINT points.
!
  call out ( alpha, beta, f, np, nprint )
!
!  Compare the computed and exact solutions.
!
  call exact ( alpha, beta, f, np, nprint, quad_num, quad_w, quad_x )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'PMETHOD'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine alpbet ( a, alpha, beta, np, quad_num, quad_w, quad_x )
!
!*******************************************************************************
!
!! ALPBET calculates the constants in the recurrence relationship.
!
!
!  Discussion:
!
!    ALPHA and BETA are the constants in the three
!    term recurrence relation for the orthogonal basis functions
!    on [-1,1].
!
!    The routine also calculates A, the square of the norm of each basis
!    function.
!
  implicit none
!
  integer np
  integer quad_num
!
  real a(0:np)
  real alpha(np)
  real beta(np)
  integer i
  integer iq
  integer k
  real pp
  real q
  real qm1
  real qm1x
  real qm2
  real qm2x
  real qq
  real quad_w(quad_num)
  real quad_x(quad_num)
  real qx
  real s
  real ss
  real su
  real sv
  real t
  real u
  real v
  real x
!
  ss = 0.0E+00
  su = 0.0E+00
  do iq = 1, quad_num
    x = quad_x(iq)
    s =  4.0*pp(x)*x**2 + qq(x)*(1.0-x*x)**2
    u =  2.0*pp(x)*x*(3.0*x**2-1.0) + x*qq(x)*(1.0-x*x)**2
    ss = ss + s*quad_w(iq)
    su = su + u*quad_w(iq)
  end do

  beta(1) = 0.0E+00
  a(0) = ss
  alpha(1) = su/ss
  do i = 2, np+1
    ss = 0.0E+00
    su = 0.0E+00
    sv = 0.0E+00
    do iq = 1, quad_num
      x = quad_x(iq)
      q = 1.0
      qm1 = 0.0
      qx = 0.0
      qm1x = 0.0
      do k = 1,i-1
        qm2 = qm1
        qm1 = q
        qm2x = qm1x
        qm1x = qx
        q = (x-alpha(k))*qm1-beta(k)*qm2
        qx = qm1+(x-alpha(k))*qm1x-beta(k)*qm2x
      end do
      t = 1.0-x**2
      s = pp(x)*(t*qx-2.0*x*q)**2+qq(x)*(t*q)**2
      u = pp(x)*(x*t*qx+(1.-3.0*x**2)*q)*(t*qx-2.0*x*q)+x*qq(x)*(t*q)**2
      v = pp(x)*(x*t*qx+(1.0-3.0*x**2)*q)*(t*qm1x-2.0*x*qm1)+x*qq(x)*t**2*q*qm1
      ss = ss+s*quad_w(iq)
      su = su+u*quad_w(iq)
      sv = sv+v*quad_w(iq)
    end do

    a(i-1) = ss

    if ( i <= np ) then
      alpha(i) = su/ss
      beta(i) = sv/a(i-2)
    end if

  end do

  return
end
subroutine exact ( alpha, beta, f, np, nprint, quad_num, quad_w, quad_x )
!
!*******************************************************************************
!
!! EXACT compares the computed and exact solutions.
!
  implicit none
!
  integer np
  integer quad_num
!
  real alpha(np)
  real beta(np)
  real big_l2
  real error
  real f(0:np)
  integer i
  integer ip
  integer j
  integer k
  integer nprint
  integer, parameter :: nsub = 10
  real phii
  real phiix
  real quad_w(quad_num)
  real quad_x(quad_num)
  real ue
  real uex
  real up
  real x
  real xl
  real xr
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Comparison of computed and exact solutions:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    X        U computed    U exact     Difference'
  write ( *, '(a)' ) ' '

  do i = 0, nprint
    x = real(2*i-nprint)/real(nprint)
    ue = uex(x)
    up = 0.0E+00
    do j = 0, np
      call phi(alpha,beta,j,np,phii,phiix,x)
      up = up+phii*f(j)
    end do
    write(*,'(f8.4,3g14.6)') x,up,ue, ue - up
  end do
!
!  Compute the big L2 error.
!
  big_l2 = 0.0E+00

  do i = 1, nsub

    xl = real ( 2 * i - nsub - 1 ) / real ( nsub )
    xr = real ( 2 * i - nsub	 ) / real ( nsub )

    do j = 1, quad_num

      x = ( xl * ( 1.0E+00 - quad_x(j) ) &
          + xr * ( 1.0E+00 + quad_x(j) ) ) / 2.0E+00

      up = 0.0E+00
      do k = 0, np
        call phi ( alpha, beta, k, np, phii, phiix, x )
        up = up + phii * f(k)
      end do

      big_l2 = big_l2 + ( up - uex(x) )**2 * quad_w(j) * ( xr - xl ) / 2.0E+00

    end do

  end do

  big_l2 = sqrt ( big_l2 )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) 'Big L2 error = ', big_l2

  return
end
function ff ( x )
!
!*******************************************************************************
!
!! FF evaluates the right hand side function F(X) at any point X.
!
  implicit none
!
  real ff
  real, parameter :: pi = 3.14159265E+00
  real x
!
!  Test problem 1
!
!     ff = 1.0+12.0*x**2-x**4
!
!  Test problem 2
!
  ff = 0.25 * pi**2 * cos ( 0.5 * pi * x )

  return
end
subroutine ortho ( a, alpha, b, beta, np, quad_num, quad_w, quad_x )
!
!*******************************************************************************
!
!! ORTHO tests the basis functions for orthogonality.
!
  implicit none
!
  integer np
  integer quad_num
!
  real a(0:np)
  real alpha(np)
  real b(0:np,0:np)
  real beta(np)
  real bij
  integer i
  integer iq
  integer j
  real phii
  real phiix
  real phij
  real phijx
  real pp
  real qq
  real quad_w(quad_num)
  real x
  real quad_x(quad_num)
!
!  Zero out the B array, so we can start summing up the dot products.
!
  b(0:np,0:np) = 0.0E+00
!
!  Approximate the integral of the product of basis function
!  I and basis function J over the interval [-1,1].
!
!  We expect to get zero, except when I and J are equal,
!  when we should get A(I).
!
  do iq = 1,quad_num
    x = quad_x(iq)
    do i = 0,np
      call phi(alpha,beta,i,np,phii,phiix,x)
      do j = 0,np
        call phi(alpha,beta,j,np,phij,phijx,x)
        bij = pp(x)*phiix*phijx+qq(x)*phii*phij
        b(i,j) = b(i,j)+bij*quad_w(iq)
      end do
    end do
  end do
!
!  Print out the results of the test.
!
  print *, ' '
  print *, 'Basis function orthogonality test:'
  print *, ' '
  print *, '   i   j     b(i,j)/a(i)'
  print *, ' '
  do i = 0, np
    print *, ' '
    do j = 0, np
      write(*,'(2i4,g14.6)') i,j,b(i,j)/a(i)
    end do
  end do

  return
end
subroutine out ( alpha, beta, f, np, nprint )
!
!*******************************************************************************
!
!! OUT prints out the computed solution.
!
  implicit none
!
  integer np
!
  real alpha(np)
  real beta(np)
  real f(0:np)
  integer i
  integer ip
  integer nprint
  real phii
  real phiix
  real up
  real x
!
  print *, ' '
  print *, 'Representation of solution:'
  print *, ' '
  print *, 'Basis function coefficients:'
  print *, ' '
  do i = 0, np
    write(*,'(i4,g14.6)')i,f(i)
  end do

  print *, ' '
  print *, ' '
  print *, '       X     Approximate Solution'
  print *, ' '
  do ip = 0, nprint
    x = real(2*ip-nprint)/real(nprint)
    up = 0.0E+00
    do i = 0, np
      call phi(alpha,beta,i,np,phii,phiix,x)
      up = up+phii*f(i)
    end do
    write(*,'(2g14.6)')x,up
  end do

  print *, ' '

  return
end
subroutine phi ( alpha, beta, i, np, phii, phiix, x )
!
!*******************************************************************************
!
!! PHI evaluates the I-th basis function at the point X.
!
!  I can be between 0 and NP.
!
  implicit none
!
  integer np
!
  real alpha(np)
  real beta(np)
  integer i
  integer j
  real phii
  real phiix
  real q
  real qm1
  real qm1x
  real qm2
  real qm2x
  real qx
  real t
  real x
!
  qm1 = 0.0E+00
  q = 1.0E+00
  qm1x = 0.0E+00
  qx = 0.0E+00
  do j = 1, i
    qm2 = qm1
    qm1 = q
    qm2x = qm1x
    qm1x = qx
    t = x-alpha(j)
    q = t*qm1-beta(j)*qm2
    qx = qm1+t*qm1x-beta(j)*qm2x
  end do

  t = 1.0-x**2
  phii = t*q
  phiix = t*qx-2.0*x*q

  return
end
function pp ( x )
!
!*******************************************************************************
!
!! PP returns the value of the coefficient function P(X).
!
  implicit none
!
  real pp
  real x
!
!  Test problem 1
!
!     pp = 1.0E+00
!
!  Test problem 2
!
  pp = 1.0E+00

  return
end
subroutine quad ( quad_num, quad_w, quad_x )
!
!*******************************************************************************
!
!! QUAD returns the abscissas and weights for gaussian quadrature on [-1,1].
!
  implicit none
!
  integer quad_num
!
  real quad_w(quad_num)
  real quad_x(quad_num)
!
!  Quadrature points on [-1,1]
!
  quad_x(1) = -0.973906528517172
  quad_x(2) = -0.865063366688985
  quad_x(3) = -0.679409568299024
  quad_x(4) = -0.433395394129247
  quad_x(5) = -0.148874338981631
  quad_x(6) =  0.148874338981631
  quad_x(7) =  0.433395394129247
  quad_x(8) =  0.679409568299024
  quad_x(9) =  0.865063366688985
  quad_x(10) = 0.973906528517172
!
!  Weight factors
!
  quad_w(1) =  0.066671344308688
  quad_w(2) =  0.149451349150581
  quad_w(3) =  0.219086362515982
  quad_w(4) =  0.269266719309996
  quad_w(5) =  0.295524224714753
  quad_w(6) =  0.295524224714753
  quad_w(7) =  0.269266719309996
  quad_w(8) =  0.219086362515982
  quad_w(9) =  0.149451349150581
  quad_w(10) = 0.066671344308688

  return
end
function qq ( x )
!
!*******************************************************************************
!
!! QQ returns the value of the coefficient function Q(X).
!
  implicit none
!
  real qq
  real x
!
!  Test problem 1
!
!     qq = 1.0
!
!  Test problem 2
!
  qq = 0.0

  return
end
subroutine sol ( a, alpha, beta, f, np, quad_num, quad_w, quad_x )
!
!*******************************************************************************
!
!! SOL solves a linear system for the finite element coefficients.
!
  implicit none
!
  integer np
  integer quad_num
!
  real a(0:np)
  real alpha(np)
  real beta(np)
  real f(0:np)
  real ff
  integer i
  integer iq
  real phii
  real phiix
  real t
  real quad_w(quad_num)
  real x
  real quad_x(quad_num)
!
  f(0:np) = 0.0E+00

  do iq = 1, quad_num
    x = quad_x(iq)
    t = ff(x) * quad_w(iq)
    do i = 0, np
      call phi(alpha,beta,i,np,phii,phiix,x)
      f(i) = f(i) + phii * t
    end do
  end do

  f(0:np) = f(0:np) / a(0:np)

  return
end
function uex ( x )
!
!*******************************************************************************
!
!! UEX returns the value of the exact solution at a point X.
!
  implicit none
!
  real, parameter :: pi = 3.14159265E+00
  real uex
  real x
!
!  Test problem 1
!
!     uex = 1.0-x**4
!
!  Test problem 2
!
  uex = cos ( 0.5 * pi * x )

  return
end
subroutine timestamp ( )
!
!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none
!
  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone
!
  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end