如何在DO循环中区分结果和前一个结果？ [英] How can I make the difference between a result and the preceding one in do loops?

问题描述

此子例程用于确定复合梯形

我要提取最终结果(积分)和前一个结果(积分-1)之间的差异，并将该差异用作重复我的间隔数的限制。

Subroutine Trapezoid(a,b,n,integration)
real,external :: f
real :: h,a,b,summ,p
real,intent(out) :: integration
integer :: n
integer :: i,j
!Here as we have the whole equation is (h/2)*[f(a)+f(b)+2*sum(Xi)
!So we calculate the first part (h/2)*[f(a)+f(b) and then calculate the anoter part
do i=1,n
    n=2**i !Double the number of interval
    h=(b-a)/n  !Calculate the delta X
    p=(h/2.)*(f(a)+f(b))
        summ=0
        do j=1,n-1
            summ=summ+h*f(a+j*h)   !h/2 *2* sum[f(Xi)
        enddo  
    if(n == 256) then      !put a limit for the number of interval 
        Stop
    end if
    integration = p + summ   !Here the sum the both parts
    print*,n,'              ',integration 
enddo
end Subroutine

因此，我希望确定差异，而不是限制为250，当差异小于10*-8时，停止 我试了很多，但没有得到我想要的。

推荐答案

我会像下面这样做(非常快地黑在一起)。请注意，对于默认种类，1e-8是一个不切实际的精确度-因此公差较低。如果您想要更高的精度，您将需要使用更高精度的REAL类型。

还请注意，我已经将它变成了一个完整的程序。在提问中，请自己动手做。从纯粹自私的角度来看，你得到有用答案的可能性要大得多。

不管怎样，下面是代码

    Program integ
    
      Implicit None
    
      Real, Parameter :: pi = 3.1415927
    
      Real :: value, delta
    
      Integer :: n_used
    
      Intrinsic :: sin
      
      Call Trapezoid( sin, 0.0, pi / 2.0, 20, n_used, value, delta )
    
      Write( *, * ) 'final result', value, ' with ', 2 ** n_used, ' intervals'
      
    Contains
      
      Subroutine Trapezoid(f,a,b,n_max,n_used,integration,delta)
        Implicit None
        Real, Parameter :: tol = 1e-4
        Interface
           Function f( x ) Result( r )
             Real :: r
             Real, Intent( In )  :: x
           End Function f
        End Interface
        Real   , Intent( In    ) :: a
        Real   , Intent( In    ) :: b
        Integer, Intent( In    ) :: n_max
        Integer, Intent(   Out ) :: n_used
        Real   , Intent(   Out ) :: integration
        Real   , Intent(   Out ) :: delta
        Real :: h,summ,p
        Real :: integration_old
        Integer :: n
        Integer :: i,j
        !Here as we have the whole equation is (h/2)*[f(a)+f(b)+2*sum(Xi)
        !So we calculate the first part (h/2)*[f(a)+f(b) and then calculate the anoter part
        delta = - Huge( delta )
        integration_old = Huge( integration_old )
        Do i=1,n_max
           n=2**i !Double the number of interval
           h=(b-a)/n  !Calculate the delta X
           p=(h/2.)*(f(a)+f(b))
           summ=0
           Do j=1,n-1
              summ=summ+h*f(a+j*h)   !h/2 *2* sum[f(Xi)
           Enddo
           integration = p + summ   !Here the sum the both parts
           If( i /= 1 ) Then
              delta = integration - integration_old
              Write( *, * ) n,'              ',integration , delta
              If( Abs( delta ) < tol ) Exit
           End If
           integration_old = integration
        Enddo
        n_used = i
      End Subroutine Trapezoid
    
    End Program

    ian@eris:~/work/stack$ gfortran --version
    GNU Fortran (Ubuntu 7.4.0-1ubuntu1~18.04.1) 7.4.0
    Copyright (C) 2017 Free Software Foundation, Inc.
    This is free software; see the source for copying conditions.  There is NO
    warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
    
    ian@eris:~/work/stack$ gfortran -Wall -Wextra -fcheck=all -O -std=f2008 integ.f90 
    ian@eris:~/work/stack$ ./a.out
               4                 0.987115800       3.90563607E-02
               8                 0.996785223       9.66942310E-03
              16                 0.999196708       2.41148472E-03
              32                 0.999799252       6.02543354E-04
              64                 0.999949872       1.50620937E-04
             128                 0.999987483       3.76105309E-05
     final result  0.999987483      with          128  intervals

这篇关于如何在DO循环中区分结果和前一个结果？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！