PROGRAM NumInt ********************************************************************** * Berechnung eines eindimensionalen Integrals mittels * Trapezverfahrens und Monte Carlo Technik ********************************************************************* REAL R(97) EXTERNAL RAN F(X) = 6.*X**2-10.*X+12. STF(X) = 2.*X**3-5*X**2+12.*X PRINT *, & 'Numerische Berechnung der Integrals ueber die Funktion' PRINT *, & 'F(X) = 6.*X**2-10.*X+12. von A=2 bis B=12' PRINT *, & 'Eingabe Anzahl der Stützpunkte und ISEED fuer den', & ' Zufallszahlengenerator:' READ *, M,ISEED A=2 B=12 RINT=0. C RAN Zufallszahl zwischen 0 und 1 DO 100, I=1,M X=A+RAN(ISEED,IX1,IX2,IX3,R)*(B-A) RINT=RINT+F(X) 100 CONTINUE RINT=(B-A)*RINT/REAL(M) C Trapezregel TINT = 0. TINT = TINT + F(A) + F(B) IF (M.GT.0) THEN DO 200, I=1,M-1 X=A+ REAL(I)*(B-A)/REAL(M) TINT=TINT+2.*F(X) 200 CONTINUE TINT = 0.5*TINT*(B-A)/REAL(M) EXAKT=STF(B)-STF(A) PRINT *, 'Wert des Integrals nach der Trapezregel und ', & 'dem Monte Carlo Verfahren:' PRINT *, TINT,RINT ELSE PRINT *, 'Anzahl der Stuetzpunkte muss groesser als Null sein' ENDIF END * * REAL FUNCTION RAN(IDUM,IX1,IX2,IX3,R) ********************************************************************** * Returns a random deviate between 0.0 and 1.0. Set IDUM to * any negative value to initialize or reinitialize the sequence. * Numerical Recipes p.196 ********************************************************************** REAL R(97) INTEGER IX1,IX2,IX3,IC1,IC2,IC3,IA1,IA2,IA3,IFF PARAMETER (M1=259200,IA1=7141,IC1=54733,RM1=1./M1) PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=1./M2) PARAMETER (M3=243000,IA3=4561,IC3=51349) DATA IFF /0/ IF (IDUM.LT.0.OR.IFF.EQ.0) THEN IFF=1 IX1=MOD(IC1-IDUM,M1) IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IX1,M2) IX1=MOD(IA1*IX1+IC1,M1) IX3=MOD(IX1,M3) DO 11, J=1,97 IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IA2*IX2+IC2,M2) R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1 11 CONTINUE IDUM=1 ENDIF IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IA2*IX2+IC2,M2) IX3=MOD(IA3*IX3+IC3,M3) J=1+(97*IX3)/M3 * IF (J.GT.97.OR.J.LT.1) PAUSE RAN=R(J) R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1 RETURN END * *