C****************************************************************************** C * C CUBIC SPLINE RAYLEIGH-RITZ METHOD * C * C****************************************************************************** C C C TO APPROXIMATE THE SOLUTION TO THE BOUNDARY-VALUE PROBLEM C C -D(P(X)Y')/DX + Q(X)Y = F(X), 0<=X<=1, Y(0)=Y(1)=0 C C WITH A SUM OF CUBIC SPLINES: C C INPUT INTEGER N C C OUTPUT COEFFICIENTS C(0),...,C(N+1) OF THE BASIS FUNCTIONS C C (NOTE THAT C(I) USED HERE IS 6 TIMES AS LARGE AS THE C(I) C USED IN THE ALGORITHM, SINCE PHI(X) IS DEFINED TO BE 1/6 C OF THE BASIS FUNCTION USED IN THE TEXT.) C C TO CHANGE PROBLEMS DO THE FOLLOWING: C 1. CHANGE FUNCTIONS P, Q AND F IN THE SUBPROGRAMS NAMED P,Q,F C 2. CHANGE N C 3. CHANGE DIMENSIONS IN MAIN PROGRAM TO VALUES GIVEN BY C DIMENSION A(N+2,N+3),X(N+2),C(N+2),CO(N+2,4,4), C DCO(N+2,4,3),AF(N+1),BF(N+1),CF(N+1),DF(N+1), C AP(N+1),BP(N+1),CP(N+1),DP(N+1), C AQ(N+1),BQ(N+1),CQ(N+1),DQ(N+1) C 4. CHANGE DIMENSIONS IN SUBROUTINE COEF TO C DIMENSION AA(N+2),BB(N+2),CC(N+2),DD(N+2), C H(N+2),XA(N+2),XL(N+2),XU(N+2),XZ(N+2) C C C GENERAL OUTLINE C C 1. NODES LABELLED X(I)=(I-1)*H, 1 <= I <= N+2, WHERE C H=1/(N+1) SO THAT ZERO SUBSCRIPTS ARE AVOIDED C 2. THE FUNCTIONS PHI(I) AND PHI'(I) ARE SHIFTED SO THAT C PHI(1) AND PHI'(1) ARE CENTERED AT X(1), PHI(2) AND PHI'(2) C ARE CENTERED AT X(2), . . . , PHI(N+2) AND C PHI'(N+2) ARE CENTERED AT X(N+2)---FOR EXAMPLE, C PHI(3) = S((X-X(3))/H) C = S(X/H + 2) C 3. THE FUNCTIONS PHI(I) ARE REPRESENTED IN TERMS OF THEIR C COEFFICIENTS IN THE FOLLOWING WAY: C (PHI(2))(X) = CO(I,K,1) + CO(I,K,2)*X + CO(I,K,3)*X**2 C CO(I,K,4)*X**3 C FOR X(J) <= X <= X(J+1) WHERE C K=1 IF J=I-2, K=2 IF J=I-1, K=3 IF J=I, K=4 IF J=I+1 C SINCE PHI(I) IS NONZERO ONLY BETWEEN X(I-2) AND X(I+2) C UNLESS I = 1, 2, N+1 OR N+2 C (SEE SUBROUTINE PHICO) C 4. THE DERIVATIVE OF PHI(I) DENOTED PHI'(I) IS REPRESENTED C AS IN 3. BY ITS COEFFICIENTS DCO(I,K,L), L = 1, 2, 3 C (SEE SUBROUTINE DPHICO). C 5. THE FUNCTIONS P,Q AND F ARE REPRESENTED BY THEIR CUBIC C SPLINE INTERPOLANTS USING CLAMPLED BOUNDARY CONDITIONS C (SEE ALGORITHM 3.3). THUS, FOR X(I) <= X <= X(I+1) WE C USE AF(I) + BF(I)*X + CF(I)*X**2 + DF(I)*X**3 TO C REPRESENT F(X). SIMILARLY, AP,BP,CP,DP ARE USED FOR P AND C AQ,BQ,CQ,DQ ARE USED FOR Q. (SEE SUBROUTINE COEF). C 6. THE INTEGRANDS IS STEPS 6 AND 9 ARE REPLACED BY PRODUCTS C OF CUBIC POLYNOMIAL APPROXIMATIONS ON EACH SUBINTERVAL OF C LENGTH H AND THE INTEGRALS OF THE RESULTING POLYNOMIALS C ARE COMPUTED EXACTLY. (SEE SUBROUTINE XINT). C C C IMPLICIT REAL*8(A-H,O-Z) EXTERNAL P,Q,F DIMENSION A(11,12),X(11),C(11),CO(11,4,4),DCO(11,4,3),AF(10) DIMENSION DF(10),AP(10),BP(10),CP(10),DP(10),AQ(10),BQ(10),DQ(10) DIMENSION BF(10),CF(10),CQ(10) CHARACTER NAME*14,NAME1*14,AA*1 INTEGER INP,OUP,FLAG LOGICAL OK COMMON N OPEN(UNIT=5,FILE='CON',ACCESS='SEQUENTIAL') OPEN(UNIT=6,FILE='CON',ACCESS='SEQUENTIAL') WRITE(6,*) 'This is the Cubic Spline Rayleigh-Ritz Method.' WRITE(6,*) 'Have the functions P, Q and F been created? ' WRITE(6,*) 'Enter Y or N ' WRITE(6,*) ' ' READ(5,*) AA IF(( AA .EQ. 'Y' ) .OR. ( AA .EQ. 'y' )) THEN OK = .FALSE. 12 IF (OK) GOTO 13 WRITE(6,*) 'Input a positive integer n where' WRITE(6,*) 'x(0) = 0, ..., x(n+1) = 1.' WRITE(6,*) ' ' READ(5,*) N IF ( N .LE. 0 ) THEN WRITE(6,*) 'Must be positive integer ' ELSE OK = .TRUE. ENDIF GOTO 12 13 WRITE(6,*) 'Input the derivative of F at 0 and at 1.' WRITE(6,*) ' ' READ(5,*) FPL, FPR WRITE(6,*) 'Input the derivative of P at 0 and at 1.' WRITE(6,*) ' ' READ(5,*) PPL, PPR WRITE(6,*) 'Input the derivative of Q at 0 and at 1.' WRITE(6,*) ' ' READ(5,*) QPL, QPR ELSE WRITE(6,*) 'The program will end so that P, Q and F ' WRITE(6,*) 'can be created ' OK = .FALSE. ENDIF IF(.NOT.OK) GOTO 400 WRITE(6,*) 'Select output destinations: ' WRITE(6,*) '1. Screen ' WRITE(6,*) '2. Text file ' WRITE(6,*) 'Enter 1 or 2 ' WRITE(6,*) ' ' READ(5,*) FLAG IF ( FLAG .EQ. 2 ) THEN WRITE(6,*) 'Input the file name in the form - ' WRITE(6,*) 'drive:name.ext' WRITE(6,*) 'with the name contained within quotes' WRITE(6,*) 'as example: ''A:OUTPUT.DTA'' ' WRITE(6,*) ' ' READ(5,*) NAME1 OUP = 3 OPEN(UNIT=OUP,FILE=NAME1,STATUS='NEW') ELSE OUP = 6 ENDIF WRITE(OUP,*) 'CUBIC SPLINE RAYLEIGH-RITZ METHOD' C STEP 1 H=1.0/(N+1) N1=N+1 N2=N+2 N3=N+3 C INITIALIZE MATRIX A AT ZERO, NOTE THAT A(I,N+3) = B(I) DO 20 I=1,N2 DO 20 J=1,N3 20 A(I,J)=0 C C STEP 2 C C X(1)=0,...,X(I)=(I-1)*H,...,X(N+2)=1 DO 30 I=1,N2 30 X(I)=(I-1)*H C C STEPS 3 AND 4 ARE IMPLEMENTED IN WHAT FOLLOWS C C INITIALIZE COEFFICIENTS CO(I,J,K) AND DCO(I,J,K) DO 40 I=1,N2 DO 50 J=1,4 DO 60 K=1,4 CO(I,J,K)=0 IF(K.NE.4) THEN DCO(I,J,K)=0 END IF 60 CONTINUE C JJ CORRESPONDS THE COEFFICIENTS OF PHI AND PHI' TO THE PROPER C INTERVAL INVOLVING J JJ=I+J-3 CALL PHICO(I,JJ,CO(I,J,1),CO(I,J,2),CO(I,J,3),CO(I,J,4)) CALL DPHICO(I,JJ,DCO(I,J,1),DCO(I,J,2),DCO(I,J,3)) 50 CONTINUE 40 CONTINUE C OUTPUT THE BASIS FUNCTIONS WRITE(OUP,1) WRITE(OUP,2) DO 70 I=1,N2 WRITE(OUP,5) I,I DO 70 J=1,4 IF(I.NE.1.OR.(J.NE.1.AND.J.NE.2)) THEN IF(I.NE.2.OR.J.NE.2) THEN IF(I.NE.N1.OR.J.NE.4) THEN IF(I.NE.N2.OR.(J.NE.3.AND.J.NE.4)) THEN JJ1=I-3+J JJ2=I-2+J WRITE(OUP,3) JJ1,JJ2,(CO(I,J,K),K=1,4),(DCO(I,J,K),K=1,3) END IF END IF END IF END IF 70 CONTINUE C C OBTAIN COEFFICIENTS FOR F, P AND Q--NOTE THAT THE 2ND AND C 3RD ARGUMENTS ARE THE DERIVATIVES AT 0 AND 1 RESP. C CALL COEF(F,FPL,FPR,X,AF,BF,CF,DF,N2,N1) CALL COEF(P,PPL,PPR,X,AP,BP,CP,DP,N2,N1) CALL COEF(Q,QPL,QPR,X,AQ,BQ,CQ,DQ,N2,N1) C C STEPS 5, 6, 7, 8, 9 ARE IMPLEMENTED IN WHAT FOLLOWS C WRITE(OUP,6) DO 80 I=1,N2 C INDICES OF LIMITS OF INTEGRATION FOR A(I,I) AND B(I) J1=MIN0(I+2,N+2) JO=MAX0(I-2,1) J2=J1-1 C INTEGRATE OVER EACH SUBINTERVAL WHERE PHI(I) NONZERO DO 90 JJ=JO,J2 C LIMITS OF INTEGRATION FOR EACH CALL XU=X(JJ+1) XL=X(JJ) C COEFFICIENTS OF BASES K=INTY(I,JJ) A1=DCO(I,K,1) B1=DCO(I,K,2) C1=DCO(I,K,3) D1=0 A2=CO(I,K,1) B2=CO(I,K,2) C2=CO(I,K,3) D2=CO(I,K,4) C CALL SUBPROGRAM FOR INTEGRATION A(I,I)=A(I,I)+XINT(XU,XL,AP(JJ),BP(JJ),CP(JJ),DP(JJ),A1,B1,C1,D1, 1A1,B1,C1,D1)+XINT(XU,XL,AQ(JJ),BQ(JJ),CQ(JJ),DQ(JJ),A2,B2,C2,D2, 1A2,B2,C2,D2) 90 A(I,N+3)=A(I,N+3)+XINT(XU,XL,AF(JJ),BF(JJ),CF(JJ),DF(JJ),A2,B2,C2, 1D2,1.0D+00,0.0D+00,0.0D+00,0.0D+00) C COMPUTE A(I,J) FOR J = I+1, . . . , MIN(I+3,N+2) K3=I+1 IF(K3.LE.N2) THEN K2=MIN0(I+3,N+2) DO 100 J=K3,K2 JO=MAX0(J-2,1) DO 110 JJ=JO,J2 XU=X(JJ+1) XL=X(JJ) K=INTY(I,JJ) A1=DCO(I,K,1) B1=DCO(I,K,2) C1=DCO(I,K,3) D1=0 A2=CO(I,K,1) B2=CO(I,K,2) C2=CO(I,K,3) D2=CO(I,K,4) K=INTY(J,JJ) A3=DCO(J,K,1) B3=DCO(J,K,2) C3=DCO(J,K,3) D3=0 A4=CO(J,K,1) B4=CO(J,K,2) C4=CO(J,K,3) D4=CO(J,K,4) 110 A(I,J)=A(I,J)+XINT(XU,XL,AP(JJ),BP(JJ),CP(JJ),DP(JJ),A1,B1,C1,D1, 1A3,B3,C3,D3)+XINT(XU,XL,AQ(JJ),BQ(JJ),CQ(JJ),DQ(JJ),A2,B2,C2,D2,A4 1,B4,C4,D4) 100 A(J,I)=A(I,J) C OUTPUT A(I,J) FOR J=I,I+1,...,MIN(I+3,N+2) AND J=N+3 WRITE(OUP,7) (I,J,A(I,J),J=I,K2),I,N3,A(I,N3) END IF 80 CONTINUE C C STEP 10 C DO 120 I=1,N1 II=I+1 DO 130 J=II,N2 CC=A(J,I)/A(I,I) DO 140 K=II,N3 140 A(J,K)=A(J,K)-CC*A(I,K) 130 A(J,I)=0 120 CONTINUE C(N2)=A(N2,N3)/A(N2,N2) DO 150 I=1,N1 J=N1-I+1 C(J)=A(J,N3) JJ=J+1 DO 160 KK=JJ,N2 160 C(J)=C(J)-A(J,KK)*C(KK) 150 C(J)=C(J)/A(J,J) C C STEP 11 C C OUTPUT THE COEFFICIENTS C(I) WRITE(OUP,8) WRITE(OUP,9) (I,C(I),I=1,N2) C OBTAIN VALUES OF THE APPROX. AT THE NODES DO 170 I=1,N2 S=0 DO 180 J=1,N2 JO=MAX0(J-2,1) J1=MIN0(J+2,N+2) SS=0 IF(I.LT.JO.OR.I.GE.J1) THEN S=S+C(J)*SS ELSE K=INTY(J,I) SS=((CO(J,K,4)*X(I)+CO(J,K,3))*X(I)+CO(J,K,2))* * X(I)+CO(J,K,1) S=S+C(J)*SS END IF 180 CONTINUE WRITE(OUP,10) I,X(I),S 170 CONTINUE C C STEP 12 C PROCEDURE IS COMPLETE 400 CLOSE(UNIT=5) CLOSE(UNIT=OUP) IF(OUP.NE.6) CLOSE(UNIT=6) STOP 1 FORMAT(1X,'BASIS FUNCTION: A + B*X + C*X**2 + D*X**3',21X, 1 'DER. OF BASIS FUNCTION: A + B*X + C*X**2 ') 2 FORMAT(1X,'A',14X,'B',15X,'C',15X,'D',18X,'A',15X,'B',14X, 1 'C') 3 FORMAT(1X,'ON (X(',I2,'),X(',I2,'))',5X,4(D13.6,2X),3X,3(D13.6, 1 2X)) 5 FORMAT(1X,'PHI(',I2,')',77X,'PHI''(',I2,')') 6 FORMAT(1X,'NONZERO ENTRIES IN MATRIX A') 7 FORMAT(1X,5('A(',I2,',',I2,') = ',D15.8)) 8 FORMAT(1X,'SOLUTION IS C(1)*PHI(1) + ... + C(N+2)*PHI(N+2) 1WHERE ') 9 FORMAT(5(1X,'C(',I2,') = ', D15.8)) 10 FORMAT(1X,'VALUE OF APPROX. AT X(',I2,') = ',D15.8,3X,'IS ',D15.8) END C SUBROUTINE COEF(F,FPO,FPN,X,A,B,C,D,N,M) C C THIS IMPLEMENTS ALGORITHM 3.3--CLAMPED BOUNDARY CUBIC SPLINE C F = FUNCTION TO BE APPROXIMATED, FPO = F'(0), FPN = F'(1), C X IS THE SET OF NODES, A,B,C,D ARE THE COEFFICIENTS TO BE C RETURNED, N = NUMBER OF NODES, M = NUMBER OF SUBINTERVALS C THE APPROX. WILL BE C A(I) + B(I)*X + C(I)*X**2 + D(I)*X**3 ON (X(I),X(I+1)) IMPLICIT REAL*8(A-H,O-Z) DIMENSION AA(11),BB(11),CC(11),DD(11) DIMENSION X(N),A(M),B(M),C(M),D(M),H(11),XA(11),XL(11),XU(11) DIMENSION XZ(11) DO 10 I=1,M AA(I)=F(X(I)) 10 H(I+1)=X(I+1)-X(I) AA(N)=F(X(N)) XA(1)=3*(AA(2)-AA(1))/H(2)-3*FPO XA(N)=3*FPN-3*(AA(N)-AA(N-1))/H(N) XL(1)=2*H(2) XU(1)=.5D+00 XZ(1)=XA(1)/XL(1) DO 20 I=2,M XA(I)=3*(AA(I+1)*H(I)-AA(I)*(X(I+1)-X(I-1))+AA(I-1)*H(I+1)) * /(H(I+1)*H(I)) XL(I)=2*(X(I+1)-X(I-1))-H(I)*XU(I-1) XU(I)=H(I+1)/XL(I) 20 XZ(I)=(XA(I)-H(I)*XZ(I-1))/XL(I) XL(N)=H(N)*(2-XU(N-1)) XZ(N)=(XA(N)-H(N)*XZ(N-1))/XL(N) CC(N)=XZ(N) DO 30 I=1,M J=N-I CC(J)=XZ(J)-XU(J)*CC(J+1) BB(J)=(AA(J+1)-AA(J))/H(J+1)-H(J+1)*(CC(J+1)+2*CC(J))/3 30 DD(J)=(CC(J+1)-CC(J))/(3*H(J+1)) DO 40 I=1,M A(I)=(((-DD(I)*X(I))+CC(I))*X(I)-BB(I))*X(I)+AA(I) B(I)=(3*DD(I)*X(I)-2*CC(I))*X(I)+BB(I) C(I)=CC(I)-3*DD(I)*X(I) 40 D(I)=DD(I) RETURN END C FUNCTION XINT (XU,XL,A1,B1,C1,D1,A2,B2,C2,D2,A3,B3,C3,D3) C C FORMS THE PRODUCT (A1+B1*X+C1*X**2+D1*X**3)*(A2+B2*X+C2*X**2+ C D2*X**3)*(A3+B3*X+C3*X**2+D3*X**3) TO OBTAIN C 10*C(1)*X**9 + 9*C(2)*X**8 + . . . + 2*C(9)*X + C(10) C WHICH IS INTEGRATED FROM XL TO XU IMPLICIT REAL*8(A-H,O-Z) DIMENSION C(10) AA=A1*A2 BB=A1*B2+A2*B1 CC=A1*C2+B1*B2+C1*A2 DD=A1*D2+B1*C2+C1*B2+D1*A2 EE=B1*D2+C1*C2+D1*B2 FF=C1*D2+D1*C2 GG=D1*D2 C(10)=AA*A3 C(9)=(AA*B3+BB*A3)/2 C(8)=(AA*C3+BB*B3+CC*A3)/3 C(7)=(AA*D3+BB*C3+CC*B3+DD*A3)/4 C(6)=(BB*D3+CC*C3+DD*B3+EE*A3)/5 C(5)=(CC*D3+DD*C3+EE*B3+FF*A3)/6 C(4)=(DD*D3+EE*C3+FF*B3+GG*A3)/7 C(3)=(EE*D3+FF*C3+GG*B3)/8 C(2)=(FF*D3+GG*C3)/9 C(1)=(GG*D3)/10 XHIGH=0 XLOW=0 DO 10 I=1,10 XHIGH=(XHIGH+C(I))*XU 10 XLOW=(XLOW+C(I))*XL XINT=XHIGH-XLOW RETURN END C SUBROUTINE PHICO(I,J,A,B,C,D) C C COMPUTES PHI(I) AS A+B*X+C*X**2+D*X**3 FOR X IN (X(J),X(J+1)) BY C MULTIPLYING AND SIMPLIFYING THE EQUATIONS IN STEP 2 IMPLICIT REAL*8(A-H,O-Z) COMMON N H=1.0D+00/(N+1) A=0 B=0 C=0 D=0 E=I-1 IF(J.LT.I-2.OR.J.GE.I+2) RETURN IF(I.EQ.1.AND.J.LT.I) RETURN IF(I.EQ.2.AND.J.LT.I-1) RETURN IF(I.EQ.N+1.AND.J.GT.N+1) RETURN IF(I.EQ.N+2.AND.J.GE.N+2) RETURN IF(J.LE.I-2) THEN A= (((-E+6)*E-12)*E+8)/24 B=((E-4)*E+4)/(8*H) C=(-E+2)/(8*H*H) D=1/(24*H**3) RETURN ELSE IF(J.GT.I-1.AND.J.GT.I) THEN A=(((E+6)*E+12)*E+8)/24 B=((-E-4)*E-4)/(8*H) C=(E+2)/(8*H*H) D=-1/(24*H**3) RETURN ELSE IF(J.GT.I-1) THEN A=((-3*E-6)*E*E+4)/24 B=(3*E+4)*E/(8*H) C=(-3*E-2)/(8*H*H) D=1/(8*H**3) IF(I.NE.1.AND.I.NE.N+1) RETURN ELSE A=((3*E-6)*E*E+4)/24 B=(-3*E+4)*E/(8*H) C=(3*E-2)/(8*H*H) D=-1/(8*H**3) IF(I.NE.2.AND.I.NE.N+2) RETURN END IF END IF END IF IF(I.LE.2) THEN AA=1.0D+00/24 BB=-1/(8*H) CC=1/(8*H*H) DD=-1/(24*H**3) IF(I.EQ.2) THEN A=A-AA B=B-BB C=C-CC D=D-DD RETURN ELSE A=A-4*AA B=B-4*BB C=C-4*CC D=D-4*DD RETURN END IF ELSE EE=N+2 AA=(((-EE+6)*EE-12)*EE+8)/24 BB=((EE-4)*EE+4)/(8*H) CC=(-EE+2)/(8*H*H) DD=1/(24*H**3) IF(I.EQ.N+1) THEN A=A-AA B=B-BB C=C-CC D=D-DD RETURN ELSE A=A-4*AA B=B-4*BB C=C-4*CC D=D-4*DD END IF END IF RETURN END C SUBROUTINE DPHICO(I,J,A,B,C) C C SAME AS PHICO EXCEPT THE COEFFICIENTS ARE FOR PHI'(I) IMPLICIT REAL*8(A-H,O-Z) COMMON N H=1.0D+00/(N+1) A=0 B=0 C=0 E=I-1 IF(J.LT.I-2 .OR. J.GE.I+2) RETURN IF(I.EQ.1 .AND. J.LT.I) RETURN IF(I.EQ.2 .AND. J.LT.I-1) RETURN IF(I.EQ.N+1 .AND. J.GT.N+1) RETURN IF(I.EQ.N+2 .AND. J.GE.N+2) RETURN IF(J.LE.I-2) THEN A=((E-4)*E+4)/(8*H) B=(-E+2)/(4*H*H) C=1/(8*H**3) RETURN ELSE IF(J.GT.I-1.AND.J.GT.I) THEN A=((-E-4)*E-4)/(8*H) B=(E+2)/(4*H*H) C=-1/(8*H**3) RETURN ELSE IF(J.GT.I-1) THEN A=(3*E+4)*E/(8*H) B=(-3*E-2)/(4*H*H) C=3/(8*H**3) IF(I.NE.1 .AND.I.NE.N+1) RETURN ELSE A=(-3*E+4)*E/(8*H) B=(3*E-2)/(4*H*H) C=-3/(8*H**3) IF(I.NE.2 .AND.I.NE.N+2) RETURN END IF END IF END IF IF(I.LE.2) THEN AA=-1/(8*H) BB=1/(4*H*H) CC=-1/(8*H**3) IF(I.EQ.2) THEN A=A-AA B=B-BB C=C-CC RETURN ELSE A=A-4*AA B=B-4*BB C=C-4*CC RETURN END IF ELSE EE=N+2 AA=((EE-4)*EE+4)/(8*H) BB=(-EE+2)/(4*H*H) CC=1/(8*H**3) IF(I.EQ.N+1) THEN A=A-AA B=B-BB C=C-CC RETURN ELSE A=A-4*AA B=B-4*BB C=C-4*CC END IF END IF RETURN END C FUNCTION INTY(J,JJ) C C CORRESPONDS THE INTERVAL (X(JJ),X(JJ+1)) TO PROPER INDEX C K IN CO(J,K,L),L=1,...,4 AND IN DCO(J,K,L),L=1,2,3 IMPLICIT REAL*8(A-H,O-Z) IF(JJ.LE.J-3 .OR. JJ.GE.J+2) THEN WRITE(6,2) 2 FORMAT(1X,'ERROR IN INTY') STOP ELSE IF(JJ.EQ.J-2) INTY=1 IF(JJ.EQ.J-1) INTY=2 IF(JJ .EQ. J) INTY=3 IF(JJ.EQ.J+1) INTY=4 END IF RETURN END C FUNCTION P(X) IMPLICIT REAL*8(A-H,O-Z) P=1.0D+00 RETURN END C FUNCTION Q(X) IMPLICIT REAL*8(A-H,O-Z) Q=3.141593*3.141593 RETURN END C FUNCTION F(X) IMPLICIT REAL*8(A-H,O-Z) F=2*(3.141593)**2*DSIN(3.141593*X) RETURN END