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 |
Partager