chapter xiii, example 1. a large norm amplifies noise. c from * a practical guide to splines * by c. de boor calls splint(banfac/slv,bsplvb),bsplpp(bsplvb*),ppvalu(interv),round c an initially uniform data point distribution of n points is c changed i t e r m x times by moving the j c l o s e -th data point c toward its left neighbor, cutting the distance between the two by a c factor of r a t e each time . to demonstrate the corresponding c increase in the norm of cubic spline interpolation at these data c points, the data are taken from a cubic polynomial (which would be c interpolated exactly) but with noise of size s i z e added. the re- c sulting noise or error in the interpolant (compared to the cubic) c gives the norm or noise amplification factor and is printed out c together with the diminishing distance h between the two data c points. integer i,iflag,istep,iter,itermx,j,jclose,l,n,nmax,nm1 parameter (nmax=200) real amax,bcoef(nmax),break(nmax),coef(nmax*4),dx,fltnm1,fx * ,gtau(nmax),h,rate,scrtch(nmax*7),size,step,t(nmax+4),tau(nmax),x C real amax,bcoef(200),break(200),coef(800),dx,fltnm1,fx C * ,gtau(200),h,rate,scrtch(1400),size,step,t(204),tau(200),x common /rount/ size data step, istep / 20., 20 / c function to be interpolated . g(x) = 1.+x*(1.+x*(1.+x)) read 500,n,itermx,jclose,size,rate 500 format(3i3/e10.3/e10.3) print 600,size 600 format(16h size of noise =,e8.2// * 25h h max.error) c start with uniform data points . nm1 = n - 1 fltnm1 = float(nm1) do 10 i=1,n 10 tau(i) = float(i-1)/fltnm1 c set up knot sequence for not-a-knot end condition . do 11 i=1,4 t(i) = tau(1) 11 t(n+i) = tau(n) do 12 i=5,n 12 t(i) = tau(i-2) c do 100 iter=1,itermx do 21 i=1,n 21 gtau(i) = round(g(tau(i))) call splint ( tau, gtau, t, n, 4, scrtch, bcoef, iflag ) go to (24,23),iflag 23 print 623 623 format(27h something wrong in splint.) stop 24 call bsplpp ( t, bcoef, n, 4, scrtch, break, coef, l ) c calculate max.interpolation error . amax = 0. do 30 i=4,n dx = (break(i-2) - break(i-3))/step do 25 j=2,istep x = break(i-2) - dx*float(j-1) fx = ppvalu(break,coef,l,4,x,0) 25 amax = amax1(amax,abs(fx-g(x))) 30 continue h = tau(jclose) - tau(jclose-1) print 630,h,amax 630 format(e9.2,e15.3) c move tau(jclose) toward its left neighbor so as to cut c their distance by a factor of rate . tau(jclose) = (tau(jclose) + (rate-1.)*tau(jclose-1))/rate 100 continue stop end