subroutine colloc(aleft,aright,lbegin,iorder,ntimes,addbrk,relerr)
c from * a practical guide to splines * by c. de boor
chapter xv, example. solution of an ode by collocation.
calls colpnt, difequ(ppvalu(interv)), knots, eqblok(putit(difequ*,
c bsplvd(bsplvb)))), slvblk(various subprograms), bsplpp(bsplvb*),
c newnot
c
c****** i n p u t ******
c aleft, aright endpoints of interval of approximation
c lbegin initial number of polynomial pieces in the approximation.
c a uniform breakpoint sequence is chosen.
c iorder order of polynomial pieces in the approximation
c ntimes number of passes through n e w n o t to be made
c addbrk the number (possibly fractional) of breaks to be added per
c pass through newnot. e.g., if addbrk = .33334, then a break-
c point will be added at every third pass through newnot.
c relerr a tolerance. newton iteration is stopped if the difference
c between the b-coeffs of two successive iterates is no more
c than relerr*(absol.largest b-coefficient).
c
c****** p r i n t e d o u t p u t ******
c consists of the pp-representation of the approximate solution,
c and of the error at selected points.
c
c****** m e t h o d ******
c the m-th order ordinary differential equation with m side condit-
c ions, to be specified in subroutine d i f e q u , is solved approx-
c imately by collocation.
c the approximation f to the solution g is pp of order k+m with
c l pieces and m-1 continuous derivatives. f is determined by the
c requirement that it satisfy the d.e. at k points per interval (to
c be specified in c o l p n t ) and the m side conditions.
c this usually nonlinear system of equations for f is solved by
c newton's method. the resulting linear system for the b-coeffs of an
c iterate is constructed appropriately in e q b l o k and then solved
c in s l v b l k , a program designed to solve a l m o s t b l o c k
c d i a g o n a l linear systems efficiently.
c there is an opportunity to attempt improvement of the breakpoint
c sequence (both in number and location) through use of n e w n o t .
c
integer npiece
parameter (npiece=100)
integer iorder,lbegin,ntimes, i,iflag,ii,integs(3,npiece),iside
* ,iter,itermx,k,kmax,kpm,l,lenblk,lnew,m,n,nbloks
* ,ndim,ncoef,nncoef,nt
parameter (ndim=200,kmax=20,ncoef=npiece*kmax,lenblk=ncoef)
real addbrk,aleft,aright,relerr, a(ndim),amax,asave(ndim)
* ,b(ndim),bloks(lenblk),break,coef,dx,err,rho,t(ndim)
* ,templ(lenblk),temps(ndim),xside
equivalence (bloks,templ)
common /approx/ break(npiece), coef(ncoef), l,kpm
common /side/ m, iside, xside(10)
common /other/ itermx,k,rho(kmax-1)
c
kpm = iorder
if (lbegin*kpm .gt. ncoef) go to 999
c *** set the various parameters concerning the particular dif.equ.
c including a first approx. in case the de is to be solved by
c iteration ( itermx .gt. 0) .
call difequ (1, temps(1), temps )
c *** obtain the k collocation points for the standard interval.
k = kpm - m
call colpnt ( k, rho )
c *** the following five statements could be replaced by a read in or-
c der to obtain a specific (nonuniform) spacing of the breakpnts.
dx = (aright - aleft)/float(lbegin)
temps(1) = aleft
do 4 i=2,lbegin
4 temps(i) = temps(i-1) + dx
temps(lbegin+1) = aright
c *** generate, in knots, the required knots t(1),...,t(n+kpm).
call knots ( temps, lbegin, kpm, t, n )
nt = 1
c *** generate the almost block diagonal coefficient matrix bloks and
c right side b from collocation equations and side conditions.
c then solve via slvblk , obtaining the b-representation of the ap-
c proximation in t , a , n , kpm .
10 call eqblok(t,n,kpm,temps,a,bloks,lenblk,integs,nbloks,b)
call slvblk(bloks,integs,nbloks,b,temps,a,iflag)
iter = 1
if (itermx .le. 1) go to 30
c *** save b-spline coeff. of current approx. in asave , then get new
c approx. and compare with old. if coeff. are more than relerr
c apart (relatively) or if no. of iterations is less than itermx ,
c continue iterating.
20 call bsplpp(t,a,n,kpm,templ,break,coef,l)
do 25 i=1,n
25 asave(i) = a(i)
call eqblok(t,n,kpm,temps,a,bloks,lenblk,integs,nbloks,b)
call slvblk(bloks,integs,nbloks,b,temps,a,iflag)
err = 0.
amax = 0.
do 26 i=1,n
amax = amax1(amax,abs(a(i)))
26 err = amax1(err,abs(a(i)-asave(i)))
if (err .le. relerr*amax) go to 30
iter = iter+1
if (iter .lt. itermx) go to 20
c *** iteration (if any) completed. print out approx. based on current
c breakpoint sequence, then try to improve the sequence.
30 print 630,kpm,l,n,(break(i),i=2,l)
630 format(47h approximation from a space of splines of order,i3
* ,4h on ,i3,11h intervals,/13h of dimension,i4
* ,16h. breakpoints -/(5e20.10))
if (itermx .gt. 0) print 635,iter,itermx
635 format(6h after,i3,3h of,i3,20h allowed iterations,)
call bsplpp(t,a,n,kpm,templ,break,coef,l)
print 637
637 format(46h the pp representation of the approximation is)
do 38 i=1,l
ii = (i-1)*kpm
38 print 638, break(i),(coef(ii+j),j=1,kpm)
638 format(f9.3,e13.6,10e11.3)
c *** the following call is provided here for possible further analysis
c of the approximation specific to the problem being solved.
c it is, of course, easily omitted.
call difequ ( 4, temps(1), temps )
c
if (nt .gt. ntimes) return
c *** from the pp-rep. of the current approx., obtain in newnot a new
c (and possibly better) sequence of breakpoints, adding (on the ave-
c rage) a d d b r k breakpoints per pass through newnot.
lnew = lbegin + ifix(float(nt)*addbrk)
if (lnew*kpm .gt. ncoef) go to 999
call newnot(break,coef,l,kpm,temps,lnew,templ)
call knots(temps,lnew,kpm,t,n)
nt = nt + 1
go to 10
999 nncoef = ncoef
print 699,nncoef
699 format(11h **********/23h the assigned dimension,i5
* ,25h for coef is too small.)
return
end