# I need help with Fortran 90 : Simpson's rule

1. Apr 29, 2008

### fluidistic

I must calculate $$\int_0^1 e^{-x} dx$$ 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 $$\frac{s-I}{r-I}$$, 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'
Write(*,*)'Enter n'
Write(*,*)'Enter 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

2. Apr 29, 2008

### stingray78

3. Apr 29, 2008

### 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.

4. Apr 29, 2008

### alphysicist

Hi fluidistic,

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

Code (Text):

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

I believe this should just be

Code (Text):

s=(f(x_i)+f(x_f))

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.

5. Apr 29, 2008

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

6. Apr 29, 2008

### fluidistic

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

7. Apr 29, 2008

### 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.

8. Apr 29, 2008

### ice109

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

9. Apr 29, 2008

### 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)))

10. Apr 29, 2008

### 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.

11. Apr 29, 2008

### ice109

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

12. Apr 29, 2008

### 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.

13. Apr 30, 2008

### 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'

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

14. May 8, 2008

### surfernj

can't get this to work

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

15. May 8, 2008

### 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.

16. May 8, 2008

### 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.

17. May 8, 2008

### 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
File size:
1 KB
Views:
95
18. May 8, 2008

### surfernj

thanks i'll try it...i don't have class again till monday so i'll post my results then.

19. May 9, 2008

### alphysicist

Hi surfernj,

I think the Lahey compiler is strict on what it allows. I don't have access to a Lahey compiler but you might try the following things to correct the errors messages you mentioned:

put the statement

RETURN

before every END FUNCTION and END SUBROUTINE (3 places)

Also put:

STOP

before the line with CONTAINS

In the subroutines, you'll find two lines that say:

Integer::m

and

Integer::n

Change these to:

Integer,Intent(in)::m

and

Integer,Intent(in)::n

One of my compilers did not accept kind=8 reals; I did a find and replace to convert every 8 to a 2 and it compiled fine.

20. May 9, 2008

### surfernj

ok its running now however i'm trying to find the area from the function f(x)= e^(x/2)^2 from 0 to 2
and i keep getting 2.9.... using 20 or 50 divisions when the answer should be 4.7....

heres the code:
11 IMPLICIT NONE
12 INTEGER, PARAMETER:: DOUBLE=kind(1.0D0)
13 REAL(KIND=DOUBLE) :: e,x,s,x_i,x_f,r,q
14 Integer :: k,n,m
15 !
16 ! Open Printer
17 OPEN (UNIT= 2, FILE='iii1')
18 !
19 ! Get value of e
20 e = EXP(1.0)
21 !
22 !
23 Write(UNIT=*,FMT=*)'Enter x_i y x_f'
25 Write(UNIT=*,FMT=*)'Enter n'
27 Write(UNIT=*,FMT=*)'Enter m'
29 !
30 Call Simp(x_i,x_f,s,n)
31 Write(UNIT=*,FMT=*)'When n=20, s is worth',s
32 !
33 Call Simp2(x_i,x_f,r,m)
34 Write(UNIT=*,FMT=*)'When n=50, s is worth',r
35 !
36 !
37 Q=(s-(1-exp(-1.)))/(r-(1-exp(-1.)))
38 Write(UNIT=*,FMT=*)'The quotient of precision is',Q
39 !
40 !
41 stop
42 Contains
43 Subroutine Simp(x_i,x_f,s,n)
44 implicit none
45 Real(KIND=DOUBLE), intent(in) :: x_i,x_f
46 Real(KIND=DOUBLE), Intent(out) :: s
47 Real(KIND=DOUBLE) :: dx
48 INTEGER, INTENT(IN) :: n
49 !
50 !
51 dx=(x_f-x_i)/n
52 s=(f(x_i)+f(x_f))
53 !
54 !
55 Do k=2,n-2,2
56 x=x_i+k*dx
57 s=s+2*f(x)
58 end do
59 !
60 !
61 Do k=1,n-1,2
62 x=x_i+k*dx
63 s=s+4*f(x)
64 END DO
65 !
66 s=(dx/3.)*s
67 !
68 return
69 END Subroutine Simp
70 !
71 !
72 Subroutine Simp2(x_i,x_f,r,m)
73 implicit none
74 Real(KIND=DOUBLE), Intent(in) :: x_i,x_f
75 Real(KIND=DOUBLE), Intent(out) :: r
76 Real(KIND=DOUBLE) :: dx
77 INTEGER, INTENT(IN) :: m
78 !
79 !
80 dx=(x_f-x_i)/m
81 r=(f(x_i)+f(x_f))
82 !
83 !
84 DO k=2,m-2,2
85 x=x_i+k*dx
86 r=r+2*f(x)
87 END Do
88 !
89 !
90 DO k=1,m-1,2
91 x=x_i+k*dx
92 r=r+4*f(x)
93 END DO
94 r=(dx/3.)*r
95 return
96 end Subroutine Simp2
97 !
98 ! Defining the internal function
99 FUNCTION F(x)
100 REAL (KIND = double) :: F
101 REAL (KIND = double), INTENT(IN) :: X
102 ! Defining the equation
103 f= e**((x/2)**2)
104 return
105 END function F

Last edited: May 9, 2008