1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
|
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,x
* ,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
|
