(Beginning of README.TXT File) This disk contains the following six files: README.TXT - the file you are currently reading TRI_DVR.FOR - driver test program for double integration over a triangle F_UV.FOR - integrand for the double integral LG.FOR - subroutine to compute Gauss-Legendre abscissas/weights TRI.FOR - subroutine for Algorithm 7, computation of double integrals over a triangle TRI_DVR.EXE - IBM PC compiled version of previous four files (End of README.TXT File) ------------------------------------------------------------------------------- (Beginning of TRI.FOR File) C C------------------------------DOCUMENTATION---------------------------- C SUBROUTINE DTRIA C C-----PURPOSE C C DTRIA COMPUTES AN APPROXIMATION TO THE DOUBLE INTEGRAL OF C F(U,V) OVER A TRIANGLE IN THE UV-PLANE BY USING AN N**2 POINT, C PRECISION 2*N-2, GENERALIZED GAUSS-LEGENDRE PRODUCT RULE C C-----USAGE C C CALL DTRIA(N,U1,V1,U2,V2,U3,V3,ABSC,WGHT,APROX,F) C C-----DESCRIPTION OF INPUT PARAMETERS C C N---------NUMBER OF POINTS IN THE GENERATING N-POINT C GAUSS-LEGENDRE QUADRATURE RULE ON (-1,1) C U1,V1----- C U2,V2----- C U3,V3-----DOUBLE PRECISION COORDINATES OF THE VERTICES OF THE C TRIANGLE C ABSC------ C WGHT------DOUBLE PRECISION ARRAYS CONTAINING THE NONNEGATIVE C ABSCISSAS AND CORRESPONDING WEIGHTS FOR THE N-POINT C GAUSS-LEGENDRE QUADRATURE RULE ON (-1,1) C C-----DESCRIPTION OF OUTPUT VARIABLES C C APPROX----DOUBLE PRECISION APPROXIMATION TO THE DOUBLE C INTEGRAL OVER THE TRIANGLE C C-----RESTRICTIONS C C THE ARRAYS ABSC AND WGHT MUST BE DIMENSIONED TO AT LEAST THE C INTEGER PART OF (N+1)/2 IN THE CALLING PROGRAM. IF N IS ODD, C THE ZERO ABSCISSA AND CORRESPONDING WEIGHT OF THE N-POINT C GAUSS-LEGENDRE RULE MUST BE IN THE (N+1)/2 POSITION OF THE C ARRAYS ABSC AND WGHT C C-----FUNCTIONS REQUIRED C C F---------USER SUPPLIED DOUBLE PRECISION FUNCTION F(U,V) OF C THE DOUBLE PRECISION VARIABLES U AND V. F MUST BE C LISTED IN AN EXTERNAL STATEMENT IN THE CALLING C PROGRAM C C-----REFERENCE C C LETHER, F.G., COMPUTATION OF DOUBLE INTEGRALS OVER A TRIANGLE, C ALGORITHM 007, J. COMP. APPL. MATH. 2(1976), 219-224. C C----------------------------------------------------------------------- C SUBROUTINE DTRIA(N,U1,V1,U2,V2,U3,V3,ABSC,WGHT,APROX,F) C DOUBLE PRECISION U1,U2,U3,V1,V2,V3,U21,U31,V21,V31,SAVE1,SAVE2, * SAVE3,SAVE4,SAVE5,SAVE6,X,Y,XX,YY,ABSJAC,APROX, * TEMP,HOLD1,HOLD2,ABSC(1),WGHT(1),F C C-----COMPUTE COEFFICIENTS FOR AFFINE TRANSFORMATION AND ABSOLUTE C VALUE OF JACOBIAN C U21 = U2 - U1 V21 = V2 - V1 U31 = U3 - U1 V31 = V3 - V1 ABSJAC = ABS(U21*V31 - V21*U31) C C-----INITIALIZE PARAMETERS C U21 = U21 * 0.5D0 V21 = V21 * 0.5D0 U31 = U31 * 0.25D0 V31 = V31 * 0.25D0 C NN = N NNP1 = NN + 1 ITEST = NNP1/2 APROX = 0.D0 C C-----CALCULATE CUBATURE SUM C DO 600 L = 1, NN C ---HOLD1 AND WGHT(I) CONTAIN THE VALUE OF ABSC(L) AND WGHT(L) C ---SAVE CONSTANTS TO AVOID REPETITIOUS CALCULATION C IF (L.LE.ITEST) GO TO 100 I = NNP1 - L HOLD1 = -ABSC(I) GO TO 200 100 I = L HOLD1 = ABSC(I) 200 SAVE1 = 1.D0 + HOLD1 SAVE2 = 1.D0 - HOLD1 SAVE3 = U31*SAVE1 SAVE4 = V31*SAVE1 SAVE5 = SAVE1*WGHT(I) XX = U21*SAVE2 + U1 YY = V21*SAVE2 + V1 TEMP = 0.D0 C DO 500 M=1,NN C C ---HOLD2 AND WGHT(J) CONTAIN THE VALUE OF ABSC(M) AND WGHT(M) C IF (M.LE.ITEST) GO TO 300 J = NNP1 - M HOLD2 = - ABSC(J) GO TO 400 300 J = M HOLD2 = ABSC(J) 400 SAVE6 = 1.D0 - HOLD2 X = XX + SAVE3*SAVE6 Y = YY + SAVE4*SAVE6 TEMP = TEMP + WGHT(J)*F(X,Y) 500 CONTINUE C APROX = APROX + SAVE5*TEMP C 600 CONTINUE C APROX = APROX*ABSJAC*0.125D0 C RETURN END (End of TRI.FOR File) ------------------------------------------------------------------------------- (Beginning of TRI_DVR.FOR File) C C Driver program to test SUBROUTINE DTRIA by approximating C the double integral of cos(u-v) over the triangle with C vertices at (0,0), (1,1) and (0,1). The exact value of C this test double integral is 1 - cos 1 C EXTERNAL F DOUBLE PRECISION X(25),W(25),EPSMCH,U1,V1,U2,V2,U3,V3, * APROX,EXACT,F REAL DIFFX,DIFFW INTEGER N,NN,ITEST C DATA ITEST / 1 / DATA U1,V1,U2,V2,U3,V3 / 0.D0,0.D0,1.D0,1.D0,0.D0,1.D0 / C C-----Determine the machine epsilon for SUBROUTINE GAUSSL C EPSMCH = 1.D0 10 EPSMCH = EPSMCH / 2.D0 IF ((1.0D0 + EPSMCH) .GT. 1.D0 ) GOTO 10 EPSMCH = 2.D0*EPSMCH C C-----Enter desired number of points for the n-point Gauss-Legendre C quadrature rule C 20 WRITE(*,*) ' ' WRITE(*,*) 'ENTER N .LE. 50 (0 to stop) :' READ(*,*) N IF (N .EQ. 0) GOTO 40 C C-----Compute nonnegative abscissas and corresponding weights for C the n-point Gauss-Legendre quadrature rule on (-1,1) C CALL GAUSSL (N,EPSMCH,ITEST,X,W,DIFFX,DIFFW) C C-----Approximate the double integral using a generalized C product cubature rule C CALL DTRIA(N,U1,V1,U2,V2,U3,V3,X,W,APROX,F) C C-----Write approximation and exact answer for test case C WRITE(*,*) ' ' WRITE(*,*) 'APROX = ',APROX EXACT = 1.D0 - DCOS(1.D0) WRITE(*,*) 'EXACT = ',EXACT C GOTO 20 C 40 STOP END (End of TRI_DVR.FOR File) ------------------------------------------------------------------------------- (Beginning of F_UV.FOR File) C C------------------------------------------------------------------- C DOUBLE PRECISION FUNCTION F(U,V) C C-----Integrand for the double integral C DOUBLE PRECISION U,V F = COS(U -V) RETURN END (End of F_UV.FOR File) ------------------------------------------------------------------------------- (Beginning of LG.FOR File) C C--------------------------- DOCUMENTATION --------------------------- C SUBROUTINE GAUSSL C C-----PURPOSE C C LET C NN=INTEGER PART OF (N+1)/2 C GAUSSL COMPUTES THE NN NONNEGATIVE ABSCISSAS AND CORRESPONDING C WEIGHTS FOR THE N-POINT GAUSS-LEGENDRE QUADRATURE RULE ON C (-1,1) C C-----USAGE C C CALL GAUSSL(N,EPSMCH,ITEST,X,W,DIFFX,DIFFW) C C-----DESCRIPTION OF INPUT PARAMETERS C C N---------NUMBER OF POINTS IN THE GAUSS-LEGENDRE RULE ON (-1,1) C EPSMCH----UNIT ROUND-OFF OF THE MACHINE. EPSMCH IS A MACHINE C DEPENDENT PARAMETER WHICH MAY BE TAKEN AS THE C SMALLEST DOUBLE PRECISION NUMBER FOR WHICH C 1.D0 + EPSMCH .GT. 1.D0 C IS TRUE IN THE DOUBLE PRECISION ARITHMETIC OF THE C COMPUTER BEING USED C ITEST-----USER SPECIFIED TEST PARAMETER C IF C ITEST=0 , C GAUSSL COMPUTES THE RELATIVE ERROR DIFFX IN THE SUM C OF SQUARES OF THE NONNEGATIVE ABSCISSAS, AND THE C RELATIVE ERROR DIFFW IN THE SUM OF THE CORRESPONDING C WEIGHTS C IF C ITEST=1 , C DIFFX AND DIFFW ARE NOT COMPUTED C C-----DESCRIPTION OF OUTPUT PARAMETERS C C X--------- C W---------DOUBLE PRECISION ARRAYS CONTAINING THE NONNEGATIVE C ABSCISSAS C X(NN) < ... < X(2) < X(1) < 1 C AND CORRESPONDING WEIGHTS C W(NN) , ... , W(2) , W(1) C DIFFX-----(COMPUTED) RELATIVE ERROR IN C X(1)**2 + X(2)**2 +...+ X(NN)**2 C DIFFW-----(COMPUTED) RELATIVE ERROR IN C W(1) + W(2) +...+ W(NN) . C DIFFX AND DIFFW INDICATE (APPROXIMATELY) THE ACCURACY C TO WHICH THE ABSCISSAS X(K) AND WEIGHTS W(K) HAVE C BEEN COMPUTED BY GAUSSL C DIFFX AND DIFFW ARE SINGLE PRECISION VARIABLES C C-----RESTRICTIONS C C THE DOUBLE PRECISION ARRAYS X AND W MUST BE DIMENSIONED AT C LEAST TO THE INTEGER PART OF (N+1)/2 IN THE CALLING PROGRAM C C-----REFERENCE C C LETHER, F.G., ON THE CONSTRUCTION OF GAUSS-LEGENDRE QUADRATURE C RULES, J. COMP. APPL. MATH., 4(1978), 47-52 C C----------------------------------------------------------------------- C SUBROUTINE GAUSSL(N,EPSMCH,ITEST,X,W,DIFFX,DIFFW) C INTEGER ITEST,ITR,K,KK,KLAST,M,MLAST,MZERO,N,NM1,NPOINT, * NSAVE,MOD REAL ALPHAN,ANGLE,BZERO,CNST1,CNST2,DIFFW,DIFFX,EPSLON, * FACTOR,VNM DOUBLE PRECISION C1,C2,C3,C4,C5,DNOM,EPSMCH,FLOAT1,FLOAT2, * ONEMX,PN,PNM1, PNM2,PROD,SUMW,SUMXSQ,TEMPX, * TERM1,TERM2,U,V,W(1),X(1),XNM,XNMSQ,XPNM1, * DATAN,DSQRT,DBLE,DFLT C C-----INITIALIZATION FOR THE FIRST POSITIVE ZERO 2.404825557695772769 OF C THE BESSEL FUNCTION OF THE FIRST KIND OF ORDER ZERO. BZERO SHOULD C BE GIVEN TO THE SINGLE PRECISION ACCURACY OF THE COMPUTER BEING C USED. C DATA BZERO/2.4048256/ C C-----DEFINE NON-ANSI FORTRAN DOUBLE PRECISION FLOAT FUNCTION. C DFLT(INT)=FLOAT(INT) C C-----IF N=1,2,3 COMPUTE NONNEGATIVE ABSCISSAS AND WEIGHTS, THEN RETURN. C IF(N.GE.4) GO TO 4 DIFFX=0.E0 DIFFW=0.E0 GO TO (1,2,3),N C 1 X(1)=0.D0 W(1)=2.D0 GO TO 800 C 2 X(1)=SQRT(3.D0)/3.D0 W(1)=1.D0 GO TO 800 C 3 X(1)=SQRT(15.D0)/5.D0 X(2)=0.D0 W(1)=5.D0/9.D0 W(2)=8.D0/9.D0 GO TO 800 C C-----INITIALIZATION. C C INTEGER PARAMETERS. C 4 MLAST=N/2 NM1=N-1 NPOINT=N/3 C C DOUBLE PRECISION PARAMETERS. C FLOAT1=DFLT(N*N+N) FLOAT2=5.D0*FLOAT1 C1=(3.D0+FLOAT1)/3.D0 C2=(1.D0+FLOAT1)/3.D0 C3=(6.D0+FLOAT2)/6.D0 C4=(4.D0+FLOAT2)/6.D0 C5=4.D0*ATAN(1.D0)/DFLT(4*N+2) C C SINGLE PRECISION PARAMETERS. C CNST1=1.D0-0.125D0*(1.D0-1.D0/DFLT(N))/DFLT(N*N) CNST2=1.D0/DFLT(N)**4/384.D0 C C-----COMPUTE (GATTESCHI) INITIAL GUESS FOR LARGEST ABSCISSA X(1). C ALPHAN=FLOAT1+1.D0/3.D0 ANGLE=BZERO/SQRT(ALPHAN) ANGLE=ANGLE*(1.D0-(BZERO*BZERO-2.D0)/(360.D0*ALPHAN*ALPHAN)) XNM=COS(ANGLE) C C MAIN LOOP TO COMPUTE POSITIVE ABSCISSAS AND WEIGHTS. C ****************************************************************** DO 400 M=1,MLAST C ITR=0 IF(M.EQ.1) GO TO 100 C ---COMPUTE (TRICOMI) INITIAL GUESS FOR M-TH ABSCISSA IF M > 1 . C VNM=FLOAT(4*M-1)*C5 FACTOR=0.D0 IF(M.LE.NPOINT) FACTOR=39.D0-28.D0/SIN(VNM)**2 XNM=(CNST1-CNST2*FACTOR)*COS(VNM) C C ---CALCULATE LEGENDRE POLNOMIALS BY THREE TERM RECURRENCE C RELATION. ALSO COMPUTE DENOMINATOR NEEDED FOR THE M-TH WEIGHT. C 100 PNM2=1.D0 PNM1=XNM XNMSQ=XNM*XNM DNOM=1.D0+3.D0*XNMSQ DO 200 K=2,NM1 XPNM1=XNM*PNM1 PN=XPNM1+XPNM1-PNM2-(XPNM1-PNM2)/DFLT(K) DNOM=DNOM+DFLT(K+K+1)*PN*PN PNM2=PNM1 PNM1=PN 200 CONTINUE XPNM1=XNM*PNM1 PN=XPNM1+XPNM1-PNM2-(XPNM1-PNM2)/DFLT(N) C C ---CHECK ACCURACY OF XNM. C TRY TO AVOID LOSS OF ACCURACY IN COMPUTING 1-X . C V=PN/(PNM1-XNM*PN)/DFLT(N) ONEMX=1.D0-XNM TEMPX=(1.D0-ONEMX)-XNM ONEMX=ONEMX+TEMPX U=ONEMX*(1.D0+XNM)*V EPSLON=2.D0*ABS(XNM)*EPSMCH IF(ABS(U).LE.EPSLON) GO TO 300 C C ---LIMIT NUMBER OF ITERATIONS TO AT MOST 2. C ITR=ITR+1 IF(ITR.GT.2) GO TO 300 C C ---COMPUTE 5-TH ORDER ITERATION. C TERM1=C1*XNMSQ-C2 TERM2=(C3*XNMSQ-C4)*XNM XNM=XNM-U*(1.D0+V*(XNM+V*(TERM1+V*TERM2))) GO TO 100 C C ---LOAD M-TH POSITIVE ABSCISSA AND DETERMINE CORRESPONDING WEIGHT. C 300 X(M)=XNM W(M)=2.D0/DNOM C WRITE(5,*)M,X(M),W(M) C 400 CONTINUE C ****************************************************************** C END MAIN LOOP FOR COMPUTING POSITIVE ABSCISSAS AND WEIGHTS. C C-----IF N IS ODD COMPUTE THE WEIGHT CORRESPONDING TO THE ZERO ABSCISSA. C NSAVE=MOD(N,2) IF(NSAVE.EQ.0) GO TO 600 MZERO=MLAST+1 X(MZERO)=0.D0 KLAST=MLAST-1 PROD=2.D0/DFLT(N) DO 500 K=1,KLAST PROD=PROD*DFLT(K+K+2)/DFLT(K+K+1) 500 CONTINUE W(MZERO)=2.D0*PROD*PROD C C-----IF REQUESTED, COMPUTE RELATIVE ERRORS IN TEST SUMS FOR NONNEGATIVE C ABSCISSAS AND CORRESPONDING WEIGHTS. SUM TERMS IN INCREASING ORDER C OF MAGNITUDE. C 600 IF(ITEST.EQ.1) GO TO 800 SUMXSQ=0.D0 SUMW=0.D0 NPOINT=MLAST+1 DO 700 K=1,MLAST KK=NPOINT-K SUMXSQ=SUMXSQ+X(KK)*X(KK) 700 SUMW=SUMW+W(K) IF(NSAVE.GT.0) SUMW=SUMW+0.5D0*W(MZERO) DIFFW=ABS(SUMW-1.D0) DIFFX=ABS(DFLT(4*N-2)/DFLT(N*N-N)*SUMXSQ-1.D0) 800 RETURN END (End of LG.FOR File)