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
|
chapter ii. runge example
c from * a practical guide to splines * by c. de boor
integer i,istep,j,k,n,nmk,nm1
real aloger,algerp,d(20),decay,dx,errmax,g,h,pnatx,step,tau(20),x
data step, istep /20., 20/
g(x) = 1./(1.+(5.*x)**2)
print 600
600 format(28h n max.error decay exp.//)
decay = 0.
do 40 n=2,20,2
c choose interpolation points tau(1), ..., tau(n) , equally
c spaced in (-1,1), and set d(i) = g(tau(i)), i=1,...,n.
nm1 = n-1
h = 2./float(nm1)
do 10 i=1,n
tau(i) = float(i-1)*h - 1.
10 d(i) = g(tau(i))
c calculate the divided differences for the newton form.
c
do 20 k=1,nm1
nmk = n-k
do 20 i=1,nmk
20 d(i) = (d(i+1)-d(i))/(tau(i+k)-tau(i))
c
c estimate max.interpolation error on (-1,1).
errmax = 0.
do 30 i=2,n
dx = (tau(i)-tau(i-1))/step
do 30 j=1,istep
x = tau(i-1) + float(j)*dx
c evaluate interp.pol. by nested multiplication
c
pnatx = d(1)
do 29 k=2,n
29 pnatx = d(k) + (x-tau(k))*pnatx
c
30 errmax = amax1(errmax,abs(g(x)-pnatx))
aloger = alog(errmax)
if (n .gt. 2) decay =
* (aloger - algerp)/alog(float(n)/float(n-2))
algerp = aloger
40 print 640,n,errmax,decay
640 format(i3,e12.4,f11.2)
stop
end
|
