使用Fortran 90进行数值积分

Question

我正在尝试使用Simpson的规则来评估sqrt（1-x ^ 2）在-1到1之间的积分。但是，在我开发的代码中，变量“s”表示的总和不是总而言之，我已经完全融入了pi 2.我绝对是新手和编程的新手，所以请耐心等待。我究竟做错了什么？

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

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

由于这个积分只是半径为1的圆的面积的一半，所以它应该收敛到超过2的pi，但是它大大超过了这个值。我考虑过让步骤更小以获得精确度，但这不是问题，因为当我尝试这个时它超过了预期的值。

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

请注意，上面的循环嵌套是以非常通用的方式编写的，并非特定于Simpson规则。我只是假设h是恒定的，但即便如此也可以轻易地进行。对于辛普森的规则，人们可以轻松地优化它。您当然希望每个间隔只进行两次功能评估。如果是学校作业，你需要将这些点视为奇数 - 偶数而不是中心 - 边界，你必须自己调整，这很容易。