我需要对我的 Simpson's Rule Fortran 代码进行哪些更改才能获得正确的结果？答案

【问题标题】：What do I need to change to my Simpson's Rule Fortran code to get the correct results?我需要对我的 Simpson's Rule Fortran 代码进行哪些更改才能获得正确的结果？
【发布时间】：2019-07-12 18:29:17
【问题描述】：

我正在尝试实现一段代码，通过使用Simpson's Rule 获取积分来获取曲线下的面积。

!(file:///D:/1-%20TUD/Semester%201/Numerical%20Methods%20BIWO-04/Lectures/simpson's%20rule.JPG)

我已经使用 MatchCAD 进行了尝试，并且得到了正确的结果函数：f(x)= x**5+(x-2)*sin(x)+(x-1)

program simpsons

  implicit none

  real a, b, h
  real integ, fa, fb
  integer i, m

  write(*,*) 'enter the lower boundary'
  read(*,*) a
  write(*,*) 'enter the upper boundary'
  read(*,*) b

  do while(a.ge.b)
    write(*,*) 'reenter the lower boundary'
    read(*,*) a
    write(*,*) 'reenter the upper boundary'
    read(*,*) b
  enddo

  write(*,*) 'enter the intervals number'
  read(*,*) m

  h=(b-a)/2.0

  fa=a**5.0+(a-2.0)*sin(a)+(a-1.0)
  fb=b**5.0+(b-2.0)*sin(b)+(b-1.0)

  integ=0

  do i=1, m/2

    integ=integ+4*((a+(2*i-1)*h)**5.0+((a+(2*i-1)*h)-2.0)*sin((a+(2*i-1)*h))+ ((a+(2*i-1)*h)-1.0))

    if (i.le.((m/2)-1)) then

      integ=integ+2*((a+2*i*h)**5.0+((a+2*i*h)-2.0)*sin((a+2*i*h))+((a+2*i*h)-1.0))

    endif
  enddo

  integ=(fa+fb+integ)*(h/3.0)

  write(*,*) 'integration = ', integ
end

当我输入 a=-1 和 b=1 和 m=20 时，我应该得到 Integration=-1.398

当我输入 a=-1 和 b=1 和 m=40 时，我应该得到 Integration=-1.398

但不知何故我得到了整合 =7015869.0

【问题讨论】：

欢迎您，请拨打tour。我建议你使用缩进来使你的代码结构清晰。有关一些基本知识，请参阅我的编辑。
最好为您集成的数学函数创建一个 Fortran 函数。代码将更具可读性并且更容易检查错误。
h=(b-a)/2.0，你确定这是正确的吗？
感谢大家的帮助。 @IanBush 应该是 h=(b-a)/m

标签： fortran simpsons-rule

【解决方案1】：

我前段时间使用过这段代码。正如@VladimirF 所指出的，最好将函数f(x) 与算法simpson(f,a,b,integral,n) 分开。

program main
  implicit none
  double precision a, b, integral
  integer n

  a = -1; b = 1; n = 20
  call simpson(f,a,b,integral,n)
  write (*,*) 'integration = ', integral

contains 

function f(x)
  implicit none
  double precision f, x
  f = x**5 + (x-2)*sin(x) + (x-1)
end

Subroutine simpson(f,a,b,integral,n)
  implicit none
  double precision f, a, b, integral, s
  double precision h, x
  integer n, i

  ! if n is odd we add +1 to make it even
  if((n/2)*2.ne.n) n = n+1

  ! loop over n (number of intervals)
  s = 0.0
  h = (b-a)/dble(n)
  do i = 2, n-2, 2
     x = a + i*h
     s = s + 2.0*f(x) + 4.0*f(x+h)
  end do
  integral = (s + f(a) + f(b) + 4.0*f(a+h))*h/3.0
end subroutine simpson

end program main

使用正确的输出： integration = -1.39766605364858

【讨论】：