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.

You should certainly set `h` before the loop. It is not the main problem here, but it is important for the future. Also, you really should tell us what exact results you are getting and which exact results you are expecting intead. — Vladimir F Героям слава, Apr 06 '19 at 20:52
When we speak about convergence in numerical integration, we usually mean the convergence for h->0. You only have h=0.1 now, you should go to lower and lower values. Also, since you claim to be a beginner, note thatr `real*8` is not standard Fortran https://stackoverflow.com/questions/838310/fortran-90-kind-parameter https://stackoverflow.com/questions/3170239/fortran-integer4-vs-integer4-vs-integerkind-4 — Vladimir F Героям слава, Apr 06 '19 at 21:05
The step is being set twice because I need it to start at zero in order to evaluate the odd terms correctly. The size of the step is not the problem here, something else is not working. — Insight, Apr 06 '19 at 21:10
Yes, I *did* already told you that it is not the main problem here. But, beleive me, the step size `h` really should be set before the loop and let fixed. Please do not link to any pastebin. Put the important information here right into your own question. Questions here should be available for others even when external resources cease to work in the future. — Vladimir F Героям слава, Apr 06 '19 at 21:56

Vladimir F Героям слава · Accepted Answer · 2019-04-08T06:27:07.727

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.

Thanks for your help, but it doesn't seem like it's working though. I substituted the value of 2 for interval_length, 10 for n and for the functions calculated at the points I just wrote sqrt(1-(xa)**(2)) for example. But it returns something close to 15 for the value of s, and if I try to take smaller steps by changing n it returns even greater values. — Insight, Apr 08 '19 at 02:26
@Insight I forgot to divide by n, it should have been fairly obvious, I wrote the code without any testing,that is your task. — Vladimir F Героям слава, Apr 08 '19 at 08:22

Numerical Integration Using Fortran 90

1 Answers1