d9chm.f
1.47 KB
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
double precision function d9chm (a, c, x)
double precision a, c, x, aa(4), bb(4), ct1, ct2, eps,
1 g1, g2, g3, x2, xn0, xn1, xn2, xn3, d1mach
external d1mach
double precision bmax
c
data eps /0.0d0/
c
if (eps.eq.0.0d0) eps = 4.0d0*d1mach(4)
c
10 aa(1) = 1.0d0
bb(1) = 1.0d0
c
x2 = -0.5d0*x
bb(2) = 1.0d0 + (1.0d0+a)*x2/c
aa(2) = bb(2) + a*x/c
c
ct2 = x2/(1.0d0+c)
bb(3) = 1.0d0 + (2.0d0+bb(2))*(2.0d0+a)*ct2/3.0d0
aa(3) = bb(3) + (1.0d0+ct2)*a*x/c
bmax = dmax1 (bb(1), bb(2), bb(3))
c
do 30 i=3,300
xn0 = i
xn1 = i - 1
xn2 = i - 2
xn3 = i - 3
ct1 = 2*i - 3
c
ct2 = x2/(ct1*(c+xn1))
g1 = 1.0d0 + ct2*(xn2-a)
ct2 = ct2*(a+xn1)/(c+xn2)
g2 = ct2*(c-xn1+(a+xn0)*x2/(ct1+2.0d0))
g3 = ct2*x2*x2*(a+xn2)*(a-xn2)/(ct1*(ct1-2.0d0)*(c+xn3))
c
bb(4) = g1*bb(3) + g2*bb(2) + g3*bb(1)
aa(4) = g1*aa(3) + g2*aa(2) + g3*aa(1)
bmax = dmax1 (bb(4), bmax)
c
d9chm = aa(4)/bb(4)
if (dabs(d9chm-aa(1)/bb(1)).lt.eps*dabs(d9chm)) go to 40
c
do 20 j=1,3
aa(j) = aa(j+1)
bb(j) = bb(j+1)
20 continue
30 continue
call seteru (40hd9chm no convergence in 300 iterations, 40,
1 1, 2)
c
40 continue
errest = d1mach(4)*dmax1(dabs(d9chm), 1.0d0)
print 500, i, aa(4), bb(4), d9chm, bmax, errest
500 format (1x, i3, 3d15.5, 2e15.5)
c
return
end