【发布时间】：2014-12-22 22:51:07
【问题描述】：
我是 Fortran90 编程的菜鸟。我将 NR 方法用于 Numerical Recipes 中的非线性方程组，并将使用 GFortran 编译时不会产生任何错误的代码放在一起。问题是，它也不会为输出生成任何值。 可能是因为我最初猜测的根值与实际根值相差甚远，还是我在这段代码中出错了？
我们将非常感谢您就此事提供任何帮助/建议。
program main
implicit real*8(a-h,o-z)
parameter(n=4)
logical check
real*8 x(n),fvec(n),fjac(n)
open(20,file="output1.txt",status="unknown")
do i=1,n
x(i)= 4.
enddo
call mnewt(ntrial,x,n,tolx,tolf)
call usrfun(x,n,fvec,fjac)
do i=1,n
write(20,*) 'x(',i,')=',x(i),fvec(i)
enddo
end
subroutine mnewt(ntrial,x,n,tolx,tolf)
integer n,ntrial,np
real*8 tolf,tolx,x(n)
parameter (np=15)
!uses lubksb, ludcmp, usrfun
! Given an initial guess x for a root in n dimensions, take ntrial Newton-Raphson steps to
! improve the root. Stop if the root converges in either summed absolute variable increments
! tolx or summed absolute function values tolf.
integer i,k,indx(np)
real*8 d,errf,errx,fjac(np,np),fvec(np),p(np)
do 14 k=1,ntrial
call usrfun(x,n,fvec,fjac)
errf=0.
do 11 i=1,n
errf=errf+abs(fvec(i))
11 continue
if(errf.le.tolf)return
do 12 i=1,n
p(i)=-fvec(i)
12 continue
call ludcmp(fjac,n,np,indx,d)
call lubksb(fjac,n,np,indx,p)
errx=0.
do 13 i=1,n
errx=errx+abs(p(i))
x(i)=x(i)+p(i)
13 continue
if(errx.le.tolx)return
14 continue
return
end
subroutine usrfun(x,n,fvec,fjac)
implicit none
integer n
real*8 x(n),fvec(n),fjac(n,n), hl, ul, br, bl
hl=1.00
ul=1.00
br=0.20
bl=0.00
! Initial guesses
x(1)=0.0
x(2)=1.5
x(3)=0.5
x(4)=0.5
fvec(1)=(x(2))+(2*sqrt((x(1))))-ul-(2*(sqrt(hl)))
fvec(2)=((x(3))*(x(4)))-((x(1))*(x(2)))
fvec(3)=((x(3))*(x(4))*(x(4)))+(0.5*(x(3))*(x(3)))-((x(1))*(x(2))*(x(2)))-(0.5*(x(1))*(x(1)))+(0.5*(br-bl)*x(1)+x(3))
fvec(4)=(x(4))-sqrt((x(3)))
fjac(1,1)=((x(1))**(-0.5))
fjac(1,2)=1
fjac(1,3)=0
fjac(1,4)=0
fjac(2,1)=(-x(2))
fjac(2,2)=(-x(1))
fjac(2,3)=x(4)
fjac(2,4)=x(3)
fjac(3,1)=((x(2))**2)-(x(1))+(0.5)*(br-bl)
fjac(3,2)=-2*((x(1))*(x(2)))
fjac(3,3)=((x(4))*(x(4)))+(x(3))+(0.5)*(br-bl)*(x(3))
fjac(3,4)=2*((x(3))*(x(4)))
fjac(4,1)=0
fjac(4,2)=0
fjac(4,3)=-0.5*((x(3))**(-0.5))
fjac(4,4)=1
end subroutine usrfun
subroutine ludcmp(a,n,np,indx,d) !fjac=a
integer n,np,indx(n),nmax
real*8 d,a(np,np),tiny
parameter (nmax=2500,tiny=1.0e-20)
integer i,imax,j,k
real*8 aamax,dum,sum,vv(nmax)
d=1.
do 12 i=1,n
aamax=0.
do 11 j=1,n
if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
! print*,a(21,j)
11 continue
! print*,i,aamax
! pause
if (aamax.eq.0.) pause 'singular matrix in ludcmp'
vv(i)=1./aamax
12 continue
do 19 j=1,n
do 14 i=1,j-1
sum=a(i,j)
do 13 k=1,i-1
sum=sum-a(i,k)*a(k,j)
13 continue
a(i,j)=sum
14 continue
aamax=0.
do 16 i=j,n
sum=a(i,j)
do 15 k=1,j-1
sum=sum-a(i,k)*a(k,j)
15 continue
a(i,j)=sum
dum=vv(i)*abs(sum)
if (dum.ge.aamax) then
imax=i
aamax=dum
endif
16 continue
if (j.ne.imax)then
do 17 k=1,n
dum=a(imax,k)
a(imax,k)=a(j,k)
a(j,k)=dum
17 continue
d=-d
vv(imax)=vv(j)
endif
indx(j)=imax
if(a(j,j).eq.0.)a(j,j)=tiny
if(j.ne.n)then
dum=1./a(j,j)
do 18 i=j+1,n
a(i,j)=a(i,j)*dum
18 continue
endif
19 continue
return
end
!lubksb
subroutine lubksb(a,n,np,indx,b)
integer n,np,indx(n)
real*8 a(np,np),b(n)
integer i,ii,j,ll
real*8 sum
ii=0
do 12 i=1,n
ll=indx(i)
sum=b(ll)
b(ll)=b(i)
if (ii.ne.0)then
do 11 j=ii,i-1
sum=sum-a(i,j)*b(j)
11 continue
else if (sum.ne.0.) then
ii=i
endif
b(i)=sum
12 continue
do 14 i=n,1,-1
sum=b(i)
do 13 j=i+1,n
sum=sum-a(i,j)*b(j)
13 continue
b(i)=sum/a(i,i)
14 continue
return
end
【问题讨论】：
-
如果您必须从 NR 复制代码，您还必须 (1) 在所有编译单元中插入
implicit none（以防止声明新实体的拼写错误）； (2) 将您的子例程放入module并使用关联它们（编译器将检查调用是否匹配声明）； (3) 设置每个过程参数的intent（主要是为了确保你不写你不应该写的参数）。就我个人而言，我什至不会考虑查看不包含这些必要的安全编程实践的代码，你也不应该这样做。 -
谢谢。但我怀疑你的任何观点是否是实际问题。因为我尝试在没有模块或意图语句的情况下编写全局收敛 NR 方法，并且效果很好。还是谢谢。
-
不，这不是实际问题，而是本世纪理智和安全编程的必要条件。最重要的是，它会告诉你，当你在调用过程时使用了错误的参数。