colloc Subroutine

subroutine colloc ( aleft, aright, lbegin, iorder, ntimes, addbrk, relerr )

  !*****************************************************************************80
  !
  !! COLLOC solves an ordinary differential equation by collocation.
  !
  !  Method:
  !
  !    The M-th order ordinary differential equation with M side 
  !    conditions, to be specified in subroutine DIFEQU, is solved 
  !    approximately by collocation.
  !
  !    The approximation F to the solution G is piecewise polynomial of order 
  !    K+M with L pieces and M-1 continuous derivatives.   F is determined by 
  !    the requirement that it satisfy the differential equation at K points 
  !    per interval (to be specified in COLPNT ) and the M side conditions.
  !
  !    This usually nonlinear system of equations for F is solved by
  !    Newton's method. the resulting linear system for the B-coefficients of an
  !    iterate is constructed appropriately in EQBLOK and then solved
  !    in SLVBLK, a program designed to solve almost block
  !    diagonal linear systems efficiently.
  !
  !    There is an opportunity to attempt improvement of the breakpoint
  !    sequence, both in number and location, through the use of NEWNOT.
  !
  !    Printed output consists of the piecewise polynomial representation 
  !    of the approximate solution, and of the error at selected points.
  !
  !  Modified:
  !
  !    14 February 2007
  !
  !  Author:
  !
  !    Carl DeBoor
  !
  !  Reference:
  !
  !    Carl DeBoor,
  !    A Practical Guide to Splines,
  !    Springer, 2001,
  !    ISBN: 0387953663,
  !    LC: QA1.A647.v27.
  !
  !  Parameters: 
  !
  !    Input, real ( kind = 8 ) ALEFT, ARIGHT, the endpoints of the interval.
  !
  !    Input, integer ( kind = 4 ) LBEGIN, the initial number of polynomial 
  !    pieces in the approximation.  A uniform breakpoint sequence will be chosen.
  !
  !    Input, integer ( kind = 4 ) IORDER, the order of the polynomial pieces 
  !    to be used in the approximation
  !
  !    Input, integer ( kind = 4 ) NTIMES, the number of passes to be made 
  !    through NEWNOT.
  !
  !    Input, real ( kind = 8 ) ADDBRK, the number, possibly fractional, of 
  !    breaks to be added per pass through NEWNOT.  For instance, if 
  !    ADDBRK = 0.33334, then a breakpoint will be added at every third pass 
  !    through NEWNOT.
  !
  !    Input, real ( kind = 8 ) RELERR, a tolerance.  Newton iteration is 
  !    stopped if the difference between the B-coefficients of two successive 
  !    iterates is no more than RELERR*(absolute largest B-coefficient).
  !
  implicit none

  integer ( kind = 4 ), parameter :: npiece = 100
  integer ( kind = 4 ), parameter :: ndim = 200
  integer ( kind = 4 ), parameter :: ncoef = 2000
  integer ( kind = 4 ), parameter :: lenblk = 2000

  real ( kind = 8 ) a(ndim)
  real ( kind = 8 ) addbrk
  real ( kind = 8 ) aleft
  real ( kind = 8 ) amax
  real ( kind = 8 ) aright
  real ( kind = 8 ) asave(ndim)
  real ( kind = 8 ) b(ndim)
  real ( kind = 8 ) bloks(lenblk)
  real ( kind = 8 ) break
  real ( kind = 8 ) coef
  real ( kind = 8 ) dx
  real ( kind = 8 ) err
  integer ( kind = 4 ) i
  integer ( kind = 4 ) iflag
  integer ( kind = 4 ) ii
  integer ( kind = 4 ) integs(3,npiece)
  integer ( kind = 4 ) iorder
  integer ( kind = 4 ) iside
  integer ( kind = 4 ) itemps(ndim)
  integer ( kind = 4 ) iter
  integer ( kind = 4 ) itermx
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) kpm
  integer ( kind = 4 ) l
  integer ( kind = 4 ) lbegin
  integer ( kind = 4 ) lnew
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nbloks
  integer ( kind = 4 ) nt
  integer ( kind = 4 ) ntimes
  real ( kind = 8 ) relerr
  real ( kind = 8 ) rho
  real ( kind = 8 ) t(ndim)
  real ( kind = 8 ) templ(lenblk)
  real ( kind = 8 ) temps(ndim)
  real ( kind = 8 ) xside

  equivalence ( bloks, templ )

  save / approx /
  save / other /
  save / side /

  common / approx / break(npiece), coef(ncoef), l, kpm
  common / other / itermx, k, rho(19)
  common / side / m, iside, xside(10)

  kpm = iorder

  if ( ncoef < lbegin * kpm ) then
     write ( *, '(a)' ) ' '
     write ( *, '(a)' ) 'COLLOC - Fatal error!'
     write ( *, '(a)' ) '  The assigned dimension for COEF is too small.'
     stop
  end if
  !
  !  Set the various parameters concerning the particular differential
  !  equation, including a first approximation in case the differential
  !  equation is to be solved by iteration ( 0 < ITERMX ).
  !
  call difequ ( 1, temps(1), temps )
  !
  !  Obtain the K collocation points for the standard interval.
  !
  k = kpm - m
  call colpnt ( k, rho )
  !
  !  The following five statements could be replaced by a read in
  !  order to obtain a nonuniform spacing of the breakpoints.
  !
  dx = ( aright - aleft ) / real ( lbegin, kind = 8 )

  temps(1) = aleft
  do i = 2, lbegin
     temps(i) = temps(i-1) + dx
  end do
  temps(lbegin+1) = aright
  !
  !  Generate the required knots T(1:N+KPM).
  !
  call knots ( temps, lbegin, kpm, m, t, n )
  nt = 1
  !
  !  Generate the almost block diagonal coefficient matrix BLOKS and
  !  right hand side B from collocation equations and side conditions.
  !
  !  Then solve via SLVBLK, obtaining the B-representation of the 
  !  approximation in T, A, N, KPM.
  !
  do

     call eqblok ( t, n, kpm, temps, a, bloks, lenblk, integs, nbloks, b )

     call slvblk ( bloks, integs, nbloks, b, itemps, a, iflag )
     !
     !  Save B-spline coefficients of current approximation in ASAVE, then 
     !  get new approximation and compare with old. 
     !
     !  If coefficients are more than RELERR apart (relatively) or if number 
     !  of iterations is less than ITERMX, continue iterating.
     !
     do iter = 1, itermx

        call bsplpp ( t, a, n, kpm, templ, break, coef, l )

        asave(1:n) = a(1:n)

        call eqblok ( t, n, kpm, temps, a, bloks, lenblk, integs, nbloks, b )

        call slvblk ( bloks, integs, nbloks, b, itemps, a, iflag )

        amax = maxval ( abs ( a(1:n) ) )
        err = maxval ( abs ( a(1:n) - asave(1:n) ) )

        if ( err <= relerr * amax ) then
           exit
        end if

     end do
     !
     !  Iteration (if any) completed.  Print out approximation based on current
     !  breakpoint sequence, then try to improve the sequence.
     !
     write ( *, '(a)' ) ' '
     write ( *,'(a,i3,a,i3,a)' ) &
          '  Approximation from a space of splines of order ', kpm, &
          ' on ', l, ' intervals'
     write ( *, '(a,i4)' ) '  of dimension ', n
     write ( *, '(a)' ) ' '
     write ( *, '(a)' ) '  Breakpoints:'
     write ( *, '(a)' ) ' '
     write ( *, '(5g14.6)' ) break(2:l)

     if ( 0 < itermx ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a,i8)' ) '  Results on interation ', iter
     end if

     call bsplpp ( t, a, n, kpm, templ, break, coef, l )

     write ( *, '(a)' ) ' '
     write ( *, '(a)' ) &
          '  The piecewise polynomial representation of the approximation:'
     write ( *, '(a)' ) ' '

     do i = 1, l
        ii = ( i - 1 ) * kpm
        write ( *, '(f9.3,2x,e12.4,10e11.3)' ) break(i), coef(ii+1:ii+kpm)
     end do
     !
     !  The following call is provided here for possible further analysis
     !  of the approximation specific to the problem being solved.
     !  It is, of course, easily omitted.
     !
     call difequ ( 4, temps(1), temps )

     if ( ntimes < nt ) then
        exit
     end if
     !
     !  From the piecewise polynomial representation of the current approximation, 
     !  obtain in NEWNOT a new, and possibly better, sequence of breakpoints, 
     !  adding, on average, ADDBRK breakpoints per pass through NEWNOT.
     !
     lnew = lbegin + int ( real ( nt, kind = 8 ) * addbrk )

     if ( ncoef < lnew * kpm ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'COLLOC - Fatal error!'
        write ( *, '(a)' ) '  The assigned dimension for COEF is too small.'
        stop
     end if

     call newnot ( break, coef, l, kpm, temps, lnew, templ )

     call knots ( temps, lnew, kpm, m, t, n )
     nt = nt + 1

  end do

  return  
end subroutine colloc