subroutine myquad(n,type,alpha,beta,g,w,valid) implicit none include 'max.inc' c c Purpose: to perform the computation of quadrature abscissas and weights c c Inputs: c + n: quadrature order c + type: quadrature type c 1: Legendre quadrature c 2: Legendre-Radau quadrature c 3: Legendre-Lobatto quadrature c 4: Jacobi quadrature c + alpha,beta: values used in the Jacobi polynomials c c Outputs: c + g: quadrature abscissas c + w: quadrature weights c + valid: 1 if valid solution was found; 0 otherwise c c I/O integer n,type,valid double precision alpha,beta double precision g(1:Nqmx) double precision w(1:Nqmx) c temp integer rootf double precision roots(1:Nqmx) integer n0,n1 double precision dif1(1:Nqmx) double precision dif2(1:Nqmx) double precision dif3(1:Nqmx) c label integer strlen character*(Nchar_mx) label label='subroutine myquad' valid=1 if ((type.ge.1).and.(type.le.3)) then ! Legendre quadrature call findroots(n,type,roots,rootf) else if (type.eq.4) then ! jacobi quadrature n0=0 n1=0 call jacobi_roots(Nqmx,n,n0,n1,alpha,beta,dif1,dif2,dif3,roots) rootf=1 endif if (rootf.eq.0) then valid=0 goto 111 endif if (type.eq.1) then call gauss_legendre(n,roots,g,w) else if (type.eq.2) then call gauss_legendre_radau(n,roots,g,w) else if (type.eq.3) then call gauss_legendre_lobatto(n,roots,g,w) else if (type.eq.4) then call jacobi_weights(n,alpha,beta,roots,g,w) endif 111 continue return end double precision function dk_jacobi(n,k,alpha,beta,x) implicit none c c Purpose: to evaluate the k-th derivative of the specified Jacobi polynomial c c Inputs: c + n: order of the polynomial c + k: order of the derivative c + alpha, beta: alpha and beta in the definition of Pn(alpha,beta)(x) c + x: value of x c c Outputs: c + dk_jacobi: d**k(Pn(alpha,beta)(x))/dx**k c c I/O integer n,k double precision alpha,beta,x c temp double precision dgammaf double precision jacobi dk_jacobi=dgammaf(alpha+beta+n+1.0D+0+k) & /2.0D+0**k & /dgammaf(alpha+beta+n+1.0D+0) & *jacobi(n-k,alpha+k,beta+k,x) return end double precision function jacobi(n,alpha,beta,x) implicit none c c Purpose: to compute the value of the specified Jacobi polynomial c c Inputs: c + n: order of the polynomial c + alpha, beta: alpha and beta in the definition of Pn(alpha,beta)(x) c + x: value of x c c Outputs: c + jacobi: Pn(alpha,beta)(x) c c I/O integer n double precision alpha,beta,x c temp double precision jcb,jcb_rec integer order,max_order double precision Pnm1,Pn,Pnp1 double precision a,b,c,d max_order=12 if (n.le.max_order) then jacobi=jcb(n,alpha,beta,x) else Pn=jcb(max_order,alpha,beta,x) Pnm1=jcb(max_order-1,alpha,beta,x) do order=max_order,n-1 Pnp1=jcb_rec(order,alpha,beta,Pn,Pnm1,x) Pnm1=Pn Pn=Pnp1 enddo jacobi=Pnp1 endif return end double precision function jcb(n,alpha,beta,x) implicit none c c Purpose: to compute the value of the specified Jacobi polynomial c using direct analytical function c c Inputs: c + n: order of the polynomial c + alpha, beta: alpha and beta in the definition of Pn(alpha,beta)(x) c + x: value of x c c Outputs: c + jcb: Pn(alpha,beta)(x) c c I/O integer n double precision alpha,beta,x c temp double precision binomial_coefficient double precision dgammaf double precision elem,sum integer fact integer m sum=0.0D+0 do m=0,n elem=binomial_coefficient(dble(n),dble(m)) & *dgammaf(alpha+beta+n+m+1.0D+0) & /dgammaf(alpha+m+1.0D+0) & *((x-1.0D+0)/2.0D+0)**m sum=sum+elem enddo jcb=sum & *dgammaf(alpha+n+1.0D+0) & /dgammaf(alpha+beta+n+1.0D+0) & /fact(n) return end double precision function jcb_rec(n,alpha,beta,Pn,Pnm1,x) implicit none c c Purpose: to compute the value of the specified Jacobi polynomial c using recurrence property c c Inputs: c + n: order of the polynomial c + alpha, beta: alpha and beta in the definition of Pn(alpha,beta)(x) c + x: value of x c + Pn: Pn(alpha,beta)(x) c + Pnm1: P(n-1)(alpha,beta)(x) c c Outputs: c + jcb_rec: P(n+1)(alpha,beta)(x) c c I/O integer n double precision alpha,beta,Pn,Pnm1,x c temp double precision pochhammer jcb_rec=( & ((2*n+alpha+beta+1.0D+0)*(alpha**2.0D+0-beta**2.0D+0) & +pochhammer(2*n+alpha+beta,3)*x)*Pn & -2.0D+0*(n+alpha)*(n+beta)*(2*n+alpha+beta+2.0D+0)*Pnm1 & ) & /(2.0D+0*(n+1)*(n+alpha+beta+1.0D+0)*(2*n+alpha+beta)) return end double precision function pochhammer(x,n) implicit none c c Purpose: to compute (x)n=x*(x+1)*...*(x+n-1) c double precision x integer n integer i double precision tmp tmp=1.0D+0 do i=0,n-1 tmp=tmp*(x+i) enddo pochhammer=tmp return end double precision function binomial_coefficient(x,y) implicit none double precision x,y double precision dgammaf binomial_coefficient=dgammaf(x+1.0D+0)/ & ( & dgammaf(y+1.0D+0)* & dgammaf(x-y+1.0D+0) & ) return end **************************************************************** * * The following routines (JCOBI, DIF, DFOPR, INTRP, AND RADAU) * are the same as found in Villadsen, J. and M.L. Michelsen, * Solution of Differential Equation Models by Polynomial * Approximation, Prentice-Hall (1978) pages 418-420. * * Cosmetic changes (elimination of arithmetic IF statements, most * GO TO statements, and indentation of program blocks) made by: * * John W. Eaton * Department of Chemical Engineering * The University of Texas at Austin * Austin, Texas 78712 * * June 6, 1987 * * Some error checking additions also made on June 7, 1987 * * Further cosmetic changes made August 20, 1987 * ************************************************************************ * subroutine jacobi_roots(nd,n,n0,n1,alpha,beta, & dif1,dif2,dif3,roots) implicit none C INTEGER + + ND, N, N0, N1 C DOUBLE PRECISION + + ALPHA, BETA, DIF1(ND), DIF2(ND), DIF3(ND), ROOT(ND) double precision roots(nd) C C*********************************************************************** C C VILLADSEN AND MICHELSEN, PAGES 131-132, 418 C C THIS SUBROUTINE COMPUTES THE ZEROS OF THE JACOBI POLYNOMIAL C C (ALPHA,BETA) C P (X) C N C C USE DIF (GIVEN BELOW) TO COMPUTE THE DERIVATIVES OF THE NODE C POLYNOMIAL C C N0 (ALPHA,BETA) N1 C P (X) = (X) * P (X) * (1 - X) C NT N C C AT THE INTERPOLATION POINTS. C C INPUT PARAMETERS: C C ND : THE DIMENSION OF THE VECTORS DIF1, DIF2, DIF3, AND ROOT C C N : THE DEGREE OF THE JACOBI POLYNOMIAL, (i.e. THE NUMBER C OF INTERIOR INTERPOLATION POINTS) C C N0 : DETERMINES WHETHER X = 0 IS INCLUDED AS AN C INTERPOLATION POINT C C N0 = 0 ==> X = 0 IS NOT INCLUDED C N0 = 1 ==> X = 0 IS INCLUDED C C N1 : DETERMINES WHETHER X = 1 IS INCLUDED AS AN C INTERPOLATION POINT C C N1 = 0 ==> X = 1 IS NOT INCLUDED C N1 = 1 ==> X = 1 IS INCLUDED C C ALPHA : THE VALUE OF ALPHA IN THE DESCRIPTION OF THE JACOBI C POLYNOMIAL C C BETA : THE VALUE OF BETA IN THE DESCRIPTION OF THE JACOBI C POLYNOMIAL C C FOR A MORE COMPLETE EXPLANATION OF ALPHA AN BETA, SEE VILLADSEN C AND MICHELSEN, PAGES 57 TO 59 C C OUTPUT PARAMETERS: C C ROOTS : ONE DIMENSIONAL VECTOR CONTAINING ON EXIT THE C N + N0 + N1 ZEROS OF THE NODE POLYNOMIAL USED IN THE C INTERPOLATION ROUTINE C C DIF1 : ONE DIMENSIONAL VECTOR CONTAINING THE FIRST DERIVATIVE C OF THE NODE POLYNOMIAL AT THE ZEROS C C DIF2 : ONE DIMENSIONAL VECTOR CONTAINING THE SECOND DERIVATIVE C OF THE NODE POLYNOMIAL AT THE ZEROS C C DIF3 : ONE DIMENSIONAL VECTOR CONTAINING THE THIRD DERIVATIVE C OF THE NODE POLYNOMIAL AT THE ZEROS C C COMMON BLOCKS: NONE C C REQUIRED ROUTINES: VILERR, DIF C C*********************************************************************** C INTEGER I,J,NT,IER DOUBLE PRECISION AB,AD,AP,Z1,Z,Y,X,XD,XN,XD1,XN1,XP,XP1,ZC DOUBLE PRECISION ZERO,ONE,TWO LOGICAL LSTOP C PARAMETER ( ZERO = 0.0D+00, ONE = 1.0D+00, TWO = 2.0D+00 ) C C -- ERROR CHECKING C IF ((N0 .NE. 0) .AND. (N0 .NE. 1)) THEN IER = 1 LSTOP = .TRUE. CALL VILERR(IER,LSTOP) ELSE END IF C IF ((N1 .NE. 0) .AND. (N1 .NE. 1)) THEN IER = 2 LSTOP = .TRUE. CALL VILERR(IER,LSTOP) ELSE END IF C IF (ND .LT. (N + N0 + N1)) THEN IER = 3 LSTOP = .TRUE. CALL VILERR(IER,LSTOP) ELSE END IF C IF ((N + N0 + N1) .LT. 1) THEN IER = 7 LSTOP = .TRUE. CALL VILERR(IER,LSTOP) ELSE END IF C C -- FIRST EVALUATION OF COEFFICIENTS IN RECURSION FORMULAS. C -- RECURSION COEFFICIENTS ARE STORED IN DIF1 AND DIF2. C AB = ALPHA + BETA AD = BETA - ALPHA AP = BETA*ALPHA DIF1(1) = (AD/(AB + TWO) + ONE)/TWO DIF2(1) = ZERO C IF(N .GE. 2) THEN DO 10 I = 2,N C Z1 = DBLE(I) - ONE Z = AB + 2*Z1 DIF1(I) = (AB*AD/Z/(Z + TWO) + ONE)/TWO C IF (I .EQ. 2) THEN DIF2(I) = (AB + AP + Z1)/Z/Z/(Z + ONE) ELSE Z = Z*Z Y = Z1*(AB + Z1) Y = Y*(AP + Y) DIF2(I) = Y/Z/(Z - ONE) END IF C 10 CONTINUE ELSE END IF C C -- ROOT DETERMINATION BY NEWTON METHOD WITH SUPPRESSION OF C -- PREVIOUSLY DETERMINED ROOTS C X = ZERO C DO 20 I = 1,N C 25 CONTINUE XD = ZERO XN = ONE XD1 = ZERO XN1 = ZERO C DO 30 J = 1,N XP = (DIF1(J) - X)*XN - DIF2(J)*XD XP1 = (DIF1(J) - X)*XN1 - DIF2(J)*XD1 - XN XD = XN XD1 = XN1 XN = XP XN1 = XP1 30 CONTINUE C ZC = ONE Z = XN/XN1 C IF (I .NE. 1) THEN DO 22 J = 2,I ZC = ZC - Z/(X - ROOT(J-1)) 22 CONTINUE ELSE END IF C Z = Z/ZC X = X - Z C IF (DABS(Z) .GT. 1.D-09) THEN C C -- BACKWARD BRANCH C GO TO 25 ELSE END IF C ROOT(I) = X X = X + 0.0001D0 C 20 CONTINUE C C -- ADD INTERPOLATION POINTS AT X = 0 AND/OR X = 1 C NT = N + N0 + N1 C IF (N0 .NE. 0) THEN DO 31 I = 1,N J = N + 1 - I ROOT(J+1) = ROOT(J) 31 CONTINUE ROOT(1) = ZERO ELSE END IF C IF (N1 .EQ. 1) THEN ROOT(NT) = ONE ELSE END IF C CALL DIF ( NT, ROOT, DIF1, DIF2, DIF3 ) do i=1,n roots(i)=2.0D+0*root(i)-1.0D+0 enddo C RETURN END SUBROUTINE DIF ( NT, ROOT, DIF1, DIF2, DIF3 ) implicit none C INTEGER NT DOUBLE PRECISION ROOT(NT), DIF1(NT), DIF2(NT), DIF3(NT) C C*********************************************************************** C C SUBROUTINE DIF C C THIS ROUTINE IS NOT GIVEN SEPARATELY BY VILLADSEN AND MICHELSEN C BUT AS PART OF JCOBI C C DIF COMPUTES THE FIRST THREE DERIVATIVES OF THE NODE POLYNOMIAL C C N0 (ALPHA,BETA) N1 C P (X) = (X) * P (X) * (1 - X) C NT N C C AT THE INTERPOLATION POINTS. EACH OF THE PARAMETERS N0 AND N1 C MAY BE GIVEN THE VALUE 0 OR 1. NT = N + N0 + N1 C C THE VALUES OF ROOT MUST BE KNOWN BEFORE A CALL TO DIF IS POSSIBLE. C THEY MAY BE COMPUTED USING JCOBI. C C PARAMETER LIST: SEE THE SUBROUTINE JCOBI C C COMMON BLOCKS: NONE C C REQUIRED ROUTINES: VILERR C C*********************************************************************** C INTEGER I,J,IER DOUBLE PRECISION X,Y DOUBLE PRECISION ZERO,ONE,TWO,THREE LOGICAL LSTOP C PARAMETER ( ZERO = 0.0D+00, ONE = 1.0D+00, + TWO = 2.0D+00, THREE = 3.0D+00 ) C C -- ERROR CHECKING C IF (NT .LT. 1) THEN IER = 7 LSTOP = .TRUE. CALL VILERR(IER,LSTOP) ELSE END IF C C -- EVALUATE DERIVATIVES OF NODE POLYNOMIAL USING RECURSION FORMULAS C DO 40 I = 1,NT C X = ROOT(I) DIF1(I) = ONE DIF2(I) = ZERO DIF3(I) = ZERO C DO 30 J = 1,NT C IF (J .NE. I) THEN Y = X - ROOT(J) DIF3(I) = Y*DIF3(I) + THREE*DIF2(I) DIF2(I) = Y*DIF2(I) + TWO *DIF1(I) DIF1(I) = Y*DIF1(I) ELSE END IF C 30 CONTINUE 40 CONTINUE C RETURN END SUBROUTINE VILERR ( IER, LSTOP ) implicit none C INTEGER IER LOGICAL LSTOP C C*********************************************************************** C C THIS SUBROUTINE HANDLES ERRORS FOR THE SUBROUTINES JCOBI, DFOPR, C INTRP, AND RADAU GIVEN BY VILLADSEN AND MICHELSEN. C C PARAMETER LIST: C C IER : ERROR NUMBER C LSTOP : LOGICAL FLAG C C LSTOP = .TRUE. ==> FATAL ERROR, PROGRAM TERMINATION C LSTOP = .FALSE. ==> WARNING ERROR, NORMAL RETURN C C COMMON BLOCKS: NONE C C REQUIRED ROUTINES: NONE C C*********************************************************************** C C -- BEGIN C IF ( IER .EQ. 1) THEN C WRITE(*,*) '** VILERR : Illegal value for N0 ' C ELSE IF ( IER .EQ. 2) THEN C WRITE(*,*) '** VILERR : Illegal value for N1 ' C ELSE IF ( IER .EQ. 3 ) THEN C WRITE(*,*) '** VILERR : Insufficient dimension for problem ' C ELSE IF ( IER .EQ. 4 ) THEN C WRITE(*,*) '** VILERR : Index less than zero in DFOPR ' C ELSE IF ( IER .EQ. 5 ) THEN C WRITE(*,*) '** VILERR : Index greater than NTOTAL in DFOPR ' C ELSE IF ( IER .EQ. 6 ) THEN C WRITE(*,*) '** VILERR : Illegal ID in DFOPR ' C ELSE IF ( IER .EQ. 7 ) THEN C WRITE(*,*) '** VILERR : Number of interpolation points ' WRITE(*,*) ' less than 1 ' C ELSE IF ( IER .EQ. 8 ) THEN C WRITE(*,*) '** VILERR : Illegal ID in RADAU ' C ELSE IF ( IER .EQ. 9 ) THEN C WRITE(*,*) '** VILERR : ID = 1 but N1 not equal to 1 in RADAU ' C ELSE IF ( IER .EQ. 10 ) THEN C WRITE(*,*) '** VILERR : ID = 2 but N0 not equal to 1 in RADAU ' C ELSE IF ( IER .EQ. 11 ) THEN C WRITE(*,*) '** VILERR : ID = 3 but N0 not equal to 1 or ' WRITE(*,*) ' N1 not equal to 1 in RADAU ' C ELSE C WRITE(*,*) 'UNRECOGNIZED ERROR FLAG SET FOR VILERR ' C END IF C IF ( LSTOP ) THEN C C -- PROGRAM EXECUTION TERMINATES HERE C c CALL XSTOPX (' ') write(*,*) 'Error from vilerr' stop C ELSE END IF C RETURN END subroutine jacobi_weights(n,alpha,beta,roots,g,w) implicit none include 'max.inc' c c Purpose: to compute the quadrature abscissas and weights that are associated to c a Jacobi quadrature c c Inputs: c + n: quadrature order c + alpha,beta: coefficients for Jacobi polynomial c + roots: array of the n roots c c Outputs: c + g: array of the n quadrature abscissas c + w: array of the n quadrature weights c c I/O integer n double precision alpha,beta double precision roots(1:Nqmx) double precision g(1:Nqmx) double precision w(1:Nqmx) c temp double precision dgammaf double precision An,sgamma double precision jacobi,dk_jacobi double precision sum_w,sumt integer fact integer i c label integer strlen character*(Nchar_mx) label label='subroutine jacobi_weights' sum_w=0.0D+0 do i=1,n g(i)=(roots(i)+1.0D+0)/2.0D+0 w(i)=An(n,alpha,beta)/An(n-1,alpha,beta) & *sgamma(n-1,alpha,beta) & /( & jacobi(n-1,alpha,beta,roots(i)) & *dk_jacobi(n,1,alpha,beta,roots(i)) & ) c w(i)=-(2*n+alpha+beta+2.0D+0) c & /(n+alpha+beta+1.0D+0) c & *dgammaf(n+alpha+1.0D+0) c & *dgammaf(n+beta+1.0D+0) c & /dgammaf(n+alpha+beta+1.0D+0) c & /fact(n+1) c & *2.0D+0**(alpha+beta) c & /dk_jacobi(n,1,alpha,beta,roots(i)) c & /jacobi(n+1,alpha,beta,roots(i)) sum_w=sum_w+w(i) enddo c Normalization is performed here since the sum of weights is unknown do i=1,n w(i)=w(i)/sum_w enddo return end double precision function An(n,alpha,beta) implicit none cc Purpose: to compute coefficient An required for Jacobi quadrature weights c c Inputs: c + n: order of the Jacobi polynomial c + alpha,beta: coefficients of the Jacobi polynomial c c Ouputs: c + An c c I/O integer n double precision alpha,beta c temp integer fact double precision dgammaf An=dgammaf(2*n+alpha+beta+1.0D+0)/ & ( & 2.0D+0**n & *fact(n) & *dgammaf(n+alpha+beta+1.0D+0) & ) return end double precision function sgamma(n,alpha,beta) implicit none integer n double precision alpha,beta integer fact double precision dgammaf sgamma=2.0D+0**(2*n+alpha+beta+1.0D+0) & *fact(n) & *dgammaf(n+alpha+1.0D+0) & *dgammaf(n+beta+1.0D+0) & /( & 2.0D+0**(2*n) & *(fact(n))**2.0D+0 & *(2*n+alpha+beta+1.0D+0) & *dgammaf(n+alpha+beta+1.0D+0) & ) return end DOUBLE PRECISION FUNCTION DGAMMAF(X) implicit none C---------------------------------------------------------------------- C C This routine calculates the GAMMA function for a real argument X. C Computation is based on an algorithm outlined in reference 1. C The program uses rational functions that approximate the GAMMA C function to at least 20 significant decimal digits. Coefficients C for the approximation over the interval (1,2) are unpublished. C Those for the approximation for X .GE. 12 are from reference 2. C The accuracy achieved depends on the arithmetic system, the C compiler, the intrinsic functions, and proper selection of the C machine-dependent constants. C C C******************************************************************* C******************************************************************* C C Explanation of machine-dependent constants C C beta - radix for the floating-point representation C maxexp - the smallest positive power of beta that overflows C XBIG - the largest argument for which GAMMA(X) is representable C in the machine, i.e., the solution to the equation C GAMMA(XBIG) = beta**maxexp C XINF - the largest machine representable floating-point number; C approximately beta**maxexp C EPS - the smallest positive floating-point number such that C 1.0+EPS .GT. 1.0 C XMININ - the smallest positive floating-point number such that C 1/XMININ is machine representable C C Approximate values for some important machines are: C C beta maxexp XBIG C C CRAY-1 (S.P.) 2 8191 966.961 C Cyber 180/855 C under NOS (S.P.) 2 1070 177.803 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 2 128 35.040 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 2 1024 171.624 C IBM 3033 (D.P.) 16 63 57.574 C VAX D-Format (D.P.) 2 127 34.844 C VAX G-Format (D.P.) 2 1023 171.489 C C XINF EPS XMININ C C CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466 C Cyber 180/855 C under NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 3.40E+38 1.19E-7 1.18E-38 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 1.79D+308 2.22D-16 2.23D-308 C IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76 C VAX D-Format (D.P.) 1.70D+38 1.39D-17 5.88D-39 C VAX G-Format (D.P.) 8.98D+307 1.11D-16 1.12D-308 C C******************************************************************* C******************************************************************* C C Error returns C C The program returns the value XINF for singularities or C when overflow would occur. The computation is believed C to be free of underflow and overflow. C C C Intrinsic functions required are: C C INT, DBLE, EXP, LOG, REAL, SIN C C C References: "An Overview of Software Development for Special C Functions", W. J. Cody, Lecture Notes in Mathematics, C 506, Numerical Analysis Dundee, 1975, G. A. Watson C (ed.), Springer Verlag, Berlin, 1976. C C Computer Approximations, Hart, Et. Al., Wiley and C sons, New York, 1968. C C Latest modification: October 12, 1989 C C Authors: W. J. Cody and L. Stoltz C Applied Mathematics Division C Argonne National Laboratory C Argonne, IL 60439 C C---------------------------------------------------------------------- INTEGER I,N LOGICAL PARITY DOUBLE PRECISION 1 C,CONV,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI,SUM,TWELVE, 2 TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO DIMENSION C(7),P(8),Q(8) C---------------------------------------------------------------------- C Mathematical constants C---------------------------------------------------------------------- DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/, 1 SQRTPI/0.9189385332046727417803297D0/, 2 PI/3.1415926535897932384626434D0/ C---------------------------------------------------------------------- C Machine dependent parameters C---------------------------------------------------------------------- DATA XBIG,XMININ,EPS/171.624D0,2.23D-308,2.22D-16/, 1 XINF/1.79D308/ C---------------------------------------------------------------------- C Numerator and denominator coefficients for rational minimax C approximation over (1,2). C---------------------------------------------------------------------- DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1, 1 -3.79804256470945635097577D+2,6.29331155312818442661052D+2, 2 8.66966202790413211295064D+2,-3.14512729688483675254357D+4, 3 -3.61444134186911729807069D+4,6.64561438202405440627855D+4/ DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2, 1 -1.01515636749021914166146D+3,-3.10777167157231109440444D+3, 2 2.25381184209801510330112D+4,4.75584627752788110767815D+3, 3 -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/ C---------------------------------------------------------------------- C Coefficients for minimax approximation over (12, INF). C---------------------------------------------------------------------- DATA C/-1.910444077728D-03,8.4171387781295D-04, 1 -5.952379913043012D-04,7.93650793500350248D-04, 2 -2.777777777777681622553D-03,8.333333333333333331554247D-02, 3 5.7083835261D-03/ C---------------------------------------------------------------------- C Statement functions for conversion between integer and float C---------------------------------------------------------------------- CONV(I) = DBLE(I) PARITY = .FALSE. FACT = ONE N = 0 Y = X IF (Y .LE. ZERO) THEN C---------------------------------------------------------------------- C Argument is negative C---------------------------------------------------------------------- Y = -X Y1 = AINT(Y) RES = Y - Y1 IF (RES .NE. ZERO) THEN IF (Y1 .NE. AINT(Y1*HALF)*TWO) PARITY = .TRUE. FACT = -PI / SIN(PI*RES) Y = Y + ONE ELSE RES = XINF GO TO 900 END IF END IF C---------------------------------------------------------------------- C Argument is positive C---------------------------------------------------------------------- IF (Y .LT. EPS) THEN C---------------------------------------------------------------------- C Argument .LT. EPS C---------------------------------------------------------------------- IF (Y .GE. XMININ) THEN RES = ONE / Y ELSE RES = XINF GO TO 900 END IF ELSE IF (Y .LT. TWELVE) THEN Y1 = Y IF (Y .LT. ONE) THEN C---------------------------------------------------------------------- C 0.0 .LT. argument .LT. 1.0 C---------------------------------------------------------------------- Z = Y Y = Y + ONE ELSE C---------------------------------------------------------------------- C 1.0 .LT. argument .LT. 12.0, reduce argument if necessary C---------------------------------------------------------------------- N = INT(Y) - 1 Y = Y - CONV(N) Z = Y - ONE END IF C---------------------------------------------------------------------- C Evaluate approximation for 1.0 .LT. argument .LT. 2.0 C---------------------------------------------------------------------- XNUM = ZERO XDEN = ONE DO 260 I = 1, 8 XNUM = (XNUM + P(I)) * Z XDEN = XDEN * Z + Q(I) 260 CONTINUE RES = XNUM / XDEN + ONE IF (Y1 .LT. Y) THEN C---------------------------------------------------------------------- C Adjust result for case 0.0 .LT. argument .LT. 1.0 C---------------------------------------------------------------------- RES = RES / Y1 ELSE IF (Y1 .GT. Y) THEN C---------------------------------------------------------------------- C Adjust result for case 2.0 .LT. argument .LT. 12.0 C---------------------------------------------------------------------- DO 290 I = 1, N RES = RES * Y Y = Y + ONE 290 CONTINUE END IF ELSE C---------------------------------------------------------------------- C Evaluate for argument .GE. 12.0, C---------------------------------------------------------------------- IF (Y .LE. XBIG) THEN YSQ = Y * Y SUM = C(7) DO 350 I = 1, 6 SUM = SUM / YSQ + C(I) 350 CONTINUE SUM = SUM/Y - Y + SQRTPI SUM = SUM + (Y-HALF)*LOG(Y) RES = EXP(SUM) ELSE RES = XINF GO TO 900 END IF END IF C---------------------------------------------------------------------- C Final adjustments and return C---------------------------------------------------------------------- IF (PARITY) RES = -RES IF (FACT .NE. ONE) RES = FACT / RES 900 DGAMMAF = RES RETURN C ---------- Last line of GAMMA ---------- END subroutine gauss_legendre(n,roots,g,w) implicit none include 'max.inc' c c Purpose: to find the n quadrature weights w and abscissas g c for a n-order Legendre quadrature c c Inputs: c + n: order of the quadrature c + roots: array of the n roots c c Outputs: c + g: array of the n quadrature abscissas c + w: array of the n quadrature weights c c I/O integer n double precision roots(1:Nqmx) double precision g(1:Nqmx) double precision w(1:Nqmx) c temp integer i double precision legendre c label integer strlen character*(Nchar_mx) label label='subroutine gauss_legendre' do i=1,n g(i)=(roots(i)+1.0D+0)/2.0D+0 w(i)=(1.0D+0-roots(i)**2.0D+0)/ & ( & (n+1)**2.0D+0* & legendre(n+1,roots(i))**2.0D+0 & ) enddo 666 continue return end subroutine gauss_legendre_radau(n,roots,g,w) implicit none include 'max.inc' c c Purpose: to find the n quadrature weights w and abscissas g c for a n-order Legendre-Radau quadrature c c Inputs: c + n: order of the quadrature c + roots: array of the n roots c c Outputs: c + g: array of the n quadrature abscissas c + w: array of the n quadrature weights c + valid: 1 if a valid solution was found; 0 otherwise c c I/O integer n double precision roots(1:Nqmx) double precision g(1:Nqmx) double precision w(1:Nqmx) c temp integer i double precision legendre double precision legendre_derivative c label integer strlen character*(Nchar_mx) label label='subroutine gauss_legendre_radau' c first abscissa g(1)=0.0D+0 w(1)=1.0D+0/(n**2.0D+0) c n-1 remaining abscissas do i=1,n-1 g(i+1)=(roots(i)+1.0D+0)/2.0D+0 w(i+1)=1.0D+0/ & ( & 2.0D+0*(1.0D+0-roots(i))* & (legendre_derivative(n-1,roots(i))**2.0D+0) & ) enddo 666 continue return end subroutine gauss_legendre_lobatto(n,roots,g,w) implicit none include 'max.inc' c c Purpose: to find the n quadrature weights w and abscissas g c for a n-order Legendre-Lobatto quadrature c c Inputs: c + n: order of the quadrature c + roots: array of the n roots c c Outputs: c + g: array of the n quadrature abscissas c + w: array of the n quadrature weights c c I/O integer n double precision roots(1:Nqmx) double precision g(1:Nqmx) double precision w(1:Nqmx) c temp integer i double precision legendre c label integer strlen character*(Nchar_mx) label label='subroutine gauss_legendre_lobatto' c first abscissa g(1)=0.0D+0 w(1)=1.0D+0/(n*(n-1)) c n-2 following abscissa do i=1,n-2 g(i+1)=(roots(i)+1.0D+0)/2.0D+0 w(i+1)=1.0D+0/ & ( & n*(n-1)*(legendre(n-1,roots(i))**2.0D+0) & ) enddo c last abscissa g(n)=1.0D+0 w(n)=1.0D+0/(n*(n-1)) 666 continue return end subroutine findroots(n,type,roots,rootf) implicit none include 'max.inc' c c Purpose: to find the n roots of the n-th Legendre polynomia c c Inputs: c + n: order of the Legendre polynomia c + type: type of quadrature (1: Legendre, 2: Legendre-Radau, 3: Legendre-Lobatto) c c Outputs: c + roots: array of the n roots c + rootf: if 0, root could not be found; 1 otherwise c c I/O integer n,type,rootf double precision roots(1:Nqmx) c temp integer method integer i,Nlast c label integer strlen character*(Nchar_mx) label label='subroutine findroots' if (n.lt.1) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'wrong input parameter:' write(*,*) 'n=',n stop endif if (type.eq.1) then c Legendre quadrature Nlast=n else if (type.eq.2) then c Legendre-Radau quadrature Nlast=n-1 else if (type.eq.3) then c Legendre-Lobatto quadrature Nlast=n-2 endif do i=1,Nlast c initialization of numerical method to use if (n.gt.16) then method=0 ! method=0 means use dichotomy else method=1 ! method=1 means use newton endif 111 continue if (method.eq.1) then call newton(n,i,type,roots(i),rootf) if (rootf.eq.0) then c then try dichotomy method=0 goto 111 endif else call dichotomy_roots(n,i,type,roots(i),rootf) endif if (rootf.eq.0) then write(*,*) 'Warning from ',label(1:strlen(label)),' :' write(*,*) 'root index:',i,' could not be found' write(*,*) 'for polynomia index:',n stop endif enddo return end subroutine newton(n,i,type,root,rootf) implicit none include 'max.inc' c c Purpose: to find the i-th root of the n-th Legendre polynomia c c Inputs: c + n: order of the Legendre polynomia c + i: index of the root to find c + type: type of quadrature (1: Legendre, 2: Legendre-Radau, 3: Legendre-Lobatto) c c Outputs: c + root: value of the i-th root c + rootf: if 0, root could not be found; 1 otherwise c c I/O integer n,i,type double precision root integer rootf c temp double precision polynom double precision dpolynom_dx double precision f_init,dfdx_init double precision x_init,x0,x1,x,dx,eps,err,err_init double precision xa,xb double precision legendre double precision legendre_derivative double precision legendre_second_derivative double precision f,dfdx integer Niter_mx,iter,x0f,solf parameter(Niter_mx=500) parameter(eps=1.0D-10) ! accuracy over f(x) parameter(dx=1.0D-5) c label integer strlen character*(Nchar_mx) label label='subroutine newton' if (i.gt.n) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'wrong input parameters:' write(*,*) 'n:',n write(*,*) 'i:',i stop endif c x0 is an estimation of the initial guess for the i-th root c of the n-th Legendre polynomia if (type.eq.1) then c analytical formula for x0 values, in the case of Legendre polynomia only c call guess_x0_v2(n,i,x_init,x0f) call guess_x0(n,i,type,x_init,xa,xb,x0f) else c Legendre-Radau or Legendre-Lobatto quadratures call guess_x0(n,i,type,x_init,xa,xb,x0f) endif x0=x_init if (x0f.eq.0) then rootf=0 goto 666 endif f_init=polynom(n,type,x_init) dfdx_init=dpolynom_dx(n,type,x_init) err_init=dabs(f_init) err=err_init iter=0 do while (err.gt.eps) x=x0 f=polynom(n,type,x) dfdx=dpolynom_dx(n,type,x) if (dfdx.eq.0.0D+0) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'dfdx(x=',x,')=',dfdx stop else x1=x0-f/dfdx err=dabs(polynom(n,type,x1)) x0=x1 endif iter=iter+1 if (iter.gt.Niter_mx) then rootf=0 goto 666 endif enddo root=x0 rootf=1 666 continue return end subroutine guess_x0(n,i,type,x0,x1,x2,x0f) implicit none include 'max.inc' c c Purpose: to guess a value for initial x0 c c Inputs: c + n: order of the Legendre polynomia c + i: index of the i-th root for which x0 must be guessed c + type: type of quadrature (1: Legendre, 2: Legendre-Radau, 3: Legendre-Lobatto) c c Outputs: c + x0 c + x1,x2: interval of x in which x0 has been identified, when possible c + x0f: if 0, x0 could not be found; 1 otherwise c c I/O integer n,i,type double precision x0 integer x0f c temp integer neven double precision legendre double precision legendre_derivative double precision polynom double precision a,b double precision dx,x1,x2,f1,f2 integer p,k,nint,inv,invif c label integer strlen character*(Nchar_mx) label label='subroutine guess_x0' 20 format(4e17.9) a=-1.0D+0 b=1.0D+0 p=1000 nint=p*n call is_even(n,neven) dx=(b-a)/nint inv=0 invif=0 do k=1,nint x1=a+(k-1)*dx x2=a+k*dx f1=polynom(n,type,x1) f2=polynom(n,type,x2) c Debug c if (x1.gt.0.8D+0) then c write(*,20) x1,x2,f1,f2 c endif c Debug if (f1*f2.le.0.0D+0) then ! f1 and f2 are of opposite signs if ((neven.eq.0).and.(x1.eq.0.0D+0)) then c nothing to do: no supplementary inversion has been found else inv=inv+1 ! index of the inversion that has been found if (inv.eq.i) then invif=1 x0=(x1+x2)/2.0D+0 x0f=1 goto 666 endif endif endif enddo if (invif.eq.0) then write(*,*) 'Warning from ',label(1:strlen(label)),' :' write(*,*) 'x0 could not be guessed for:' write(*,*) 'n:',n,' i:',i c Debug c do k=1,nint c x1=a+(k-1)*dx c x2=a+k*dx c f1=(legendre(n-1,x1)+legendre(n,x1))/(1.0D+0+x1) c f2=(legendre(n-1,x2)+legendre(n,x2))/(1.0D+0+x2) c write(*,20) x1,x2,f1,f2 c enddo c stop c Debug x0f=0 endif 666 continue return end subroutine guess_x0_v2(n,i,x0,x0f) implicit none c c Purpose: to guess a value for initial x0 c c Inputs: c + n: order of the Legendre polynomia c + i: index of the i-th root for which x0 must be guessed c c Outputs: c + x0 c + x0f: if 0, x0 could not be found; 1 otherwise c c I/O integer n,i double precision x0 integer x0f c temp double precision pi pi=4.0D+0*datan(1.0D+0) x0=-dcos(pi*(4.0D+0*(i-1)+3.0D+0)/(4.0D+0*n+2.0D+0)) x0f=1 return end subroutine is_even(n,even) implicit none include 'max.inc' c c Purpose: to check if n is even or odd c c Inputs: c + n: n c c Outputs: c + even: 1 if n is even; 0 is n is odd c c I/O integer n,even c temp c label integer strlen character*(Nchar_mx) label label='subroutine is_even' if (dble(n)/2.0D+0.eq.int(dble(n)/2.0D+0)) then even=1 else even=0 endif return end double precision function polynom(n,type,x) implicit none include 'max.inc' c c Purpose: to compute the value of the polynom used for the quadrature c c Inputs: c + n: quadrature order c + type: type of quadrature (1: Legendre, 2: Legendre-Radau, 3: Legendre-Lobatto) c + x: value of x c c Outputs: c + polynom: f(x) c c I/O integer n,type double precision x c temp double precision legendre double precision legendre_derivative if (type.eq.1) then c Legendre quadrature polynom=legendre(n,x) else if (type.eq.2) then c Legendre-Radau quadrature polynom=(legendre(n-1,x)+legendre(n,x))/(1.0D+0+x) else if (type.eq.3) then c Legendre-Lobatto quadrature polynom=legendre_derivative(n-1,x) endif return end double precision function dpolynom_dx(n,type,x) implicit none include 'max.inc' c c Purpose: to compute the value of the derivative of the polynom used for the quadrature c c Inputs: c + n: quadrature order c + type: type of quadrature (1: Legendre, 2: Legendre-Radau, 3: Legendre-Lobatto) c + x: value of x c c Outputs: c + dpolynom_dx: d(f(x))/d(x) c c I/O integer n,type double precision x c temp double precision legendre double precision legendre_derivative double precision legendre_second_derivative if (type.eq.1) then c Legendre quadrature dpolynom_dx=legendre_derivative(n,x) else if (type.eq.2) then c Legendre-Radau quadrature dpolynom_dx=( & (x+1.0D+0)* & (legendre_derivative(n-1,x)+legendre_derivative(n,x)) & -(legendre(n-1,x)+legendre(n,x)) & )/ & ((1.0D+0+x)**2.0D+0) else if (type.eq.3) then c Legendre-Lobatto quadrature dpolynom_dx=legendre_second_derivative(n-1,x) endif return end subroutine dichotomy_roots(n,i,type,root,rootf) implicit none include 'max.inc' c c Purpose: to find a root using a dichotomy procedure c c Inputs: c + n: quadrature order c + i: index of the root to find c + type: type of quadrature (1: Legendre, 2: Legendre-Radau, 3: Legendre-Lobatto) c c Outputs: c + root: x solution c + rootf: 1 when a solution has been found; 0 otherwise c c I/O integer n,i,type double precision root integer rootf c temp double precision x_init,xa,xb integer iter,Niter_mx,x0f parameter(Niter_mx=1000) double precision polynom double precision err,eps parameter(eps=1.0D-8) double precision x0,x1,x2,f1,f2,f0 c label integer strlen character*(Nchar_mx) label label='subroutine dichotomy' if (i.gt.n) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'wrong input parameters:' write(*,*) 'n:',n write(*,*) 'i:',i stop endif c x0 is an estimation of the initial guess for the i-th root c of the n-th Legendre polynomia if (type.eq.1) then c analytical formula for x0 values, in the case of Legendre polynomia only c call guess_x0_v2(n,i,x_init,x0f) call guess_x0(n,i,type,x_init,xa,xb,x0f) else c Legendre-Radau or Legendre-Lobatto quadratures call guess_x0(n,i,type,x_init,xa,xb,x0f) c Debug c write(*,*) 'n=',n,' i=',i c write(*,*) 'xa=',xa c write(*,*) 'x_init=',x_init c write(*,*) 'xb=',xb c write(*,*) 'x0f=',x0f c Debug endif x0=x_init if (x0f.eq.0) then rootf=0 goto 666 endif x0=x_init x1=xa x2=xb f0=polynom(n,type,x0) f1=polynom(n,type,x1) f2=polynom(n,type,x2) err=dabs(f0) iter=0 do while (err.gt.eps) iter=iter+1 if (f1*f0.lt.0.0D+0) then x2=x0 else x1=x0 endif x0=(x1+x2)/2.0D+0 f0=polynom(n,type,x0) f1=polynom(n,type,x1) f2=polynom(n,type,x2) err=dabs(f0) if (iter.gt.Niter_mx) then rootf=0 goto 666 endif enddo rootf=1 root=x0 666 continue return end double precision function legendre(n,x) implicit none include 'max.inc' c c Purpose: to compute the value of the Nth Legendre polynom c for a given value of x c c Inputs: c + n: order of the polynom c + x: value of x c c Outputs: c + legendre c c I/O integer n double precision x c temp double precision P,Prec integer max_order,order double precision Pnp1,Pn,Pnm1 c label integer strlen character*(Nchar_mx) label label='function legendre' max_order=12 if (n.le.max_order) then legendre=P(n,x) else c we need to use the recurence property over Lengendre polynomials: c (n+1)P(n+1,x)=(2n+1)*x*P(n,x)-n*P(n-1,x) c since evaluation of P(n,x) is NOT accurate for n>max_order Pn=P(max_order,x) Pnm1=P(max_order-1,x) do order=max_order,n-1 Pnp1=Prec(Pnm1,Pn,order,x) Pnm1=Pn Pn=Pnp1 enddo legendre=Pnp1 endif return end double precision function legendre_derivative(n,x) implicit none include 'max.inc' c c Purpose: to compute the value of the derivative of the Nth Legendre polynom c for a given value of x c c Inputs: c + n: order of the polynom c + x: value of x c c Outputs: c + legendre_derivative: d(Pn(x))/d(x) c c I/O integer n double precision x c temp integer order double precision legendre double precision dPndx,dPndx_rec double precision Pn,dPn_dx,dPnm1_dx,dPnp1_dx integer max_order c label integer strlen character*(Nchar_mx) label label='function legendre_derivative' max_order=12 if (n.le.max_order) then legendre_derivative=dPndx(n,x) else dPn_dx=dPndx(max_order,x) dPnm1_dx=dPndx(max_order-1,x) do order=max_order,n-1 Pn=legendre(order,x) dPnp1_dx=dPndx_rec(Pn,dPn_dx,dPnm1_dx,order,x) dPnm1_dx=dPn_dx dPn_dx=dPnp1_dx enddo legendre_derivative=dPnp1_dx endif return end double precision function legendre_second_derivative(n,x) implicit none include 'max.inc' c c Purpose: to compute the value of the second derivative of the Nth Legendre polynom c for a given value of x c c Inputs: c + n: order of the polynom c + x: value of x c c Outputs: c + legendre_second_derivative: d²(Pn(x))/d(x²) c c I/O integer n double precision x c temp integer order double precision legendre_derivative double precision d2Pndx2,d2Pndx2_rec double precision dPn_dx,d2Pn_dx2,d2Pnm1_dx2,d2Pnp1_dx2 integer max_order c label integer strlen character*(Nchar_mx) label label='function legendre_second_derivative' max_order=12 if (n.le.max_order) then legendre_second_derivative=d2Pndx2(n,x) else d2Pn_dx2=d2Pndx2(max_order,x) d2Pnm1_dx2=d2Pndx2(max_order-1,x) do order=max_order,n-1 dPn_dx=legendre_derivative(n,x) d2Pnp1_dx2=d2Pndx2_rec(dPn_dx,d2Pn_dx2,d2Pnm1_dx2,order,x) d2Pnm1_dx2=d2Pn_dx2 d2Pn_dx2=d2Pnp1_dx2 enddo legendre_second_derivative=d2Pnp1_dx2 endif return end double precision function Prec(Pnm1,Pn,n,x) implicit none c c Purpose: to evaluate P(n+1,x) from P(n-1,x) and P(n,x) using the recurence property: c (n+1)P(n+1,x)=(2n+1)*x*P(n,x)-n*P(n-1,x) c c Inputs: c + Pnm1: P(n-1,x) c + Pn: P(n,x) c + n: n c + x: x c c Outputs: c + Prec: P(n+1,x) c c I/O double precision Pnm1,Pn,x integer n Prec=((2.0D+0*n+1.0D+0)*x*Pn-n*Pnm1)/(n+1.0D+0) return end double precision function P(n,x) implicit none include 'max.inc' c c Purpose: to compute the value of the Nth Legendre polynom c for a given value of x, and n<13 c c Inputs: c + n: order of the polynom c + x: value of x c c Outputs: c + P: Pn(x) c c I/O integer n double precision x c temp integer k integer c double precision kelem,sum c label integer strlen character*(Nchar_mx) label label='function P' sum=0.0D+0 do k=0,n kelem=(c(n,k)**2.0D+0)* & ((x-1.0D+0)**(n-k))* & (x+1.0D+0)**k sum=sum+kelem enddo P=sum/(2.0D+0)**n return end double precision function dPndx_rec(Pn,dPndx,dPnm1dx,n,x) implicit none c c Purpose: to compute the derivative of the n-order Legendre polynomia c This function has to be used instead of function dPndx for n>12. c c Inputs: c + Pn: P(n,x) c + dPndx: d(P(n,x))/d(x) c + dPnm1dx: d(P(n-1,x))/d(x) c + n: order of the polynomia c + x: value of x c c Outputs: c + dPndx_rec= d(Pn(x))/d(x) c integer n double precision x double precision Pn,dPndx,dPnm1dx dPndx_rec=( & (2*n+1)*(Pn+x*dPndx) & -n*dPnm1dx & ) & /(n+1) return end double precision function dPndx(n,x) implicit none c c Purpose: to compute the derivative of the n-order Legendre polynomia c This function has to be used for n<13 only. Use recurrence c property above. c c Inputs: c + n: order of the polynomia c + x: value of x c c Outputs: c + dPndx= d(Pn(x))/d(x) c c I/O integer n double precision x c temp integer c double precision kelem double precision sum integer k sum=0.0D+0 do k=0,n kelem=(c(n,k)**2.0D+0)* & ( & (n-k)*(x-1.0D+0)**(n-k-1)*(x+1.0D+0)**k & +k*(x+1.0D+0)**(k-1)*(x-1.0D+0)**(n-k) & ) sum=sum+kelem enddo dPndx=sum/(2.0D+0)**n return end double precision function d2Pndx2_rec(dPndx,d2Pndx2,d2Pnm1dx2,n,x) implicit none c c Purpose: to compute the second derivative of the n-order Legendre polynomia c This function has to be used instead of function dPndx for n>12. c c Inputs: c + dPndx: d(P(n,x))/d(x) c + d2Pndx2: d²(P(n,x))/d(x²) c + d2Pnm1dx2: d²(P(n-1,x))/d(x²) c + n: order of the polynomia c + x: value of x c c Outputs: c + d2Pndx2_rec= d²(Pn(x))/d(x²) c integer n double precision x double precision dPndx,d2Pndx2,d2Pnm1dx2 d2Pndx2_rec=( & (2*n+1)*(2.0D+0*dPndx+x*d2Pndx2) & -n*d2Pnm1dx2 & ) & /(n+1) return end double precision function d2Pndx2(n,x) implicit none c c Purpose: to compute the second derivative of the n-order Legendre polynomia c This function has to be used for n<13 only. Use recurrence c property above. c c Inputs: c + n: order of the polynomia c + x: value of x c c Outputs: c + d2Pndx2= d²(Pn(x))/d(x²) c c I/O integer n double precision x c temp integer c double precision kelem double precision sum integer k sum=0.0D+0 do k=0,n kelem=(c(n,k)**2.0D+0)* & ( & (n-k)*(n-k-1)*(x-1.0D+0)**(n-k-2)*(x+1.0D+0)**k & +k*(k-1)*(x+1.0D+0)**(k-2)*(x-1.0D+0)**(n-k) & +2.0D+0*k*(n-k)*(x-1.0D+0)**(n-k-1)*(x+1.0D+0)**(k-1) & ) sum=sum+kelem enddo d2Pndx2=sum/(2.0D+0)**n return end integer function c(n,p) implicit none include 'max.inc' c c Purpose: to compute C(n,p) c c Inputs: c + n,p: integers for which C(n,p) has to be computed c c Outputs: c + c: value of C(n,p) c c I/O integer n,p c temp integer fact c label integer strlen character*(Nchar_mx) label label='function c' if (n.lt.p) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'n:',n,' < p:',p stop else c=fact(n)/(fact(p)*fact(n-p)) endif return end subroutine symetry(n,w,sym) implicit none include 'max.inc' c c Purpose: to check symetry of an array of values, i.e.: c w(1)=w(n) c w(2)=w(n-1) c etc... c c Inputs: c + n: size of the array c + w: array of values c c Outputs: c + sym: 1 if array "w" is symetric; 0 otherwise c c I/O integer n double precision w(1:Nqmx) integer sym c temp double precision relative_difference double precision err,eps integer neven,last,i parameter(eps=1.0D-2) ! maximal relative error allowed c label integer strlen character*(Nchar_mx) label label='subroutine symetry' if (n.lt.1) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'wrong input parameter n:',n stop endif if (n.eq.1) then sym=1 goto 666 else call is_even(n,neven) if (neven.eq.1) then last=n/2 else last=(n-1)/2 endif sym=1 do i=1,last err=relative_difference(w(i),w(n-i+1)) if (err.gt.eps) then sym=0 c Debug c write(*,*) w(i),w(n-i+1),err c Debug goto 666 endif enddo ! i endif 666 continue c if (sym.eq.0) then c write(*,*) 'Warning from ',label(1:strlen(label)),' :' c write(*,*) 'symetry was not found' c endif return end double precision function relative_difference(x1,x2) implicit none c c Purpose: to compute the relative difference between two given values c c Inputs: c + x1,x2: values whose relative difference has to be computed c c Outputs: c + relative_difference c c I/O double precision x1,x2 relative_difference=dabs((x1-x2)/x1) return end integer function fact(n) implicit none include 'max.inc' c c Purpose: to compute the factorial of a given integer c c Inputs: c + n: integer whose factorial has to be evaluated c c Outputs: c + fact=n! c c I/O integer n c temp integer i,f c label integer strlen character*(Nchar_mx) label label='function fact' if (n.lt.0) then write(*,*) 'Error from ',label(1:strlen(label)),' :' write(*,*) 'n=',n,' is negative' stop endif if (n.eq.0) then fact=1 else if (n.eq.1) then fact=1 else f=1 do i=1,n f=f*i enddo fact=f endif return end subroutine write_result(file,n,type,g,w) implicit none include 'max.inc' c c Purpose: to write results to file c c Inputs: c + file: name of the file to generate c + n: quadrature order c + type: quadrature type c 1: Legendre quadrature c 2: Legendre-Radau quadrature c 3: Legendre-Lobatto quadrature c + g: quadrature abscissas c + w: quadrature weights c c I/O character*(Nchar_mx) file integer n integer type double precision g(1:Nqmx) double precision w(1:Nqmx) double precision w_norm(1:Nqmx) c temp double precision sum_w integer i,sym c label integer strlen character*(Nchar_mx) label label='subroutine write_result' c Computation of the sum of quadrature weights sum_w=0.0D+0 do i=1,n sum_w=sum_w+w(i) enddo c Normalization do i=1,n w_norm(i)=w(i)/sum_w enddo c Symetry check call symetry(n,w,sym) open(12,file=file(1:strlen(file))) write(12,*) 'Quadrature order:' write(12,*) n write(12,*) 'Quadrature type:' if (type.eq.1) then write(12,*) 'Legendre' else if (type.eq.2) then write(12,*) 'Legendre-Radau' else if (type.eq.3) then write(12,*) 'Legendre-Lobatto' else if (type.eq.4) then write(12,*) 'Jacobi' endif write(12,*) write(12,*) 'i / g(i) / w(i) / w(i) normalized' do i=1,n write(12,*) g(i),w(i),w_norm(i) enddo write(12,*) write(12,*) 'sum of weights (without normalization):' write(12,*) sum_w write(12,*) write(12,*) 'symetry check:' if (sym.eq.1) then write(12,*) 'ok' else write(12,*) 'no symetry' endif close(12) return end