1
subroutine dqmomo(alfa,beta,ri,rj,rg,rh,integr)
2
c***begin prologue dqmomo
3
c***date written 820101 (yymmdd)
4
c***revision date 830518 (yymmdd)
5
c***category no. h2a2a1,c3a2
6
c***keywords modified chebyshev moments
7
c***author piessens,robert,appl. math. & progr. div. - k.u.leuven
8
c de doncker,elise,appl. math. & progr. div. - k.u.leuven
9
c***purpose this routine computes modified chebsyshev moments. the k-th
10
c modified chebyshev moment is defined as the integral over
11
c (-1,1) of w(x)*t(k,x), where t(k,x) is the chebyshev
12
c polynomial of degree k.
15
c modified chebyshev moments
16
c standard fortran subroutine
17
c double precision version
20
c alfa - double precision
21
c parameter in the weight function w(x), alfa.gt.(-1)
23
c beta - double precision
24
c parameter in the weight function w(x), beta.gt.(-1)
26
c ri - double precision
27
c vector of dimension 25
28
c ri(k) is the integral over (-1,1) of
29
c (1+x)**alfa*t(k-1,x), k = 1, ..., 25.
31
c rj - double precision
32
c vector of dimension 25
33
c rj(k) is the integral over (-1,1) of
34
c (1-x)**beta*t(k-1,x), k = 1, ..., 25.
36
c rg - double precision
37
c vector of dimension 25
38
c rg(k) is the integral over (-1,1) of
39
c (1+x)**alfa*log((1+x)/2)*t(k-1,x), k = 1, ..., 25.
41
c rh - double precision
42
c vector of dimension 25
43
c rh(k) is the integral over (-1,1) of
44
c (1-x)**beta*log((1-x)/2)*t(k-1,x), k = 1, ..., 25.
47
c input parameter indicating the modified
48
c moments to be computed
49
c integr = 1 compute ri, rj
50
c = 2 compute ri, rj, rg
51
c = 3 compute ri, rj, rh
52
c = 4 compute ri, rj, rg, rh
55
c***routines called (none)
56
c***end prologue dqmomo
58
double precision alfa,alfp1,alfp2,an,anm1,beta,betp1,betp2,ralf,
62
dimension rg(25),rh(25),ri(25),rj(25)
65
c***first executable statement dqmomo
73
c compute ri, rj using a forward recurrence relation.
77
ri(2) = ri(1)*alfa/alfp2
78
rj(2) = rj(1)*beta/betp2
82
ri(i) = -(ralf+an*(an-alfp2)*ri(i-1))/(anm1*(an+alfp1))
83
rj(i) = -(rbet+an*(an-betp2)*rj(i-1))/(anm1*(an+betp1))
87
if(integr.eq.1) go to 70
88
if(integr.eq.3) go to 40
90
c compute rg using a forward recurrence relation.
93
rg(2) = -(ralf+ralf)/(alfp2*alfp2)-rg(1)
98
rg(i) = -(an*(an-alfp2)*rg(im1)-an*ri(im1)+anm1*ri(i))/
104
if(integr.eq.2) go to 70
106
c compute rh using a forward recurrence relation.
108
40 rh(1) = -rj(1)/betp1
109
rh(2) = -(rbet+rbet)/(betp2*betp2)-rh(1)
114
rh(i) = -(an*(an-betp2)*rh(im1)-an*rj(im1)+
115
* anm1*rj(i))/(anm1*(an+betp1))