函数数组和分段错误 - 无效的内存引用
Array of functions and Segmentation fault - invalid memory reference
我正在尝试将函数 f 设置为数组，但出现以下错误：
Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
Backtrace for this error:
#0 0x6f8b36e3
#1 0x6f8a2722
#2 0x402752
#3 0x747bd411
我必须对 M 的每个值求解开普勒方程：f=psi-e*sin(psi)-M。因此，如果我有一个维度为 8 的数组 M，我的程序将计算 8 个零.问题是，如果我写 f=psi-e*sin(psi)-M(1) 我会计算第一个零，如果我写 f=psi-e*sin(psi)-M(8) 我会计算最后一个 zero.But ，我的问题是如果我想计算所有零，我必须写 f=psi-e*sin(psi)-M(1:8) 并且我的程序应该键入所有零，但是，这没有发生，我得到我提到的错误 earlier.Here 是代码：
SUBROUTINE (I USED EXTERNALY):这个子程序是二分法(取零):
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.
在http://imgur.com/a/Xdsgf中你可以看到psi=0.95303344726562489。所以我刚刚计算了第一个零。但是，您也可以看到此消息 7 次 datos mal elegidos。就是说程序只能给我显示那个零（for M(1)），其他7个零都没有计算，因为我写了 f=psi-e*sin(psi)-M(1)。
我应该写什么才能得到 所有 零而不是 1 的结果，就像在这个例子中一样？
因为函数 f() 用于二分例程 bisecc()，我认为通过 DO 循环将每个输入传递给 bisecc() 会简单得多，而不是使f() 一个返回数组的函数（因为后者也需要修改 bisecc()）。我们可以通过各种方式将 M 的值传递给 f() （这几乎是 FAQ，我相信有很多 Q/A 页面）。一种简单的方法是在主程序中包含 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
答案似乎与 Kepler equation 的情节一致，其中 e = 0.2056（因此 bisecc() 可能没问题）。
以上代码还有很多需要改进的地方。特别是，将 f() 之类的函数包含到模块中（甚至将所有例程包含到模块中）通常更方便。我们还可以通过使 M 成为模块变量并从 f()（而不是使用 common 语句）或通过主机关联 use 来传递它，所以如果有兴趣请尝试一下.
MY SOLUTION:I 将为我的练习添加一个更通用的解决方案，避免提到 error.This 是一个更通用的解决方案，对于 M 的 N 个值，而不是 8:
include 'bisecc.f'
implicit real*8 (a-h,o-z)
external f
parameter (Mlong=100) !Number of elemnts of M(from 0 to 2pi)
real*8 f ,M
common M,e !to not copy them twice
kmax=100 !max number of iterations
tol=0.0001D0 !Tolerance of 0.01%
e=0.2056D0 !Mercury excentricity
pi=acos(-1.0D0)
c=sqrt((1.0D0-e)/(1.0D0+e))
open(10,file='153b.dat',status='unknown') !data will apear in a .dat file
write(*,*)' i M Theta(rad)'
write(10,*)' i M Theta(rad)'
do i=1,Mlong
xl=-1.D0 !LEFT STARTING POINT
xr=7.D0 !RIGHT POINT(psi wont be more than 2*pi)
M=2.D0*pi*i/Mlong !GENERIC M(0 TO 2PI 100STEPS)
call bisecc(f,xl,xr,kmax,tol,k,xm) !CALLING THE SUBROUTINE
write(10,*) i,M,theta ! I WILL PLOT THETA IN FUNCTION OF M
write(*,*) i,M,theta
end do
close(10)
write(*,*)
write(*,*) 'Program ENDED'
stop
end program
*我的外部函数
real*8 function f(psi)
implicit real*8 (a-h,o-z)
real*8 M
common M,e
f=psi-e*sin(psi)-M !KEPLER EQUATION
return
end function
我正在尝试将函数 f 设置为数组，但出现以下错误：
Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
Backtrace for this error:
#0 0x6f8b36e3
#1 0x6f8a2722
#2 0x402752
#3 0x747bd411
我必须对 M 的每个值求解开普勒方程：f=psi-e*sin(psi)-M。因此，如果我有一个维度为 8 的数组 M，我的程序将计算 8 个零.问题是，如果我写 f=psi-e*sin(psi)-M(1) 我会计算第一个零，如果我写 f=psi-e*sin(psi)-M(8) 我会计算最后一个 zero.But ，我的问题是如果我想计算所有零，我必须写 f=psi-e*sin(psi)-M(1:8) 并且我的程序应该键入所有零，但是，这没有发生，我得到我提到的错误 earlier.Here 是代码：
SUBROUTINE (I USED EXTERNALY):这个子程序是二分法(取零):
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.
在http://imgur.com/a/Xdsgf中你可以看到psi=0.95303344726562489。所以我刚刚计算了第一个零。但是，您也可以看到此消息 7 次 datos mal elegidos。就是说程序只能给我显示那个零（for M(1)），其他7个零都没有计算，因为我写了 f=psi-e*sin(psi)-M(1)。
我应该写什么才能得到 所有 零而不是 1 的结果，就像在这个例子中一样？
因为函数 f() 用于二分例程 bisecc()，我认为通过 DO 循环将每个输入传递给 bisecc() 会简单得多，而不是使f() 一个返回数组的函数（因为后者也需要修改 bisecc()）。我们可以通过各种方式将 M 的值传递给 f() （这几乎是 FAQ，我相信有很多 Q/A 页面）。一种简单的方法是在主程序中包含 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
答案似乎与 Kepler equation 的情节一致，其中 e = 0.2056（因此 bisecc() 可能没问题）。
以上代码还有很多需要改进的地方。特别是，将 f() 之类的函数包含到模块中（甚至将所有例程包含到模块中）通常更方便。我们还可以通过使 M 成为模块变量并从 f()（而不是使用 common 语句）或通过主机关联 use 来传递它，所以如果有兴趣请尝试一下.
MY SOLUTION:I 将为我的练习添加一个更通用的解决方案，避免提到 error.This 是一个更通用的解决方案，对于 M 的 N 个值，而不是 8:
include 'bisecc.f'
implicit real*8 (a-h,o-z)
external f
parameter (Mlong=100) !Number of elemnts of M(from 0 to 2pi)
real*8 f ,M
common M,e !to not copy them twice
kmax=100 !max number of iterations
tol=0.0001D0 !Tolerance of 0.01%
e=0.2056D0 !Mercury excentricity
pi=acos(-1.0D0)
c=sqrt((1.0D0-e)/(1.0D0+e))
open(10,file='153b.dat',status='unknown') !data will apear in a .dat file
write(*,*)' i M Theta(rad)'
write(10,*)' i M Theta(rad)'
do i=1,Mlong
xl=-1.D0 !LEFT STARTING POINT
xr=7.D0 !RIGHT POINT(psi wont be more than 2*pi)
M=2.D0*pi*i/Mlong !GENERIC M(0 TO 2PI 100STEPS)
call bisecc(f,xl,xr,kmax,tol,k,xm) !CALLING THE SUBROUTINE
write(10,*) i,M,theta ! I WILL PLOT THETA IN FUNCTION OF M
write(*,*) i,M,theta
end do
close(10)
write(*,*)
write(*,*) 'Program ENDED'
stop
end program
*我的外部函数
real*8 function f(psi)
implicit real*8 (a-h,o-z)
real*8 M
common M,e
f=psi-e*sin(psi)-M !KEPLER EQUATION
return
end function