I with Fortran 90 : Simpson's rule

  • Fortran
  • Thread starter fluidistic
  • Start date
  • Tags
    Fortran
In summary, I must calculate the integral \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 sub
  • #1
fluidistic
Gold Member
3,923
260
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
 
Technology news on Phys.org
  • #2
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 doesn't 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:
  • #3
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
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.
 
  • #5
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
Oh... you're right alphysicist! But even changed that I don't get a quotient of 16, but 1 now.
 
  • #7
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
use mine and declare functio_n to be what ever you want.
 
  • #9
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
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
fluidistic said:
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.
 
  • #12
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
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
 
  • #14
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
surfernj said:
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.
 
  • #16
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
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: 562
  • #18
thanks i'll try it...i don't have class again till monday so i'll post my results then.
 
  • #19
Hi surfernj,

surfernj said:
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.
 
  • #20
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'
24 Read(UNIT=*,FMT=*)x_i,x_f
25 Write(UNIT=*,FMT=*)'Enter n'
26 Read(UNIT=*,FMT=*)n
27 Write(UNIT=*,FMT=*)'Enter m'
28 Read(UNIT=*,FMT=*)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:
  • #21
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.
 
  • #22
from 0 to 2
 
  • #23
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.
 
  • #24
You're correct! I don't know what I was looking at but the area is correct.
 
  • #25
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!
 
  • #26
In the program it asks an input of the number of points which an even number can be enered. surly only an odd number of points can be entered as composite simpsons only works with even intervals. I am just a little confused. Any help would be greatly appreciated.
 
  • #27
I didn't quite follow the whole thread, but would you try using the next odd number, i.e. use n=n/2*2+1?
 
  • #28
bigev234 said:
In the program it asks an input of the number of points which an even number can be enered. surly only an odd number of points can be entered as composite simpsons only works with even intervals. I am just a little confused. Any help would be greatly appreciated.

Are you talking about the program in the original post? It had two errors (a formula error, and a precision error), but it looks to me like the inputs n and m are the number of subintervals, and so should be even numbers.
 
  • #29
The composite simpson's by fluidistic on the first page. So should an odd or even value be entered in THE NUMBER OF POINTS. pLEASE HELP!
 
  • #30
bigev234 said:
The composite simpson's by fluidistic on the first page. So should an odd or even value be entered in THE NUMBER OF POINTS. pLEASE HELP!

From your question, I think you might be misunderstanding the program. In fluidisitic's program, you don't enter the number of points, you enter the number of subintervals. So when it asks for n or m, if you enter an even number, then that would correspond to an odd number of points.

So when n=100, that corresponds to 101 points and 100 intervals.
 
  • #31
Is the program in this case a good one to use for simpson's composite rule because it gives a value with a very small difference when entering an ood number compared to an even number when simpson's composite needs only an even number of intervals
 
  • #32
bigev234 said:
Is the program in this case a good one to use for simpson's composite rule because it gives a value with a very small difference when entering an ood number compared to an even number when simpson's composite needs only an even number of intervals

I think I would have to disagree with that. The difference is not small at all.

The value of the integral (integrating exp(-x) from 0 to 1) should be 0.63212 to five digits.

When using n=100 intervals, you get that value; if you use n=101 intervals, you get 0.628448. That's about a half a percent difference, which I would say is a huge error for this integral.


For that matter, you can just replace the function given in the program with 1; that is, calculate the integral of 1 from 0 to 1. Of course it should be give 1, right? But the results are:

n=100 ===> integral=1.
n=101 ===> integral=0.990099


The error formula predicts zero error for this type of integral. But for an odd number of intervals it is very close to 1% error--much too large for the integral of a constant. I would say you need to use an even number of subintervals.
 
  • #33
hellloo..how can i make my compoosite simpson fortran code general.i mean...take the function as user's input and solve according..heres the codde..what needs to be changed?.

PROGRAM SIMPSON
C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994
C To accompany the text:
C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
C This free software is complements of the author.
C
C Algorithm 7.2 (Composite Simpson Rule).
C Section 7.2, Composite Trapezoidal and Simpson's Rule, Page 365
C
INTEGER M
REAL A,B,Srule
CHARACTER*60 ANS*1,DFUN,FUN
EXTERNAL F
10 CALL INPUTS(A,B,M,DFUN,FUN)
CALL SIMPRU(F,A,B,M,Srule)
CALL RESULT(A,B,M,Srule,DFUN,FUN)
WRITE(9,*)'WANT TO TRY ANOTHER INTERVAL ? <Y/N> '
READ (9,'(A)') ANS
IF (ANS.EQ.'Y' .OR. ANS.EQ.'y') GOTO 10
END

REAL FUNCTION F(X)
REAL X
F=X/(1+X*X)
RETURN
END

SUBROUTINE PRINTFUN(DFUN,FUN)
CHARACTER*(*) DFUN,FUN
FUN ='X/(1+X*X)'
DFUN='X/(1+X*X) DX'
RETURN
END

SUBROUTINE SIMPRU(F,A,B,M,Srule)
INTEGER K,M
REAL A,B,H,Sum,SumEven,SumOdd,Srule,X
EXTERNAL F
H=(B-A)/(2*M)
SumEven=0
DO K=1,(M-1)
X=A+H*2*K
SumEven=SumEven+F(X)
ENDDO
SumOdd=0
DO K=1,M
X=A+H*(2*K-1)
SumOdd=SumOdd+F(X)
ENDDO
Sum=H*(F(A)+F(B)+2*SumEven+4*SumOdd)/3
Srule=Sum
RETURN
END

SUBROUTINE XSIMPRU(F,A,B,M,Srule)
C This subroutine uses labeled DO loop(s).
INTEGER K,M
REAL A,B,H,Sum,SumEven,SumOdd,Srule,X
EXTERNAL F
H=(B-A)/(2*M)
SumEven=0
DO 10 K=1,(M-1)
X=A+H*2*K
SumEven=SumEven+F(X)
10 CONTINUE
SumOdd=0
DO 20 K=1,M
X=A+H*(2*K-1)
SumOdd=SumOdd+F(X)
20 CONTINUE
Sum=H*(F(A)+F(B)+2*SumEven+4*SumOdd)/3
Srule=Sum
RETURN
END

SUBROUTINE INPUTS(A,B,M,DFUN,FUN)
INTEGER I,M
REAL A,B
CHARACTER*(*) DFUN,FUN
CALL PRINTFUN(DFUN,FUN)
DO 10 I=1,18
WRITE(9,*)' '
10 CONTINUE
WRITE(9,*)' SIMPSON`S RULE IS USED TO COMPUTE AN APPROXIMATION'
WRITE(9,*)' '
WRITE(9,*)'FOR THE VALUE OF THE DEFINITE INTEGRAL:'
WRITE(9,*)' '
WRITE(9,*)' '
WRITE(9,*)' B'
WRITE(9,*)' /'
WRITE(9,*)' | ',DFUN
WRITE(9,*)' /'
WRITE(9,*)' A'
WRITE(9,*)' '
WRITE(9,*)' '
WRITE(9,*)'ENTER THE LEFT ENDPOINT A = '
READ(9,*) A
WRITE(9,*)' '
WRITE(9,*)'ENTER THE RIGHT ENDPOINT B = '
READ(9,*) B
WRITE(9,*)' '
WRITE(9,*)'THE NUMBER OF SUBINTERVALS USED IS 2*M'
WRITE(9,*)' ENTER THE NUMBER M = '
READ(9,*) M
WRITE(9,*)' '
RETURN
END

SUBROUTINE RESULT(A,B,M,Srule,DFUN,FUN)
INTEGER I,M
REAL A,B,Srule
CHARACTER*(*) DFUN,FUN
CALL PRINTFUN(DFUN,FUN)
DO 10 I=1,18
WRITE(9,*)' '
10 CONTINUE
WRITE(9,*)' ',B
WRITE(9,*)' /'
WRITE(9,*)' |'
WRITE(9,*)Srule,' ~ | ',DFUN
WRITE(9,*)' |'
WRITE(9,*)' /'
WRITE(9,*)' ',A
WRITE(9,*)' '
WRITE(9,*)' '
WRITE(9,*)'AN APPROXIMATE VALUE FOR THE DEFINITE INTEGRAL OF'
WRITE(9,*)' '
WRITE(9,*)' '
WRITE(9,*)'F(X) = ',FUN
WRITE(9,*)' '
WRITE(9,*)' '
WRITE(9,*)'TAKEN OVER [',A,' ,',B,' ] WAS FOUND.'
WRITE(9,*)' '
WRITE(9,*)'SINCE M = ',M,', THERE WERE ',2*M,' SUBINTERVALS.'
WRITE(9,*)' '
WRITE(9,*)'THE SIMPSON RULE APPROXIMATION IS ',Srule
WRITE(9,*)' '
RETURN
END
 

1. What is Simpson's rule and why is it important in Fortran 90?

Simpson's rule is a numerical method used for approximating the value of a definite integral. It is important in Fortran 90 because it allows for efficient and accurate calculation of integrals, which is essential in many scientific and engineering applications.

2. How does Simpson's rule work in Fortran 90?

In Fortran 90, Simpson's rule involves dividing the interval of integration into smaller subintervals and using a quadratic polynomial to approximate the function within each subinterval. The values of the function at the endpoints and midpoint of each subinterval are then used to calculate the integral.

3. What are the advantages of using Simpson's rule in Fortran 90?

One advantage of Simpson's rule in Fortran 90 is that it typically provides more accurate results compared to other numerical integration methods, especially for functions that are smooth and have no discontinuities. Additionally, it is relatively easy to implement in Fortran 90 and can handle a wide range of integrals.

4. Are there any limitations to using Simpson's rule in Fortran 90?

One limitation of Simpson's rule in Fortran 90 is that it may not be as accurate for functions with sharp corners or discontinuities. In these cases, other numerical integration methods may be more suitable. Additionally, the accuracy of Simpson's rule can be affected by the number of subintervals used and the spacing between them.

5. Can Simpson's rule be used for multidimensional integrals in Fortran 90?

Yes, Simpson's rule can be extended to handle multidimensional integrals in Fortran 90 by using a similar approach to the one-dimensional case. The interval of integration is divided into smaller subintervals in each dimension, and a polynomial is used to approximate the function within each subinterval. The values at the endpoints and midpoints of each subinterval are then used to calculate the integral in each dimension, and the results are combined to obtain the final approximation.

Similar threads

  • Programming and Computer Science
Replies
4
Views
499
  • Programming and Computer Science
Replies
8
Views
1K
  • Programming and Computer Science
Replies
1
Views
916
  • Programming and Computer Science
Replies
4
Views
1K
Replies
3
Views
1K
  • Programming and Computer Science
Replies
12
Views
931
  • Programming and Computer Science
Replies
8
Views
1K
Replies
16
Views
2K
  • Programming and Computer Science
Replies
20
Views
1K
  • Programming and Computer Science
Replies
11
Views
1K
Back
Top