I need help with Fortran 90 : Simpson's rule

• Fortran
Gold Member
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

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=jessica%20simpson%20hot&title=Jessica%2BSimpson%2BSexy%2BHot%2BPhoto%2BGallery&engine=1&id=eRakvaM_ye8&img=http%253A%252F%252Fi.ytimg.com%252Fvi%252FeRakvaM_ye8%252Fdefault.jpg [Broken]
Jejejeje maybe that works

Last edited by a moderator:
Gold Member
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
Homework Helper
Hi fluidistic,

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

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

I believe this should just be

Code:
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.

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

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

Gold Member
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.

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

alphysicist
Homework Helper
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)))

Gold Member
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.

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.

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

alphysicist
Homework Helper
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.

Gold Member
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

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?

Gold Member
I tried running this as it's written using the Lahey compilier and i get a load of errors.
What am I doing wrong?
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.

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.

Gold Member
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.

Attachments

• Simpson.txt
1 KB · Views: 502
thanks i'll try it...i don't have class again till monday so i'll post my results then.

alphysicist
Homework Helper
Hi 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.

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.

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:
Gold Member
ok its running now however i'm trying to find the area from the function f(x)= e^(x/2)^2
and i keep getting 2.9.... using 20 or 50 divisions when the answer should be 4.7....
What bounds do you take for that? I mean you are approximating a definite integral so it has 2 bounds. It's maybe an approximation at a pretty "bad" part of the function, this way it should explain why you don't get a nice approximation. Search on wikipedia the error of the composite Simpson's rule, it will tell you exactly the error committed and I hope it fits your result.

from 0 to 2

Gold Member
Well, according to Mathematica, the first 200 digits of the exact value of the integral you gave evaluated from 0 to 2 are : 2.92530349181436321760809717371397631024174019243347837132022916043752
6628581958341643799625435070321638686230019422651851777782062594018065
1700237522944901229932963456333375795347145224897561966047484
Thus your approximation is a good one.

You're correct! I don't know what I was looking at but the area is correct.

I just realized what the problem is. I coded the equation in wrong! I coded f= e**((x/2)**2)
when the equation i'm supposed to code in is f=e**(x**2/2), from 0 to 2 the area is 4.72890778561. LOL I can't believe I made such a stupid mistake!