Numerical Integration Using Fortran 90

Question

I'm trying to use Simpson's rule to evaluate the integral of sqrt(1-x^2) in the interval from -1 to 1. However, the sum represented by the variable "s", in the code that I developed, doesn't converge at all to pi over 2. I'm an absolute novice to fortran and programming in general, so please bear with me. What am I doing wrong?

PROGRAM integral

REAL*8 :: x
REAL*8 :: h
REAL*8 :: fodd
REAL*8 :: feven
REAL*8 :: simpson
REAL*8 :: s

x = -1
s = 0
simpson = 0
h = 0

   DO WHILE (x<=1)

     fodd = sqrt(1-(x+(2*h+0.1))**(2))
     feven = sqrt(1-(x+2*h)**(2))
     simpson = 4*fodd + 2*feven
     s = s + simpson*(h/3)
     WRITE(*,*) x,h, fodd, feven, simpson, s
     h = 0.1
     x = x + h
     END DO

END PROGRAM

Here's a link to the output it generates: https://pastebin.com/mW06Z6Lq

Since this integral is just half the area of a circle of radius 1 it should converge to pi over 2, but it surpasses this value by a large amount. I thought about making the step smaller for precision, but this is not the problem as it surpassed even more the expected value when I tried this.

Answer 1

You must check what happens when you get with x close to 1. You certainly cannot use DO WHILE (x<=1) because when x==1 then x+2*h is above 1 and you are computing the square root of a negative number.

I suggest not using DO WHILE at all. Just set the number of divisions, compute the step using the step size as the interval length / the number of divisions and then use

sum = 0
h = interval_length / n
x0 = -1

do i = 1, n
  xa = (i-1) * h + x0 !start of the subinterval
  xb = i * h + x0 !end of the subinterval    
  xab = (i-1) * h + h/2 + x0 !centre of the subinterval

  !compute the function values here, you can reuse f(xb) from the
  !last step as the current f(xa)


  !the integration formula you need comes here, adjust as needed
  sum = sum + fxa + 4 * fxab + fxb
end do

! final normalization, adjust to follow the integration formula above
sum = sum * h / 6

Note that the loop nest above is written in a very generic way, not specific to the Simpson's rule. I just assumed that h is constant, but even that can be chaged easily. For the Simpson's rule one can optimize it easily. You certainly want just two function evaluations per interval. If it is a school assignment and you need to treat the points as odd-even instead of centre-boundary, you must adjust that yourself, it is very easy.