grule.f 2.08 KB
 Markus Betzinger committed Apr 26, 2016 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 `````` MODULE m_grule CONTAINS SUBROUTINE grule(n,x,w) c*********************************************************************** c determines the (n+1)/2 nonnegative points x(i) and c the corresponding weights w(i) of the n-point c gauss-legendre integration rule, normalized to the c interval (-1,1). the x(i) appear in descending order. c this routine is from 'methods of numerical integration', c p.j. davis and p. rabinowitz, page 369. c m.w. c*********************************************************************** USE m_constants IMPLICIT NONE C .. C .. Arguments .. INTEGER, INTENT (IN) :: n REAL, INTENT(OUT) :: w(n/2),x(n/2) C .. C .. Locals .. INTEGER i,it,k,m REAL d1,d2pn,d3pn,d4pn,den,dp,dpn,e1,fx,h REAL p,pk,pkm1,pkp1,t,t1,u,v,x0 C .. C .. m = (n+1)/2 e1 = n* (n+1) DO i = 1,m t = (4*i-1)*pi_const/ (4*n+2) x0 = (1.- (1.-1./n)/ (8.*n*n))*cos(t) c---> iterate on the value (m.w. jan. 1982) DO it = 1,2 pkm1 = 1. pk = x0 DO k = 2,n t1 = x0*pk pkp1 = t1 - pkm1 - (t1-pkm1)/k + t1 pkm1 = pk pk = pkp1 ENDDO den = 1. - x0*x0 d1 = n* (pkm1-x0*pk) dpn = d1/den d2pn = (2.*x0*dpn-e1*pk)/den d3pn = (4.*x0*d2pn+ (2.-e1)*dpn)/den d4pn = (6.*x0*d3pn+ (6.-e1)*d2pn)/den u = pk/dpn v = d2pn/dpn h = -u* (1.+.5*u* (v+u* (v*v-u*d3pn/ (3.*dpn)))) p = pk + h* (dpn+.5*h* (d2pn+h/3.* (d3pn+.25*h*d4pn))) dp = dpn + h* (d2pn+.5*h* (d3pn+h*d4pn/3.)) h = h - p/dp x0 = x0 + h ENDDO x(i) = x0 fx = d1 - h*e1* (pk+.5*h* (dpn+h/3.* (d2pn+.25*h* (d3pn+ + .2*h*d4pn)))) w(i) = 2.* (1.-x(i)*x(i))/ (fx*fx) ENDDO IF (m+m.GT.n) x(m) = 0. RETURN END SUBROUTINE grule END MODULE M_grule``````