使用 Fortran 90 进行数值积分

Question

我正在尝试使用辛普森法则计算 sqrt(1-x^2) 在 -1 到 1 区间内的积分。但是，变量 "s" 表示的和，在我开发的代码根本不会收敛到 pi 超过 2。我是 Fortran 和一般编程的绝对新手，所以请多多包涵。我做错了什么？

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

这是它生成的输出的 link：https://pastebin.com/mW06Z6Lq

由于这个积分只是半径为 1 的圆面积的一半，它应该收敛到 pi 超过 2，但它大大超过了这个值。我考虑过将步长缩小以提高精度，但这不是问题，因为我尝试这样做时它甚至超过了预期值。

Answer 1

你必须检查当 x 接近 1 时会发生什么。你当然不能使用 DO WHILE (x<=1) 因为当 x==1 然后 x+2*h 大于 1 并且你正在计算负数的平方根。

我建议完全不要使用 DO WHILE。只需设置分割数，使用步长作为区间长度/分割数计算步长，然后使用

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

请注意，上面的循环嵌套是以非常通用的方式编写的，并非特定于辛普森规则。我只是假设 h 是常数，但即使这样也很容易改变。对于辛普森法则，可以很容易地对其进行优化。您当然希望每个时间间隔只进行两次函数评估。如果是学校作业，你需要把点数当作奇偶点而不是中心边界，你必须自己调整，这很容易。

使用 Fortran 90 进行数值积分

Numerical Integration Using Fortran 90

fortran

fortran90

numerical-integration

simpsons-rule