chapter xii, example 4. quasi-interpolant with good knots. c from * a practical guide to splines * by c. de boor calls bsplpp(bsplvb) c integer i,irate,istep,j,l,m,mmk,n,nlow,nhigh,nm1 real algerp,aloger,bcoef(22),break(20),c(4,20),decay,dg,ddg * ,dtip1,dtip2,dx,errmax,g,h,pnatx,scrtch(4,4),step,t(26),taui c istep and step = float(istep) specify point density for error det- c ermination. data step, istep /20., 20/ c g is the function to be approximated, dg is its first, and c ddg its second derivative . g(x) = sqrt(x+1.) dg(x) = .5/g(x) ddg(x) = -.5*dg(x)/(x+1.) decay = 0. c read in the exponent irate for the knot distribution and the c lower and upper limit for the number n . read 500,irate,nlow,nhigh 500 format(3i3) print 600 600 format(28h n max.error decay exp./) c loop over n = dim( spline(4,t) ) - 2 . c n is chosen as the parameter in order to afford compar- c ison with examples 2 and 3 in which cubic spline interp- c olation at n data points was used . do 40 n=nlow,nhigh,2 nm1 = n-1 h = 1./float(nm1) m = n+2 mmk = m-4 do 5 i=1,4 t(i) = -1. 5 t(m+i) = 1. c interior knots are equidistributed with respect to the c function (x + 1)**(1/irate) . do 6 i=1,mmk 6 t(i+4) = 2.*(float(i)*h)**irate - 1. c construct quasi-interpolant. c bcoef(1) = g(-1.) = 0. bcoef(1) = 0. dtip2 = t(5) - t(4) taui = t(5) c special choice of tau(2) to avoid infinite c derivatives of g at left endpoint . bcoef(2) = g(taui) - 2.*dtip2*dg(taui)/3. * + dtip2**2*ddg(taui)/6. do 15 i=3,m taui = t(i+2) dtip1 = dtip2 dtip2 = t(i+3) - t(i+2) c formula xii(30) of text is used . 15 bcoef(i) = g(taui) + (dtip2-dtip1)*dg(taui)/3. * - dtip1*dtip2*ddg(taui)/6. c convert to pp-representation . call bsplpp(t,bcoef,m,4,scrtch,break,c,l) c estimate max.interpolation error on (-1,1). errmax = 0. c loop over cubic pieces ... do 30 i=1,l dx = (break(i+1)-break(i))/step c error is calculated at istep points per piece. do 30 j=1,istep h = float(j)*dx pnatx = c(1,i)+h*(c(2,i)+h*(c(3,i)+h*c(4,i)/3.)/2.) 30 errmax = amax1(errmax,abs(g(break(i)+h)-pnatx)) c calculate decay exponent . aloger = alog(errmax) if (n .gt. nlow) decay = * (aloger - algerp)/alog(float(n)/float(n-2)) algerp = aloger c 40 print 640,n,errmax,decay 640 format(i3,e12.4,f11.2) stop end