龙贝格数值积分的fortran源程序

风花雪月

1.       SUBROUTINE ROMINT ( VAL, ERR, EPS, A, B, N, MAXE, F )
   
2. C
   
3.       IMPLICIT NONE
   
4. 
   
5.       REAL A
   
6.       REAL B
   
7.       REAL BB
   
8.       REAL EPS
   
9.       REAL ERR
   
10.       EXTERNAL F
    
11.       REAL F
    
12.       REAL H
    
13.       INTEGER J
    
14.       INTEGER K
    
15.       INTEGER K0
    
16.       INTEGER K1
    
17.       INTEGER K2
    
18.       INTEGER KK
    
19.       INTEGER KKK
    
20.       INTEGER L
    
21.       INTEGER MAXE
    
22.       INTEGER N
    
23.       INTEGER N0
    
24.       INTEGER N1
    
25.       REAL R
    
26.       REAL RM(16)
    
27.       REAL S
    
28.       REAL S0
    
29.       REAL S1
    
30.       REAL T
    
31.       REAL VAL
    
32. C
    
33. C  INITIAL TRAPEZOID RULE.
    
34. C
    
35.       T = ( B - A ) * ( F(A) + F(B) ) * 0.5E+00
    
36. C
    
37. C  INITIAL RECTANGLE VALUE.
    
38. C
    
39.       RM(1) = ( B - A ) * F ( ( A + B ) * 0.5E+00 )
    
40. 
    
41.       N = 2
    
42.       R = 4.0E+00
    
43. 
    
44.       DO 11 K = 1, MAXE
    
45. 
    
46.         BB = ( R * 0.5E+00 - 1.0E+00 ) / ( R - 1.0E+00 )
    
47. C
    
48. C  IMPROVED TRAPEZOID VALUE.
    
49. C
    
50.         T = RM(1) + BB * ( T - RM(1) )
    
51. C
    
52. C  DOUBLE NUMBER OF SUBDIVISIONS OF (A,B).
    
53. C
    
54.         N = 2 * N
    
55.         S = 0.0E+00
    
56.         H = ( B - A ) / FLOAT ( N )
    
57. C
    
58. C  CALCULATE RECTANGLE VALUE.
    
59. C
    
60.         IF ( N - 32 ) 1, 1, 2
    
61. 1       N0 = N
    
62.         GO TO 4
    
63. 2       N0 = 32
    
64. 3       IF ( N - 512 ) 4, 4, 5
    
65. 4       N1 = N
    
66.         GO TO 6
    
67. 5       N1 = 512
    
68. 6       DO 9 K2 = 1, N, 512
    
69.           S1 = 0.0E+00
    
70.           KK = K2 + N1 - 1
    
71.           DO 8 K1 = K2, KK, 32
    
72.             S0 = 0.0E+00
    
73.             KKK = K1 + N0 - 1
    
74.             DO 7 K0 = K1, KKK, 2
    
75.               S0 = S0 + F ( A + FLOAT ( K0 ) * H )
    
76. 7           CONTINUE
    
77.             S1 = S1 + S0
    
78. 8         CONTINUE
    
79.           S = S + S1
    
80. 9       CONTINUE
    
81.         RM(K+1) = 2.0E+00 * H * S
    
82. C
    
83. C  END CALCULATION OF RECTANGLE VALUE.
    
84. C
    
85.         R = 4.0E+00
    
86. C
    
87. C  FORM ROMBERG TABLE FROM RECTANGLE VALUES.
    
88. C
    
89.         DO 10 J = 1, K
    
90.           L = K + 1 - J
    
91.           RM(L) = RM(L+1) + ( RM(L+1) - RM(L) ) / ( R - 1.0E+00 )
    
92.           R = 4.0E+00 * R
    
93. 10      CONTINUE
    
94. 
    
95.         ERR = ABS ( T - RM(1) ) * 0.5E+00
    
96. C
    
97. C  CONVERGENCE TEST.
    
98. C
    
99.         IF ( ERR - EPS ) 12, 12, 11
    
100. 
     
101. 11    CONTINUE
     
102. 
     
103.       K = 0
     
104. 
     
105. 12    VAL = ( T + RM(1) ) * 0.5E+00
     
106.       N = N + 1
     
107.       MAXE = K
     
108. 
     
109.       RETURN
     
110.       END

复制代码