1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
|
module proc
implicit none
real :: x,y,z,a
contains
!! Definition des fonctions
function func(a)
implicit none
real, intent(in) :: a
real :: func, aire1
call romberg(func1,0.0,1.0,0.001,aire1)
func = aire1
end function func
function func1(z)
implicit none
real, intent(in):: z
real :: func1, aire2
call romberg(func2,0.0,z,0.001,aire2)
func1 = aire2
end function func1
function func2(y)
implicit none
real, intent(in):: y
real :: func2, aire3
call romberg(func3,0.0,y,0.001,aire3)
func2 = aire3
end function func2
function func3(x)
implicit none
real, intent(in):: x
real :: func3
!! Test function of 3 variables
func3 = a*(a**2*x+0.44*y+2.0*z)
end function func3
subroutine romberg(fn,a,b,eps,aire)
!! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !!
!! !!
!! Integration par la methode de Romberg !!
!! Dominique Lefebvre fevrier 2007 !!
!! Inspire du code de A.Garcia dans Numerical Methods for Physics !!
!! Modifiee le 30/08/09 suite a une erreur d'estimation de la precision !!
!! signalee par Maelle Nodet (Universite de Grenoble) !!
!! a = borne inferieure d'integration !!
!! b = borne superieure d'integration !!
!! eps = precision souhaitee !!
!! aire = surface retournee !!
!! !!
!! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - !!
implicit none
real , intent(in) :: a, b, eps
real , intent(out):: aire
real :: h,xi, delta, sum
real :: c, fotom, t(136)
integer :: out, pieces, i, nx(16)
integer :: nt, ii, n, nn, l, ntra, k, m, j, l2
logical :: fini
!! Interface ajoutee par saidou Diallo
interface
function fn(x)
implicit none
real, intent(in) :: x
real :: fn
end function fn
end interface
!! calcul du premier terme
h = b-a
pieces = 1
nx(1) = 1
delta = h / pieces
c = (fn(a) + fn(b)) / 2.0
sum = c
t(1) = delta * c
n = 1
nn = 2
fini =.false.
do while (.not. fini)
n = n + 1
fotom = 4.0
nx(n) = nn
pieces = pieces * 2
l = pieces - 1
l2 = (l+1) / 2
delta = h / pieces
!! evaluation par la methode des trapezes
do ii = 1, l2
i = ii * 2 - 1
xi = a + delta * i
sum = sum + fn(xi)
enddo
t(nn) = sum * delta
ntra = nx(n-1)
k = n-1
!! calcul de la nieme colonne du tableau
do m = 1, k
j = nn + m
nt = nx(n-1) + m - 1
t(j) = (fotom * t(j-1) - t(nt)) / (fotom - 1)
fotom = fotom * 4
enddo
!! estimation de la precision
if (n < 5) then
nn = j+1
else if (t(nn+1) == 0.0) then
nn = j+1
else
if (abs(t(ntra+1)-t(nn+1)) <= abs(t(nn+1)*eps)) then
fini = .true.
elseif (abs(t(nn-1) - t(j)) <= abs(t(j)*eps)) then
fini = .true.
else
nn = j+1
endif
endif
enddo
!! on retourne la valeur finale de l'aire
aire = t(j)
return
end subroutine romberg
end module proc
!! Main program
program test
use proc
implicit none
a = 1.0
print*, func(a)
end program test |
