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

by fluidistic
Tags: fortran, rule, simpson
 PF Gold P: 3,207 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' 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
 P: 23 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
 PF Gold P: 3,207 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.
 HW Helper P: 2,249 I need help with Fortran 90 : Simpson's rule Hi fluidistic, In the second executable line of your subroutine simp and simp2 you have: s=(dx/3.)/(f(x_i)+f(x_f)) I believe this should just be 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.
 P: 1,705 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
 PF Gold P: 3,207 Oh... you're right alphysicist! But even changed that I don't get a quotient of 16, but 1 now.
 PF Gold P: 3,207 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.
 P: 1,705 use mine and declare functio_n to be what ever you want.
 HW Helper P: 2,249 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)))
 PF Gold P: 3,207 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.
P: 1,705
 Quote by 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.
change that to simply real(n_points). that's the only other thing that's not portable.
 HW Helper P: 2,249 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.
 PF Gold P: 3,207 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
 P: 7 I tried running this as it's written using the Lahey compilier and i get a load of errors. What am I doing wrong?
PF Gold
P: 3,207
 Quote by 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?
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.
 P: 7 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.
PF Gold
P: 3,207
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, 16 views)
 P: 7 thanks i'll try it...i don't have class again till monday so i'll post my results then.

 Related Discussions Introductory Physics Homework 4 Calculus & Beyond Homework 7 Calculus & Beyond Homework 1 Calculus & Beyond Homework 2 Computing & Technology 2