[問題] 解boundary value problem的程式
如題，我寫的一個用到shooting和牛頓法去解四個boundary value problem 的程式
可是一直跑不出正確答案
想請求板上的高手幫我看看到底是什麼問題，我已經想好久都想不出來
Thanks
以下程式
Program shoot100 !主程式
implicit none
integer i,n2,nvar,kmax,kount,KMAXX,NMAX
parameter (n2=1,NMAX=50,KMAXX=2000)
real f(2),v(n2),x1,x2,dxsav,xp(KMAXX),yp(NMAX,KMAXX)
logical check
common /caller/ x1,x2,nvar
common /path/ kmax,kount,dxsav,xp,yp
c ================================
nvar = 4
kmax = KMAXX
dxsav = 0.01
x1=1.e-5
x2=1.e8
v(1)=100
write(*,*) 'x1=',x1,'x2=',x2
write(*,*) 'guess initial vecter v = ', v
call newt(V,n2,check)
open(UNIT=20,FILE='sol.txt',STATUS='UNKNOWN')
do i = 1, kount
write (20,*) xp(i), yp(1,i),yp(2,i),yp(3,i),yp(4,i)
enddo
close(20)
write(*,*) 'n2 is', n2
write(*,*) 'The final value of V is :', V
write(*,*) 'check of newt is (F=OK!)', check
End
c=======================================================================
SUBROUTINE newt(x,n,check) !牛頓法副程式
INTEGER n,nn,NP,MAXITS
LOGICAL check
REAL x(n),fvec,TOLF,TOLMIN,TOLX,STPMX
PARAMETER(NP=40,MAXITS=200,TOLF=1.e-4,TOLMIN=1.e-6,TOLX=1.e-7
*,STPMX=100.)
COMMON /newtv/ fvec(NP),nn
SAVE /newtv/
CU USES fdjac,fmin,lnsrch,lubksb,ludcmp
INTEGER i,its,j,indx(NP)
REAL d,den,f,fold,stpmax,sum,temp,test,fjac(NP,NP),g(NP),p(NP)
*,xold(NP),fmin
EXTERNAL fmin
c ==============================
nn=n
f=fmin(x)
test=0.
do 11 i=1,n
if (abs(fvec(i)).gt.test) test=abs(fvec(i))
11 continue
if (test.lt.0.01*TOLF) return
sum=0.
do 12 i=1,n
sum=sum+x(i)**2
12 continue
stpmax=STPMX*max(sqrt(sum),float(n))
do 21 its=1,MAXITS
call fdjac(n,x,fvec,NP,fjac)
do 14 i=1,n
sum=0.
do 13 j=1,n
sum=sum+fjac(j,i)*fvec(j)
13 continue
g(i)=sum
14 continue
do 15 i=1,n
xold(i)=x(i)
15 continue
fold=f
do 16 i=1,n
p(i)=-fvec(i)
16 continue
call ludcmp(fjac,n,NP,indx,d)
call lubksb(fjac,n,NP,indx,p)
call lnsrch(n,xold,fold,g,p,x,f,stpmax,check,fmin)
test=0.
do 17 i=1,n
if (abs(fvec(i)).gt.test) test=abs(fvec(i))
17 continue
if (test.lt.TOLF) then
check=.false.
return
endif
if (check) then
test=0.
den=max(f,.5*n)
do 18 i=1,n
temp=abs(g(i))*max(abs(x(i)),1.)/den
if (temp.gt.test) test=temp
18 continue
if (test.lt.TOLMIN) then
check=.true.
else
check=.false.
endif
return
endif
test=0.
do 19 i=1,n
temp=(abs(x(i)-xold(i)))/max(abs(x(i)),1.)
if (temp.gt.test) test=temp
19 continue
if (test.lt.TOLX) return
21 continue
pause 'MAXITS exceeded in newt'
end
c======================================================================
FUNCTION fmin(x)
INTEGER n,NP
REAL fmin,x(*),fvec
PARAMETER (NP=40)
COMMON /newtv/ fvec(NP),n
SAVE /newtv/
CU USES funcv
INTEGER i
REAL sum
c ============================
call funcv(n,x,fvec)
sum=0.
do 11 i=1,n
sum=sum+fvec(i)**2
11 continue
fmin=0.5*sum
return
end
c=======================================================================
SUBROUTINE lnsrch(n,xold,fold,g,p,x,f,stpmax,check,func)
INTEGER n
LOGICAL check
REAL f,fold,stpmax,g(n),p(n),x(n),xold(n),func,ALF,TOLX
PARAMETER (ALF=1.e-4,TOLX=1.e-7)
EXTERNAL func
CU USES func
INTEGER i
REAL a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2,slope,sum,temp,
*test,tmplam
c ================================
check=.false.
sum=0.
do 11 i=1,n
sum=sum+p(i)*p(i)
11 continue
sum=sqrt(sum)
if (sum.gt.stpmax) then
do 12 i=1,n
p(i)=p(i)*stpmax/sum
12 continue
endif
slope=0.
do 13 i=1,n
slope=slope+g(i)*p(i)
13 continue
test=0.
do 14 i=1,n
temp=abs(p(i))/max(abs(xold(i)),1.)
if (temp.gt.test) test=temp
14 continue
alamin=TOLX/test
alam=1.
1 continue
do 15 i=1,n
x(i)=xold(i)+alam*p(i)
15 continue
f=func(x)
if (alam.lt.alamin) then
do 16 i=1,n
x(i)=xold(i)
16 continue
check=.true.
return
else if (f.le.fold+ALF*alam*slope) then
return
else
if (alam.eq.1.) then
tmplam=-slope/(2.*(f-fold-slope))
else
rhs1=f-fold-alam*slope
rhs2=f2-fold2-alam2*slope
a=(rhs1/alam**2-rhs2/alam2**2)/(alam-alam2)
b=(-alam2*rhs1/alam**2+alam*rhs2/alam2**2)/(alam-alam2)
if (a.eq.0.) then
tmplam=-slope/(2.*b)
else
disc=b*b-3.*a*slope
tmplam=(-b+sqrt(disc))/(3.*a)
endif
if (tmplam.gt..5*alam) tmplam=.5*alam
endif
endif
alam2=alam
f2=f
fold2=fold
alam=max(tmplam,.1*alam)
goto 1
END
c=======================================================================
SUBROUTINE fdjac(n,x,fvec,np,df)
INTEGER n,np,NMAX
REAL df(np,np),fvec(n),x(n),EPS
PARAMETER (NMAX=40,EPS=1.e-4)
CU USES funcv
INTEGER i,j
REAL h,temp,f(NMAX)
do 12 j=1,n
temp=x(j)
h=EPS*abs(temp)
if(h.eq.0.)h=EPS
x(j)=temp+h
h=x(j)-temp
call funcv(n,x,f)
x(j)=temp
do 11 i=1,n
df(i,j)=(f(i)-fvec(i))/h
11 continue
12 continue
return
END
c=======================================================================
SUBROUTINE ludcmp(a,n,np,indx,d)
INTEGER n,np,indx(n),NMAX
REAL d,a(np,np),TINY
PARAMETER (NMAX=500,TINY=1.0e-20)
INTEGER i,imax,j,k
REAL 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))
11 continue
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
c=======================================================================
SUBROUTINE lubksb(a,n,np,indx,b)
INTEGER n,np,indx(n)
REAL a(np,np),b(n)
INTEGER i,ii,j,ll
REAL 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
c=======================================================================
SUBROUTINE odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,
*rkqs)
INTEGER nbad,nok,nvar,KMAXX,MAXSTP,NMAX
REAL eps,h1,hmin,x1,x2,ystart(nvar),TINY
EXTERNAL derivs,rkqs
PARAMETER (MAXSTP=10000,NMAX=50,KMAXX=2000,TINY=1.e-30)
INTEGER i,kmax,kount,nstp
REAL dxsav,h,hdid,hnext,x,xsav,dydx(NMAX),xp(KMAXX),y(NMAX),
*yp(NMAX,KMAXX),yscal(NMAX)
COMMON /path/ kmax,kount,dxsav,xp,yp
c =====================================
x=x1
h=sign(h1,x2-x1)
nok=0
nbad=0
kount=0
do 11 i=1,nvar
y(i)=ystart(i)
11 continue
if (kmax.gt.0) xsav=x-2.*dxsav
do 16 nstp=1,MAXSTP
call derivs(x,y,dydx)
do 12 i=1,nvar
yscal(i)=abs(y(i))+abs(h*dydx(i))+TINY
12 continue
if(kmax.gt.0)then
if(abs(x-xsav).gt.abs(dxsav)) then
if(kount.lt.kmax-1)then
kount=kount+1
xp(kount)=x
write(*,*)'xp(kount)', xp(kount)
do 13 i=1,nvar
yp(i,kount)=y(i)
13 continue
xsav=x
endif
endif
endif
if((x+h-x2)*(x+h-x1).gt.0.) h=x2-x
call rkqs(y,dydx,nvar,x,h,eps,yscal,hdid,hnext,derivs)
if(hdid.eq.h)then
nok=nok+1
else
nbad=nbad+1
endif
if((x-x2)*(x2-x1).ge.0.)then
do 14 i=1,nvar
ystart(i)=y(i)
14 continue
if(kmax.ne.0)then
kount=kount+1
xp(kount)=x
do 15 i=1,nvar
yp(i,kount)=y(i)
15 continue
endif
return
endif
if(abs(hnext).lt.hmin) pause
*'stepsize smaller than minimum in odeint'
h=hnext
16 continue
pause 'too many steps in odeint'
return
END
c=======================================================================
SUBROUTINE rkqs(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs)
INTEGER n,NMAX
REAL eps,hdid,hnext,htry,x,dydx(n),y(n),yscal(n)
EXTERNAL derivs
PARAMETER (NMAX=50)
CU USES derivs,rkck
INTEGER i
REAL errmax,h,xnew,yerr(NMAX),ytemp(NMAX),SAFETY,PGROW,PSHRNK,
*ERRCON
PARAMETER (SAFETY=0.9,PGROW=-.2,PSHRNK=-.25,ERRCON=1.89e-4)
h=htry
1 call rkck(y,dydx,n,x,h,ytemp,yerr,derivs)
errmax=0.
do 11 i=1,n
errmax=max(errmax,abs(yerr(i)/yscal(i)))
11 continue
errmax=errmax/eps
if(errmax.gt.1.)then
h=SAFETY*h*(errmax**PSHRNK)
if(h.lt.0.1*h)then
h=.1*h
endif
xnew=x+h
if(xnew.eq.x)pause 'stepsize underflow in rkqs'
goto 1
else
if(errmax.gt.ERRCON)then
hnext=SAFETY*h*(errmax**PGROW)
else
hnext=5.*h
endif
hdid=h
x=x+h
do 12 i=1,n
y(i)=ytemp(i)
12 continue
return
endif
END
c=======================================================================
SUBROUTINE rkck(y,dydx,n,x,h,yout,yerr,derivs)
INTEGER n,NMAX
REAL h,x,dydx(n),y(n),yerr(n),yout(n)
EXTERNAL derivs
PARAMETER (NMAX=50)
CU USES derivs
INTEGER i
REAL ak2(NMAX),ak3(NMAX),ak4(NMAX),ak5(NMAX),ak6(NMAX),
*ytemp(NMAX),A2,A3,A4,A5,A6,B21,B31,B32,B41,B42,B43,B51,B52,B53,
*B54,B61,B62,B63,B64,B65,C1,C3,C4,C6,DC1,DC3,DC4,DC5,DC6
PARAMETER (A2=.2,A3=.3,A4=.6,A5=1.,A6=.875,B21=.2,B31=3./40.,
*B32=9./40.,B41=.3,B42=-.9,B43=1.2,B51=-11./54.,B52=2.5,
*B53=-70./27.,B54=35./27.,B61=1631./55296.,B62=175./512.,
*B63=575./13824.,B64=44275./110592.,B65=253./4096.,C1=37./378.,
*C3=250./621.,C4=125./594.,C6=512./1771.,DC1=C1-2825./27648.,
*DC3=C3-18575./48384.,DC4=C4-13525./55296.,DC5=-277./14336.,
*DC6=C6-.25)
do 11 i=1,n
ytemp(i)=y(i)+B21*h*dydx(i)
11 continue
call derivs(x+A2*h,ytemp,ak2)
do 12 i=1,n
ytemp(i)=y(i)+h*(B31*dydx(i)+B32*ak2(i))
12 continue
call derivs(x+A3*h,ytemp,ak3)
do 13 i=1,n
ytemp(i)=y(i)+h*(B41*dydx(i)+B42*ak2(i)+B43*ak3(i))
13 continue
call derivs(x+A4*h,ytemp,ak4)
do 14 i=1,n
ytemp(i)=y(i)+h*(B51*dydx(i)+B52*ak2(i)+B53*ak3(i)+B54*ak4(i))
14 continue
call derivs(x+A5*h,ytemp,ak5)
do 15 i=1,n
ytemp(i)=y(i)+h*(B61*dydx(i)+B62*ak2(i)+B63*ak3(i)+B64*ak4(i)+
*B65*ak5(i))
15 continue
call derivs(x+A6*h,ytemp,ak6)
do 16 i=1,n
yout(i)=y(i)+h*(C1*dydx(i)+C3*ak3(i)+C4*ak4(i)+C6*ak6(i))
16 continue
do 17 i=1,n
yerr(i)=h*(DC1*dydx(i)+DC3*ak3(i)+DC4*ak4(i)+DC5*ak5(i)+DC6*
*ak6(i))
17 continue
return
END
c=======================================================================
CU SUBROUTINE shoot(n2,v,f) is named "funcv" for use with "newt"
SUBROUTINE funcv(n2,v,f)
INTEGER n2,nvar,kmax,kount,KMAXX,NMAX
REAL f(2),v(n2),x1,x2,dxsav,xp,yp,EPS
PARAMETER (NMAX=50,KMAXX=2000,EPS=1.e-6)
COMMON /caller/ x1,x2,nvar
COMMON /path/ kmax,kount,dxsav,xp(KMAXX),yp(NMAX,KMAXX)
CU USES derivs,load,odeint,rkqs,score
INTEGER nbad,nok
REAL h1,hmin,y(NMAX)
EXTERNAL derivs,rkqs
kmax=KMAXX
h1=(x2-x1)/100.
hmin=0.
call load(x1,v,y)
call odeint(y,nvar,x1,x2,EPS,h1,hmin,nok,nbad,derivs,rkqs)
cc write (*,*) 'x1 and y1 and 10', xp(1), yp(1,1), xp(10), yp(10,1)
call score(x2,y,f,v)
return
END
c=======================================================================
subroutine load(x1,V,y)
real x1,v(1),y(4)
y(1) = v(1)
y(3) =(2.e5)*y(1)
y(2) = v(1)
y(4) =(2.e5)*y(2)
end
c=======================================================================
subroutine score(x2,y,v,f)
integer n2
parameter (n2=1)
real x2,y(4),f(2),v(1)
f(1) = y(3) -(3.e8-(y(1)/1.e8))
f(2) = y(4) -(-y(2)/1.e8)
write(*,*) 'f value in subroutine score is :', f(1),f(2)
end
c=======================================================================
subroutine derivs(x,y,dydx)
real x,y(4),dydx(4),a,b
dydx(1) = y(3)
dydx(2) = y(4)
dydx(3) =-(1.e-5/(1000.*exp(25.*(x/1.e8)**2)))*y(2)+
*(2./(x**2))*y(1)
dydx(4) = (1.e-5/(1000.*exp(25.*(x/1.e8)**2)))*y(1)+
*(2./(x**2))*y(2)
end
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 140.112.218.228
推
07/19 11:07, , 1F
07/19 11:07, 1F
→
07/19 11:08, , 2F
07/19 11:08, 2F
→
07/19 11:11, , 3F
07/19 11:11, 3F
→
07/19 12:36, , 4F
07/19 12:36, 4F
推
07/19 14:16, , 5F
07/19 14:16, 5F
→
07/19 15:32, , 6F
07/19 15:32, 6F
→
07/19 15:36, , 7F
07/19 15:36, 7F
→
07/21 23:57, , 8F
07/21 23:57, 8F
→
07/21 23:59, , 9F
07/21 23:59, 9F
→
07/22 12:21, , 10F
07/22 12:21, 10F
→
07/22 12:23, , 11F
07/22 12:23, 11F
→
07/22 12:23, , 12F
07/22 12:23, 12F
→
07/22 12:42, , 13F
07/22 12:42, 13F