[問題] 解boundary value problem的程式

看板Fortran作者stsun1223 (神秘大師)時間14年前 (2010/07/18 23:18)推噓2(2推 0噓 11→)

留言13則, 5人參與討論串1/1

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