时间戳

首页 > 解决方案 > 使用 Fortran 90 进行数值积分

fortran - 使用 Fortran 90 进行数值积分

问题描述

我正在尝试使用辛普森规则来评估从 -1 到 1 区间内 sqrt(1-x^2) 的积分。但是，在我开发的代码中，由变量“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

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

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

标签： fortranfortran90numerical-integrationsimpsons-rule

解决方案

您必须检查当您x接近 1 时会发生什么。您当然不能使用DO WHILE (x<=1)，因为当x==1thenx+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是恒定的，但即使这样也可以很容易地改变。对于辛普森规则，可以轻松优化它。您当然希望每个间隔只进行两次函数评估。如果是学校作业，你需要将这些点视为奇偶而不是中心边界，你必须自己调整，这很容易。

推荐阅读