I need help with Fortran 90 : Simpson's rule

fluidistic

I must calculate [tex]\int_0^1 e^{-x} dx[/tex] using the composite Simpson's rule, i.e. the common Simpson's rule but applied on many intervals between 0 and 1. This is not all : I must divide the interval [0,1] in 100 subintervals and then in 200, to compare the value obtained of the integral. And then I must calculate the coefficient of precision, that is Q=absolute value of [tex]\frac{s-I}{r-I}[/tex], where I is the (correct) value of the integral we are searching, s its approximation by using the Simpson's rule on 100 subintervals and r its approximation by using the S' rule on 200 subintervals. And the answer, Q, is 16. While I find almost 2. I've thought a lot about my error, and I know it's in the subroutine when computing the integral, but can't figure where exactly. The bad thing is that my result of the integral is quite similar to the real one... Should be a little mistake then.
I forgot to say, my program may be a little hard to understand. When you execute it, first it will ask the bounds of the integral, just put 0,1. Then chose 100 for n and 200 for m and we're done.
Program Simpson
implicit none

Real(8) :: x,s,x_i,x_f,r,q
Integer :: k,n,m

Write(*,*)'Enter x_i y x_f'
Read(*,*)x_i,x_f
Write(*,*)'Enter n'
Read(*,*)n
Write(*,*)'Enter m'
Read(*,*)m

Call Simp(x_i,x_f,s,n)
Write(*,*)'When n=100, s is worth',s

Call Simp2(x_i,x_f,r,m)
Write(*,*)'When n=200, s is worth',r


Q=(s-(1-exp(-1.)))/(r-(1-exp(-1.)))
Write(*,*)'The quotient of precision is',Q



Contains
Subroutine Simp(x_i,x_f,s,n)
implicit none
Real(8), intent(in) :: x_i,x_f
Real(8), Intent(out) :: s
Real(8) :: dx
Integer :: n

dx=(x_f-x_i)/n
s=(dx/3.)/(f(x_i)+f(x_f))

Do k=2,n-2,2
x=x_i+k*dx
s=s+2*f(x)
end do

do k=1,n-1,2
x=x_i+k*dx
s=s+4*f(x)
end do
s=(dx/3.)*s

end subroutine

Subroutine Simp2(x_i,x_f,r,m)
implicit none
Real(8), Intent(in) :: x_i,x_f
Real(8), Intent(out) :: r
Real(8) :: dx
Integer :: m

dx=(x_f-x_i)/m
r=(dx/3.)/(f(x_i)+f(x_f))

Do k=2,m-2,2
x=x_i+k*dx
r=r+2*f(x)
end do

do k=1,m-1,2
x=x_i+k*dx
r=r+4*f(x)
end do
r=(dx/3.)*r
end subroutine


Real(8) Function f(x)
Real(8), Intent(in) :: x
f=exp(-x)
end function
end program

stingray78

try calculationg a really simple integral to test your program. If that works it's probably an error of the machine. Try placing instead of real integers like 4 use 4.0.
Also try placing prints along the program to see if everything is working ok.
If that doesnt work try this other simpson rule:
http://www.tubepolis.com/play.php?q=...52Fdefault.jpg
Jejejeje maybe that works

fluidistic

All works fine, I tested it under both linux and windows, with and without the double precision of digits and on 2 different machines. Still don't know where my subroutines are wrong with the Simpson's rule.

alphysicist

Hi fluidistic,

In the second executable line of your subroutine simp and simp2 you have:

I believe this should just be

I definitely don't think the endpoint values should be in a denominator; and also since you multiply the entire s value by (dx/3.) at the last statement of the subroutine, having it here also in this statement makes it a factor (dx/3.)^2 for the endpoints.

ice109

here's my simpson's rule integration subroutine, the (prec) has to do with the precision, just remove it.


subroutine simpson ( a,b,n_points, approx ) ! simpsons rule integration

real(prec), intent(in) :: a,b
integer, intent(in) :: n_points
integer :: k,i
real(prec), intent(out) :: approx
real(prec) :: del_x, height, area, t,x
real(prec) :: n_points_1
n_points_1 = real(n_points,DP)

del_x = abs(b-a)/n_points_1


approx = 0
t = a

do k=1, n_points-1

x = t + del_x/two
area = del_x/six * ( four*functio_n ( x ) + functio_n(t) + functio_n(t+del_x) )
t = t + del_x

approx = approx + area

end do

end subroutine

fluidistic

Oh... you're right alphysicist! But even changed that I don't get a quotient of 16, but 1 now.

fluidistic

Thanks for your answer ice109. But I don't understand all your program. I'm quite new with fortran. For example, when you first define " n_points_1 = real(n_points,DP)" real() means something special? And what is DP?
I tried without the precision once again in my program and it's the same, I don't get the 16 value I should.

ice109

use mine and declare functio_n to be what ever you want.

alphysicist

fluidistic,

You're losing your precision when you calculate Q. You have:

Q=(s-(1-exp(-1.)))/(r-(1-exp(-1.)))

which mixes type 8 real variables ( s and r) with an integer constant (1) and a default type real constant ( -1.). You need to specify that your numeric constants are type 8 reals by using:

Q=(s-(1._8-exp(-1._8)))/(r-(1._8-exp(-1._8)))

fluidistic

Alphysicist I can't believe it! It works now... I would never have thought about that. Thanks a lot.
Ice109, if you're still interested in offering me your program, you should give it to me entirely otherwise I would have to define not only functio_n but everything, and as I said, I don't understand this part : "n_points_1 = real(n_points,DP)". Thanks to all of you that pay some attention and time to my problem.

ice109

change that to simply real(n_points). that's the only other thing that's not portable.

alphysicist

Hi ice109,

I think there is a problem in one line of your subroutine. Where you have the line:

del_x = abs(b-a)/n_points_1

I think this needs to be

del_x = abs(b-a)/(n_points_1-1.)

(where the 1. needs to match the precision of the variables.) For example, if you were choosing 3 points, you would divide the interval (b-a) by 2 to find the stepsize.


When I have functio_n(x)=1, for example, for the integral from 0 to 1 with n=50, I get a result 0.9800000000000005 with a 2% error; if I make that change in the denominator, I get the result 1.0000000000000007.

fluidistic

Nicely seen Alphysicist. You're right, using the program of ice109 I got a close result, but changing it in your way I obtain exactly what I obtain using my program.
I post his correct program here :

Program Simpsons
implicit none

real, Parameter :: a=0,b=1
integer :: n_points
integer :: k,i
real :: approx
real :: del_x, height, area, t,x
real :: n_points_1

Write(*,*)'Enter number of points'
Read(*,*)n_points

n_points_1 = real(n_points)

del_x = abs(b-a)/(n_points_1-1.)


approx = 0
t = a

do k=1, n_points-1

x = t + del_x/2
area = del_x/6 * ( 4*functio_n ( x ) + functio_n(t) + functio_n(t+del_x) )
t = t + del_x

approx = approx + area

end do



Write(*,*)approx


Contains
Real Function functio_n(x)
Real, Intent(in) :: x
functio_n=exp(-x)
end function

end program

surfernj

I tried running this as it's written using the Lahey compilier and i get a load of errors.
What am I doing wrong?

fluidistic

Which program doesn't work for you?
If it's the first or the last, then it's strange. Maybe try to compile under gfortran, the program I'm using for it.

surfernj

the first program...i get errors along the lines of.."Contains" requires stop or return....also errors that state n & m must be used with "intent". In total its about 17 errors.
I don't have the compiler on my home pc or i'd post the actual error messages for you.

fluidistic

Sadly I'm not a specialist, so I cannot help you... I can try post my program on an attached file.
It's almost the same I posted here at first (but in Spanish and corrected). Download it and tell us what happens.

Attached Files

Simpson.txt (1.0 KB, 12 views)