如題，我寫的一個用到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