CODEWITHSUNDEEP
fortrangfortran
Tack_Tau
Tack_Tau

Reputation: 657

Why midpoint rule turns out more accurate than Simpson's rule when doing riemann sum approximation on Fortran

everyone.

I am just playing with the calculation of integral of x^2 from [1, 2] using both midpoint rule and Simpson's rule. And I find it out that with the same number of subintervals midpoint rule approximation seems more accurate than Simpson's rule approximation, which is really weird.

The source code of midpoint rule approximation is : 

program midpoint
implicit none                 ! Turn off implicit typing
Integer, parameter :: n=100   ! Number of subintervals
integer :: i                  ! Loop index
real :: xlow=1.0, xhi=2.0     ! Bounds of integral
real :: dx                    ! Variable to hold width of subinterval
real :: sum                   ! Variable to hold sum
real :: xi                    ! Variable to hold location of ith subinterval
real :: fi                    ! Variable to value of function at ith subinterval
dx = (xhi-xlow)/(1.0*n)       ! Calculate with of subinterval
sum = 0.0                     ! Initialize sum
xi = xlow+0.5*dx              ! Initialize value of xi
do i = 1,n,1                  ! Initiate loop
 ! xi = xlow+(0.5+1.0*i)*dx
 write(*,*) "i,xi ",i,xi      ! Print intermidiate result
 fi = xi**2                   ! Evaluate function at ith point
 sum = sum+fi*dx              ! Accumulate sum
 xi = xi+dx                   ! Increment location of ith point
end do                        ! Terminate loop
write(*,*) "sum =",sum
stop                          ! Stop execution of the program
end program midpoint

the according execution is: 

  ......                   .....      ..................  
 i,xi          100   1.99499905    
 sum =   2.33332348

The source code of Simpson's rule approximation is: 

program simpson
implicit none                 ! Turn off implicit typing
integer, parameter :: n=100   ! Number of subintervals
integer :: i=0                ! Loop index
real :: xlow=1.0, xhi=2.0     ! Bounds of integral
real :: h                     ! Variable to hold width of subinterval
real :: sum                   ! Variable to hold sum
real :: xi                    ! Variable to hold location of ith subinterval
real :: fi                    ! Variable to value of function at ith subinterval
real :: Psimp                 ! Variable of simpson polynomial of xi interval
h = (xhi-xlow)/(1.0*n)        ! Calculate width of subinterval
sum = 0.0                     ! Initialize sum
do while (xi<=xhi-h)          ! Initiate loop
 xi = xlow+i*2.0*h            ! Increment of xi
 i=i+1
 write(*,*) "i,xi ",i,xi      ! Print intermidiate result
 Psimp=xi**2+4.0*(xi+h)**2+(xi+2.0*h)**2
                              ! Evaluate function at ith point
 sum = sum+(h/3.0)*Psimp      ! Accumulate sum
end do                        ! Terminate loop
write(*,*) "sum =",sum
end program simpson

the according execution is: 

 ........                  ......    ...................  
 i,xi          101   2.00000000    
 sum =   2.37353396

To get the same precision of digits as midpoint result, I have to set the number of subintervals in Simpson's program to 100000, which is 1000 times more than the midpoint program (I initially set both of the number subintervals to 100)

I check the codes in Simpson's program and can't find whats wrong.

Simpson's rule should converge more rapid than midpoint rule if I remembered it correct.

Upvotes: 0

Views: 1000

Answers (1)

user5713492
user5713492

Reputation: 974

Craig Burley once remarked that a WHILE loop looked like as soon as the premise of the loop was violated, the loop would be exited immediately. Here the premise of the loop is violated when x=xhi but the loop doesn't break at that point, only when a whole nother iteration is completed and the test can be applied at the top of the loop. You could more consistently with Fortran idioms convert the loop into a counted DO loop with something like 

DO i = 0, n/2-1

and then comment out the

  i=i+1

line. Or simply test the loop premise immediately after modifying xi:

xi = xlow+i*2.0*h ! Increment of xi
if(xi>xhi-h) exit ! Test loop premise

Either way leads to the exact results expected for a polynomial of degree no higher than 3 for Simpson's rule.

Upvotes: 1

Related Questions