用imsl求非线性方程组时运行出错，显示结果如下

sit2018

program hello

    use imsl

        implicit none

        integer n
        parameter (n = 3)

        integer k, nout
        real fnorm, x(n), xguess(n)
        external fcn

data xguess/4.0, 4.0, 4.0/

call umach (2,nout)
call neqnf(fcn, x, XGUESS=xguess, FNORM=fnorm)

write (nout, 99999) (x(k),k=1,n), fnorm
99999 format (' the solution to the system is', /, ' x = (', 3f5.1, &
')', /, ' with fnorm = ', f5.4, //)


        end program hello

        subroutine fcn(x,f,n)
        integer n
        real x(n), f(n)
        real exp,sin
        intrinsic exp, sin

        f(1) = x(1)+exp(x(1)-1.0)+(x(2)+x(3))*(x(2)+x(3)) -27.0
        f(2) = exp(x(2)-2.0)/x(1)+x(3)*x(3)-10.0
        f(3) = x(3)+sin(x(2)-2.0)+x(2)*x(2)-7.0

        return
        end subroutine fcn

Yowai

fcode 发表于 2018-5-18 17:28
第一个参数 f 应该是一个函数，而不是数组。
[mw_shl_code=fortran,true]program kakusann
  Include 'link ...

谢谢回复。我要求解的方程组传递数组 比较复杂，所以不想用子程序来传递。能够不写子程序，而直接在主程序里就算吗？

vvt

越是复杂越应该用子程序

Yowai

谢谢您的回复。。。。。

Yowai

!!!!!第一部分。
module XXX

contains
subroutine fcn(imax,mmax,P1,P2,X,F,N)
implicit none

real(8),parameter::exp=2.71828
integer::p,N
real(8)::X(N),F(N),P1(N),P2
integer::imax,mmax

if(imax+3==mmax-1) then
F(imax+3)=(0.5-X(imax+3))-P1(imax+3)/(1+exp**(P2*X(imax+3)))

else if (imax+3==mmax-2) then
F(imax+3)=(X(imax+4)-X(imax+3))-P1(imax+3)/(1+exp**(P2*X(imax+3)))
F(mmax-1)=(0.5-X(mmax-1)-X(mmax-2))-P1(mmax-1)/(1+exp**(P2*X(mmax-1)))

else if(imax+3<mmax-2) then
F(imax+3)=(X(imax+4)-X(imax+3))-P1(imax+3)/(1+exp**(P2*X(imax+3)))
do p=imax+4,mmax-2
F(p)=(X(p+1)-2*X(p)-X(p-1))-P1(imax+3)/(1+exp**(P2*X(imax+3)))
end do
F(mmax-1)=(0.5-X(mmax-1)-X(mmax-2))-P1(mmax-1)/(1+exp**(P2*X(mmax-1)))

end if
return
end subroutine


subroutine lsjac(imax,mmax,P1,P2,X,fjac,n)
implicit none

real(8),parameter::exp=2.71828
integer::p,i1,j1,N
real(8)::X(N),fjac(N,N),P1(N),P2
integer::imax,mmax
forall(i1=imax+3:mmax-1,j1=imax+3:mmax-1,i1==j1+1)   fjac(i1,j1)=1.0
forall(i1=imax+3:mmax-1,j1=imax+3:mmax-1,i1==j1-1)   fjac(i1,j1)=1.0
forall(i1=imax+3:mmax-1,j1=imax+3:mmax-1,i1==1 .and. j1==1)      &
fjac(i1,j1)=-1-P1(i1)*P2*exp**(P2*X(i1))/(1+exp**(P2*X(i1)))**2
forall(i1=imax+3:mmax-1,j1=imax+3:mmax-1,i1==j1 .and. j1/=1 .and. i1/=1)  &
fjac(i1,j1)=-1-P1(i1)*P2*exp**(P2*X(i1))/(1+exp**(P2*X(i1)))**2
forall(i1=imax+3:mmax-1,j1=imax+3:mmax-1,i1>=j1+1 .or.i1<=j1-1)  &
fjac(i1,j1)=0
return
end subroutine lsjac

end module XXX

Yowai

!!!!!!!!!第二部分
program kakusann


Include 'link_fnl_shared.h'
USE NEQNJ_INT
USE UMACH_INT

USE XXX

implicit none

real(8),parameter::exp=2.71828
real(8),parameter::FF=96485
real(8),parameter::tau=86400
real(8),parameter::RR=8.3145
real(8),parameter::T=298.15
real(8)::dx1,dx2,dx,r1,r2,Lall,Lf,La,D1,D2,dt1,dt2,dt,L,S1,sum,sum1,sum2,Sall,Sall1,Sall2,Sall3
real(8)::Kbio,y1,y2,fe0,P2,qmaxE,KsE,Ye,binaE,bresE,qmaxC,KsC,Yc,binaC,bresC,dLf,Xf


real(8)::S(0:10000)            !変化前の基質濃度
real(8)::Sn(0:10000)           !変化後の基質濃度

real(8)::Xe(0:10000)           !発電微生物密度(mg-VS/cm^3)
real(8)::Xei(0:10000)          !不活性発電微生物密度(mg-VS/cm^3)
real(8)::Xen(0:10000)          !内生呼吸発電微生物密度(mg-VS/cm^3)
real(8)::Xes(0:10000)          !不活性と活性発電微生物密度(mg-VS/cm^3)
real(8)::dXe(0:10000)          !単位時間の発電微生物変化量
real(8)::dXei(0:10000)         !単位時間の不活性発電微生物変化量
real(8)::dXen(0:10000)         !単位時間の内生呼吸発電微生物変化量


real(8)::Xc(0:10000)           !従属微生物密度(mg-VS/cm^3)
real(8)::Xci(0:10000)          !不活性従属微生物密度(mg-VS/cm^3)
real(8)::Xcn(0:10000)          !内生呼吸従属微生物密度(mg-VS/cm^3)
real(8)::Xcs(0:10000)          !不活性と活性従属微生物密度(mg-VS/cm^3)
real(8)::dXc(0:10000)          !単位時間の従属微生物変化量
real(8)::dXci(0:10000)         !単位時間の不活性従属微生物変化量
real(8)::dXcn(0:10000)         !単位時間の内生呼吸従属微生物変化量
real(8)::dXcs(0:10000)         !単位時間の従属微生物変化量

real(8)::dXs(0:10000)         !単位時間の微生物変化量（内生呼吸微生物以外）
real(8)::Xs(0:10000)          !微生物密度変化量(mg-VS/cm^3)（内生呼吸微生物以外）


real(8)::qE(0:10000)           !発電微生物基質消費量(mg-VS/cm^3)
real(8)::qC(0:10000)           !従属微生物基質消費量(mg-VS/cm^3)
real(8)::Ac=0.0003             !最初水槽内の基質濃度(mmol/cm^3)


integer::i,j,m,i1,j1,p
integer::imax=20                !拡散層の分割数
integer::jmax                   !生物膜の分割数
integer::mmax                   !拡散層と生物層の総分割数


!itaを求める時、使われる値
!external::fcn,lsjac
integer::N                           !行列の数
integer::K,nout
real(8)::fnorm
real(8),allocatable::X(:)             !ita
real(8),allocatable::Xguess(:)        !予測値
real(8),allocatable::FJAC(:,:)       !ヤコビアン式
real(8),allocatable::F(:)           !求める方程式
real(8),allocatable::P1(:)
N=jmax-3


!拡散層の初期値
L=0.01                         !拡散層の長さ、単位(cm)
D1=0.941/24/60/60              !拡散層の拡散係数、単位(cm2/s)
dx1=L/imax

!生物層の初期値
Lf=0.002                       !生物層の長さ、単位(cm)
D2=0.753/24/60/60              !生物層の拡散係数、単位(cm2/s)
dx2=dx1
jmax=Lf/dx2
dt1=dx1**2/D1*0.3
dt2=dx2**2/D2*0.3
N=jmax-3

if(dt2>=dt1) then
dt=dt1
else
dt=dt2
end if


allocate(X(N))
allocate(Xguess(N))
allocate(FJAC(N,N))
allocate(F(N))
allocate(P1(N))

write(*,*)'dt1=',dt1,'dt2=',dt2,'dt=',dt

!----------------------------------------------------------------------------------

r1=D1*dt/dx1**2
r2=D2*dt/dx2**2


! 初期値
Lall=7                         !総長さ、単位(cm)
S(0)=Ac

La=Lall-L-Lf                   !反応層の長さ
mmax=imax+jmax

qmaxE=0.132/24/60/60           !発電微生物最大基質消費速度(mmol-Ac/mgVS/day)→(mmol-Ac/mgVS/s)
qmaxC=0.234/24/60/60           !従属微生物最大基質消費速度(mmol-Ac/mgVS/day)→(mmol-Ac/mgVS/s)
KsE=0.00003                    !発電微生物の飽和定数(mmol/cm^3)
KsC=0.0011                     !従属微生物の飽和定数(mmol/cm^3)
Ye=4.52                        !発電微生物増殖収率(mg-VS/mmol-ED)
Yc=3.0                           !従属微生物増殖収率(mg-VS/mmol-ED)
binaE=0.05/24/60/60            !発電微生物不活性に変換される速度(/day)→(/s)
bresE=0.05/24/60/60            !発電微生物内部減衰定数(内生呼吸)速度(/day)→(/s)
binaC=0.02/24/60/60            !従属微生物不活性に変換される速度(/day)→(/s)
bresC=0.07/24/60/60            !従属微生物内部減衰定数(内生呼吸)速度(/day)→(/s)

Xf=80.0                         !飽和密度


Xe=20.0                          !発電微生物活性微生物密度(mg-VS/cm^3)
Xc=80.0                          !従属微生物活性微生物密度(mg-VS/cm^3)
do p=1,imax
Xe(p)=0
Xc(p)=0
end do


!電圧の部分


Kbio=0.00005                  !電気伝導性係数
y1=8                          !酢酸の電子当量
y2=0.177                      !微生物電子当量
fe0=0.9                       !微生物の基質から得られた電子を電極に渡す割合

do p=imax+1,mmax
Xei(p)=0
Xen(p)=0
Xe(p)=20
Xes(p)=20
Xci(p)=0

end do

do p=2,mmax
S(p)=0.0
end do

S(1)=(S(0)*La)/(La+dx1)
S(0)=S(1)


!リサイクルの開始
do m=1,10000
sum1=0
do p=2, mmax
sum1=sum1+S(p)
end do

!係数の表示
do p=imax+3,mmax-1
P1(p)=FF*Xe(p)/(tau*Kbio)*(y1*fe0*qmaxE*(S(p)/(S(p)+KsE))+y2*bresE)*dx2**2
end do
P2=-FF/(RR*T)
!do p=imax+3,mmax-1
!write(*,*)P1(p),P2
!end do



!予測値
do p=imax+3,mmax-1
Xguess(p)=0.3
end do

Yowai

!!!第三部分
!-------------------------------------------------------------------------------------------------------------
!! imslを使って、itaを求める
!call fcn(imax,mmax,P1,P2,X,F,N)
!call lsjac(imax,mmax,P1,P2,X,fjac,N)

call  UMACH(2,NOUT)
call DN2QNJ (FCN, LSJAC, ERRREL, N, ITMAX, XGUESS, X, FNORM, FVEC, FJAC, R, QTF, WK)
!WRITE (NOUT,99999) (X(K),K=1,N), FNORM
!99999 FORMAT ('  The solution to the system is', /, '  X = (', 3F5.1, &
!')', /, '  with FNORM =', F5.4, //)
!-------------------------------------------------------------------------------------------------------------



!拡散始まりから生物膜手前まで
do i=2,imax-1
Sn(i)=r1*S(i+1)+(1-2*r1)*S(i)+r1*S(i-1)
end do

!拡散層と生物膜の境界の拡散層側
Sn(imax)=(Sn(imax-1)-dx1/L*S(0))/(1-dx1/L)

!拡散層と生物膜の境界の生物膜側
Sn(imax+1)=Sn(imax)+D1/D2/L*(Sn(imax)-S(0))*dx2
qE(imax+1)=qmaxE*(S(imax+1)/(S(imax+1)+KsE))*Xe(imax+1)*dt
qC(imax+1)=qmaxC*(S(imax+1)/(S(imax+1)+KsC))*Xc(imax+1)*dt
Sn(imax+1)=Sn(imax+1)-qE(imax+1)-qC(imax+1)

dXei(imax+1)=binaE*Xe(imax+1)*dt
dXen(imax+1)=bresE*Xe(imax+1)*dt
dXe(imax+1)=Ye*qE(imax+1)-dXei(imax+1)-dXen(imax+1)
dXci(imax+1)=binaC*Xc(imax+1)*dt
dXcn(imax+1)=bresC*Xc(imax+1)*dt
dXc(imax+1)=Yc*qC(imax+1)-dXci(imax+1)-dXcn(imax+1)

Xei(imax+1)=Xei(imax+1)+dXei(imax+1)
Xen(imax+1)=Xen(imax+1)+dXen(imax+1)
Xe(imax+1)=Xe(imax+1)+dXe(imax+1)
Xci(imax+1)=Xci(imax+1)+dXci(imax+1)
Xcn(imax+1)=Xcn(imax+1)+dXcn(imax+1)
Xc(imax+1)=Xc(imax+1)+dXc(imax+1)

Xes(imax+1)=Xe(imax+1)+Xei(imax+1)
Xcs(imax+1)=Xc(imax+1)+Xci(imax+1)

Xs(imax+1)=Xes(imax+1)+Xcs(imax+1)



!生物膜上の拡散 j+2
do j=2,2

Sn(j+imax)=r2*S(j+imax+1)+(1.0-2.0*r2)*S(j+imax)+r2*S(j+imax-1)
qE(j+imax)=qmaxE*S(j+imax)/(S(j+imax)+KsE)*Xe(j+imax)*1/(1+exp**(P2*X(j+imax+1)))*dt
qC(j+imax)=qmaxC*S(j+imax)/(S(j+imax)+KsC)*Xc(j+imax)*dt
Sn(j+imax)=Sn(j+imax)-qE(j+imax)-qC(j+imax)

dXei(j+imax)=binaE*Xe(j+imax)*dt
dXen(j+imax)=bresE*Xe(j+imax)*1/(1+exp**(P2*X(j+imax+1)))*dt
dXe(j+imax)=Ye*qE(j+imax)-dXei(j+imax)-dXen(j+imax)
dXci(j+imax)=binaC*Xe(j+imax)*dt
dXcn(j+imax)=bresC*Xc(j+imax)*dt
dXc(j+imax)=Yc*qC(j+imax)-dXci(j+imax)-dXcn(j+imax)

Xei(j+imax)=Xei(j+imax)+dXei(j+imax)
Xen(j+imax)=Xen(j+imax)+dXen(j+imax)
Xe(j+imax)=Xe(j+imax)+dXe(j+imax)
Xci(j+imax)=Xci(j+imax)+dXci(j+imax)
Xcn(j+imax)=Xcn(j+imax)+dXcn(j+imax)
Xc(j+imax)=Xc(j+imax)+dXc(j+imax)

Xes(j+imax)=Xe(j+imax)+Xei(j+imax)
Xcs(j+imax)=Xc(j+imax)+Xci(j+imax)

Xs(j+imax)=Xes(j+imax)+Xcs(j+imax)
end do

!生物膜上の拡散,残り分
do j=3,jmax-1

Sn(j+imax)=r2*S(j+imax+1)+(1.0-2.0*r2)*S(j+imax)+r2*S(j+imax-1)
qE(j+imax)=qmaxE*S(j+imax)/(S(j+imax)+KsE)*Xe(j+imax)*1/(1+exp**(P2*X(j+imax)))*dt
qC(j+imax)=qmaxC*S(j+imax)/(S(j+imax)+KsC)*Xc(j+imax)*dt
Sn(j+imax)=Sn(j+imax)-qE(j+imax)-qC(j+imax)

dXei(j+imax)=binaE*Xe(j+imax)*dt
dXen(j+imax)=bresE*Xe(j+imax)*1/(1+exp**(P2*X(j+imax)))*dt
dXe(j+imax)=Ye*qE(j+imax)-dXei(j+imax)-dXen(j+imax)
dXci(j+imax)=binaC*Xe(j+imax)*dt
dXcn(j+imax)=bresC*Xc(j+imax)*dt
dXc(j+imax)=Yc*qC(j+imax)-dXci(j+imax)-dXcn(j+imax)

Xei(j+imax)=Xei(j+imax)+dXei(j+imax)
Xen(j+imax)=Xen(j+imax)+dXen(j+imax)
Xe(j+imax)=Xe(j+imax)+dXe(j+imax)
Xci(j+imax)=Xci(j+imax)+dXci(j+imax)
Xcn(j+imax)=Xcn(j+imax)+dXcn(j+imax)
Xc(j+imax)=Xc(j+imax)+dXc(j+imax)

Xes(j+imax)=Xe(j+imax)+Xei(j+imax)
Xcs(j+imax)=Xc(j+imax)+Xci(j+imax)

Xs(j+imax)=Xes(j+imax)+Xcs(j+imax)
end do


!生物膜と電極の境界
Sn(mmax)=S(mmax)*(1-2*r2)+S(mmax-1)*r2
qE(mmax)=qmaxE*S(mmax)/(S(mmax)+KsE)*Xe(mmax)*1/(1+exp**(P2*0.5))*dt
qC(mmax)=qmaxC*S(mmax)/(S(mmax)+KsC)*Xc(mmax)*dt
Sn(mmax)=Sn(mmax)-qE(mmax)-qC(mmax)

dXei(mmax)=binaE*Xe(mmax)*dt
dXen(mmax)=bresE*Xe(mmax)*1/(1+exp**(P2*0.5))*dt
dXe(j+imax)=Ye*qE(mmax)-dXei(mmax)-dXen(mmax)
dXci(mmax)=binaC*Xc(mmax)*dt
dXcn(mmax)=bresC*Xc(mmax)*dt
dXc(j+imax)=Yc*qC(mmax)-dXci(mmax)-dXcn(mmax)

Xei(mmax)=Xei(mmax)+dXei(mmax)
Xen(mmax)=Xen(mmax)+dXen(mmax)
Xe(mmax)=Xe(mmax)+dXe(mmax)
Xci(mmax)=Xci(mmax)+dXci(mmax)
Xcn(mmax)=Xcn(mmax)+dXcn(mmax)
Xc(mmax)=Xc(mmax)+dXc(mmax)

Xes(mmax)=Xe(mmax)+Xei(mmax)
Xcs(mmax)=Xc(mmax)+Xci(mmax)

Xs(mmax)=Xes(mmax)+Xcs(mmax)
!Sの更新
do p=2,mmax
S(p)=Sn(p)
end do

sum2=0
do p=2, mmax
sum2=sum2+S(p)
end do

S(1)=((S(0))*(La+dx1)-(sum2-sum1)*dx1)/(La+dx1)
S(0)=S(1)

end do

Yowai

!!第四部分
do p=0,mmax
!write(*,*) 'S(',p,')=',S(p),'X(',p,')=',X(p)
end do


!-------------------------------------------------------------------------------------------------------------
! 収支の計算
! 計算後の総量
Sall1=0
do p=1,24
Sall1=Sall1+S(p)*dx1
end do
Sall2=S(0)*La
Sall=Sall1+Sall2

! 投入量
Sall3=Ac*La

!write(*,*)'Sall=',Sall/La,'Sall3=',Sall3/La
!-------------------------------------------------------------------------------------------------------------



end program

hust菜菜

sit2018 发表于 2018-5-27 22:52
Error: This is an actual argument keyword name, and not a dummy argument name.   [XGUESS]
call neqnf ...

请问你解决了吗？？我也是这个问题

胡文刚

hust菜菜 发表于 2018-7-6 10:04
请问你解决了吗？？我也是这个问题

请单独开贴，附上代码（较长的话直接传附件），告知你用的操作系统，编译器版本，错误提示。