!>
!> @file colloc.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see .
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran
!>
subroutine colloc ( aleft, aright, lbegin, iorder, ntimes, addbrk, &
relerr )
!*************************************************************************
!
!! 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 use of NEWNOT.
!
! Printed output consists of the pp-representation of the approximate
! solution, and of the error at selected points.
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! Input, real ( kind = 8 ) ALEFT, ARIGHT, the endpoints of the interval.
!
! Input, integer LBEGIN, the initial number of polynomial pieces
! in the approximation. A uniform breakpoint sequence will be chosen.
!
! Input, integer IORDER, the order of the polynomial pieces to be
! used in the approximation
!
! Input, integer NTIMES, the number of passes to be made through NEWNOT.
!
! addbrk the number (possibly fractional) of breaks to be added per
! pass through newnot. e.g., if addbrk=.33334, then a break-
! point will be added at every third pass through newnot.
!
! 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, parameter :: npiece = 100
integer, parameter :: ndim = 200
integer, parameter :: ncoef = 2000
integer, 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 i
integer iflag
integer ii
integer integs(3,npiece)
integer iorder
integer iside
integer itemps(ndim)
integer iter
integer itermx
integer j
integer k
integer kpm
integer l
integer lbegin
integer lnew
integer m
integer n
integer nbloks
integer nt
integer 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)
common /approx/ break(npiece),coef(ncoef),l,kpm
common /side/ m,iside,xside(10)
common /other/ itermx,k,rho(19)
kpm = iorder
if ( ncoef < lbegin * kpm ) then
go to 120
end if
!
! Set the various parameters concerning the particular dif.equ.
! including a first approximation in case the de 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 or-
! der to obtain a specific (nonuniform) spacing of the breakpnts.
!
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),...,t(n+kpm).
!
call knots ( temps, lbegin, kpm, t, n )
nt = 1
!
! Generate the almost block diagonal coefficient matrix bloks and
! right side b from collocation equations and side conditions.
! then solve via slvblk , obtaining the b-representation of the
! approximation in T, A, N, KPM.
!
20 continue
call eqblok ( t, n, kpm, temps, a, bloks, lenblk, integs, nbloks, b )
call slvblk ( bloks, integs, nbloks, b, itemps, a, iflag )
iter = 1
if ( itermx <= 1 ) then
go to 60
end if
!
! Save b-spline coefficients of current approx. in asave , then get new
! approx. and compare with old. if coefficients are more than relerr
! apart (relatively) or if number of iterations is less than itermx ,
! continue iterating.
!
30 continue
call bsplpp(t,a,n,kpm,templ,break,coef,l)
do i = 1, n
asave(i) = a(i)
end do
call eqblok ( t, n, kpm, temps, a, bloks, lenblk, integs, nbloks, b )
call slvblk(bloks,integs,nbloks,b,itemps,a,iflag)
err = 0.0D+00
amax = 0.0D+00
do i = 1, n
amax = max ( amax, abs ( a(i) ) )
err = max ( err, abs ( a(i)-asave(i) ) )
end do
if ( err <= relerr*amax ) then
go to 60
end if
iter = iter + 1
if ( iter < itermx ) then
go to 30
end if
!
! Iteration (if any) completed. print out approx. based on current
! breakpoint sequence, then try to improve the sequence.
!
60 continue
write(*,70)kpm,l,n,(break(i),i=2,l)
70 format (' approximation from a space of splines of order',i3, &
' on ',i3,' intervals,'/' of dimension',i4,'. breakpoints -'/ &
(5e20.10))
if ( 0 < itermx ) then
write(*,*)' '
write(*,*)'Results on interation ',iter
end if
call bsplpp(t,a,n,kpm,templ,break,coef,l)
write ( *, * ) ' '
write ( *, * ) 'The piecewise polynomial representation of the approximation:'
write ( *, * ) ' '
do i = 1, l
ii = ( i - 1 ) * kpm
write(*,'(f9.3,e13.6,10e11.3)')break(i),(coef(ii+j),j=1,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
return
end if
!
! From the pp-rep. of the current approx., obtain in NEWNOT a new
! (and possibly better) sequence of breakpoints, adding (on the
! average) ADDBRK breakpoints per pass through NEWNOT.
!
lnew = lbegin + int ( real ( nt, kind = 8 ) * addbrk )
if ( ncoef < lnew * kpm ) then
go to 120
end if
call newnot(break,coef,l,kpm,temps,lnew,templ)
call knots ( temps, lnew, kpm, t, n )
nt = nt+1
go to 20
120 continue
write(*,*)' '
write(*,*)'COLLOC - Fatal error!'
write(*,*)' The assigned dimension for COEF is ',ncoef
write(*,*)' but this is too small.'
stop
end