函数数组和分段错误 - 无效的内存引用

我正在尝试设置我的功能f作为数组，但我收到以下错误：

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.

Backtrace for this error:

#0  0x6f8b36e3
#1  0x6f8a2722
#2  0x402752
#3  0x747bd411

我必须解开普勒方程：f=psi-e*sin(psi)-M对于每个值M.所以，如果我有一个数组M对于 8 维，我的程序将计算 8 个零。问题是，如果我写f=psi-e*sin(psi)-M(1)我将计算第一个零，如果我写f=psi-e*sin(psi)-M(8)我将计算最后一个零。但是，我的问题是，如果我想一次计算所有零，我必须写f=psi-e*sin(psi)-M(1:8)我的程序应该输入全零，但是，这并没有发生，我收到了前面提到的错误。这是代码：

子程序（我在外部使用）：这个子程序是二分法（得到零）：

   subroutine bisecc(f,xl,xr,kmax,tol,k,xm)
  implicit real*8 (a-h,o-z)
  real*8 f
  fl=f(xl)
  fr=f(xr)

  if(fl*fr .gt. 0.0D0) goto 100  
  do k=1,kmax
    xm=(xr+xl)/2.0D0              
    fm=f(xm)

    dif=abs((xr-xl)/xm)        
    if(dif .lt. tol) goto 200    
    write(*,*) k,xm!,dif  
    if (fm*fr .le. 0.0D0) then
      xl=xm
      fl=fm
    else
      xr=xm
      fr=fm
    end if
  end do

  return
200   write(*,*) 'WISHED PRECISION  REACHED'
      return
100   write(*,*) 'BAD CHOICE OF DATA'
      return
      end

主要程序：

       include 'bisecc.f' 
   implicit real*8 (a-h,o-z)            
   external f     
   real*8 f    
   ! I  WRITE  THE INTERVAL OF MY 8 ZEROS(left and right point)
   b=0.1D0


   xl1=-0.5D0                
   xr1=0.D0                 
   xl2=xr1+b
   xr2=1.D0
   xl3=xr2+b
   xr3=2.D0
   xl4=xr3+b
   xr4=3.D0
   xl5=xr4+b
   xr5=4.D0
   xl6=xr5+b
   xr6=5.D0
   xl7=xr6+b
   xr7=6.D0
   xl8=xr7+b
   xr8=7.D0        
   kmax=100                            
   tol=0.0001D0                               
   call bisecc(f,xl1,xr1,kmax,tol,k,xm1)         
   call bisecc(f,xl2,xr2,kmax,tol,k,xm2)          
   call bisecc(f,xl3,xr3,kmax,tol,k,xm3)         
   call bisecc(f,xl4,xr4,kmax,tol,k,xm4)         
   call bisecc(f,xl5,xr5,kmax,tol,k,xm5)         
   call bisecc(f,xl6,xr6,kmax,tol,k,xm6)         
   call bisecc(f,xl7,xr7,kmax,tol,k,xm7)          
   call bisecc(f,xl8,xr8,kmax,tol,k,xm8)                     
   write(*,*) 'Program ended'
   stop
   end program

   real*8 function f(psi)
   implicit real*8 (a-h,o-z)
   real*8 M(8)
   dimension f(8)
      e=0.2056D0                          
      pi=acos(-1.0D0)                        
      M=(/pi/4.D0,pi/2.D0,3.D0/4.D0*pi,pi,5.D0/4.D0*pi,3.D0*
 &    pi/2.D0,7.D0/4.D0*pi,2.D0*pi/)
      c=sqrt((1.0D0-e)/(1.0D0+e))                
   f=psi-e*sin(psi)-M(1:8)             !KEPLER EQUATION

   return
   end function   

例子： 这里我想计算M的第一个值的psi值，M(1)=pi/4.

In https://i.stack.imgur.com/qxEv2.jpg https://i.stack.imgur.com/qxEv2.jpg你可以看到psi=0.95303344726562489。所以我刚刚计算了第一个零。但是，您还可以看到此消息 7 次datos mal elegidos。这意味着该程序只能向我显示零（对于M(1)），其余7个零不计算，因为我写了 f=psi-e*sin(psi)-M(1)。 我应该写什么才能得到结果all零而不是 1，就像这个例子中那样？

因为函数f()用于二分例程中bisecc()，我认为将每个输入传递给会更简单bisecc()通过 DO 循环，而不是使f()返回数组的函数（因为后者需要修改bisecc()还）。我们可以传递的值M to f()以各种方式（这几乎是常见问题解答，我相信有很多问答页面）。一种简单的方法是包含f()在主程序中并使用主机关联M。所以简化的代码可能看起来像

program main
    implicit none
    integer  kmax, kiter, i
    real*8   xl( 8 ), xr( 8 ), xans( 8 ), tol, M( 8 ), b, pi

    pi   = acos(-1.0D0)                 
    kmax = 100
    tol  = 1.0d-8

    M = [ pi/4.D0,      pi/2.D0,      3.D0/4.D0*pi, pi, &
          5.D0/4.D0*pi, 3.D0*pi/2.D0, 7.D0/4.D0*pi, 2.D0*pi ]
    ! or M = [( i, i=1,8 )] * pi/4.0D0

    ! Use a fixed interval for simplicity.
    xl   =  0.0d0
    xr   = 10.0d0
    xans =  0.0d0

    do i = 1, 8
        call bisecc( f, xl( i ), xr( i ), kmax, tol, kiter, xans( i ) )
        ! print *, "check: f(xans(i)) = ", f( xans( i ) )
    enddo

contains

function f( psi ) result( res )
    implicit none
    real*8  psi, e, res
    e = 0.2056D0
    res = psi - e * sin( psi ) - M( i )   !<-- this "M(i)" refers to that defined above
end function 

end program

与外部bisecc例程（稍加修改以免使用 GOTO）

subroutine bisecc( f, xl, xr, kmax, tol, k, xm )
    implicit none
    real*8  f, xl, xr, tol, xm
    external f
    integer kmax, k
    real*8  fl, fr, fm, dif

    fl = f( xl )
    fr = f( xr )
    if( fl * fr > 0.0D0 ) then
        write(*,*) "bad input data (xl,xr)"
        return
    endif

    do k = 1, kmax
        xm = (xr + xl) / 2.0D0              
        fm = f( xm )

        dif = abs( (xr-xl) / xm )
        if ( dif < tol ) then
            write(*,*) "bisection converged: k=", k, "xm=", xm
            return
        endif

        if ( fm * fr <= 0.0D0 ) then
            xl = xm
            fl = fm
        else
            xr = xm
            fr = fm
        end if
    end do  !! iteration

    write(*,*) "bisection did not converge: k=", k, "xm=", xm
end

这使

 bisection converged: k=          31 xm=  0.95299366395920515     
 bisection converged: k=          31 xm=   1.7722388869151473     
 bisection converged: k=          30 xm=   2.4821592587977648     
 bisection converged: k=          30 xm=   3.1415926571935415     
 bisection converged: k=          29 xm=   3.8010260276496410     
 bisection converged: k=          29 xm=   4.5109464414417744     
 bisection converged: k=          29 xm=   5.3301916457712650     
 bisection converged: k=          29 xm=   6.2831853143870831

答案似乎很符合剧情开普勒方程 https://en.wikipedia.org/wiki/Kepler%27s_equation其中 e = 0.2056 （所以bisecc()可能没问题）。

上面的代码还有很多需要改进的地方。特别是，包含这样的函数通常会更方便f()到一个模块中（或者甚至将所有例程包含到一个模块中）。我们还可以通过M通过使其成为模块变量并且use它来自f()（而不是使用common声明）或通过主办协会，所以如果有兴趣请尝试一下。