Fortran中五阶Runge-Kutta方法的自适应渐进方法

-2

我想用自适应步长方法使用五阶Runge-Kutta方法来求解一组方程。我找到了一个由Taner Akgun编写的有用代码。这里是代码：Fortran中五阶Runge-Kutta方法的自适应渐进方法

c 
c Adaptive Size Method for 5th Order Runge-Kutta Method 
c (Based on Numerical Recipes.) 
c 
c Taner Akgun 
c June, 2002 
c 
c Read on for various definitions. 
c (For more information consult the book.) 
c 
    program main 
    implicit none 
    integer nvar,nok,nbad 
    double precision x,y,dydx 
    double precision ystart,x1,x2,eps,h1,hmin 
c Number of derivatives to be integrated: 
c (In general, we can specify a set of equations.) 
    parameter(nvar=1) 
    external derivs,rkqs 
    open(1,file='test') 
c Integration boundaries and initial values: 
    x1 = 0d0 
    x2 = 2d0 
    ystart = 1d0 
c Stepsize definitions: 
c (h1 - initial guess; hmin - minimum stepsize) 
    h1 = 1d-1 
    hmin = 1d-9 
c write(1,*)'(Initial) Stepsize:',h1 
c Allowable error for the adaptive size method: 
     eps = 1d-9 
c Calculations: 
c Adaptive Size Method: 
    y = ystart 
    call odeint(y,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,rkqs) 
    stop 
    end 
c  
c Subroutine for the differential equation to be integrated. 
c Calculates derivative dydx at point x for a function y. 
c 
    subroutine derivs(x,y,dydx) 
     implicit none 
    double precision x,y,dydx 
    dydx = dexp(x) 
    return 
    end 
c   
c Subroutine for the Adaptive Size Method. 
c See Numerical Recipes for further information. 
c 
    subroutine odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs, 
    * rkqs) 
    implicit none 
    integer nbad,nok,nvar,KMAXX,MAXSTP,NMAX 
    double precision eps,h1,hmin,x1,x2,ystart(nvar),TINY 
    external derivs,rkqs 
    parameter(MAXSTP=10000,NMAX=50,KMAXX=200,TINY=1.e-30) 
    integer i,kmax,kount,nstp 
    double precision dxsav,h,hdid,hnext,x,xsav,dydx(NMAX) 
    double precision xp(KMAXX),y(NMAX),yp(NMAX,KMAXX),yscal(NMAX) 
    common /path/ kmax,kount,dxsav,xp,yp 
    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.d0*dxsav 
    do 16 nstp=1,MAXSTP 
     call derivs(x,y,dydx) 
     do 12 i=1,nvar 
      yscal(i)=dabs(y(i))+dabs(h*dydx(i))+TINY 
12  continue 
     if(kmax.gt.0)then 
      if(abs(x-xsav).gt.dabs(dxsav))then 
      if(kount.lt.kmax-1)then 
       kount=kount+1 
       xp(kount)=x 
       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.d0) 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.d0)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 
c Subroutine for the Adaptive Size Method. 
c See Numerical Recipes for further information. 
c Uses 'derivs' and 'rkck'. 
c 
    subroutine rkqs(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs) 
    implicit none 
    integer n,NMAX 
    double precision eps,hdid,hnext,htry,x,dydx(n),y(n),yscal(n) 
    external derivs 
    parameter(NMAX=50) 
    integer i 
    double precision errmax,h,htemp,xnew,yerr(NMAX),ytemp(NMAX) 
    double precision 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=0d0 
    do 11 i=1,n 
     errmax=max(errmax,dabs(yerr(i)/yscal(i))) 
11 continue 
    errmax=errmax/eps 
    if(errmax.gt.1d0)then 
     htemp=SAFETY*h*(errmax**PSHRNK) 
     h=sign(max(dabs(htemp),0.1d0*dabs(h)),h) 
     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=5d0*h 
     endif 
     hdid=h 
     x=x+h 
     do 12 i=1,n 
      y(i)=ytemp(i) 
12  continue 
     return 
    endif 
    end 
c 
c Subroutine for the Adaptive Size Method. 
c See Numerical Recipes for further information. 
c 
    subroutine rkck(y,dydx,n,x,h,yout,yerr,derivs) 
    implicit none 
    integer n,NMAX 
    double precision h,x,dydx(n),y(n),yerr(n),yout(n) 
    external derivs 
    parameter(NMAX=50) 
    integer i 
    double precision ak2(NMAX),ak3(NMAX),ak4(NMAX) 
    double precision ak5(NMAX),ak6(NMAX),ytemp(NMAX) 
    double precision A2,A3,A4,A5,A6 
    double precision B21,B31,B32,B41,B42,B43,B51,B52,B53, 
    * B54,B61,B62,B63,B64,B65 
    double precision 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

不幸的是，我对Fortran根本不熟悉。我将使用此代码解下面的一组方程。

DY/DX = -x
DY/DX = -1

里面的代码，它说NVAR变量方程的数量和在此代码，它被设置为1.如果我想将其更改为1以外，应如何更改代码？

此外，我想要保存输出文件中的所有x和y的值。我怎么能这样做？

来源

2016-06-07 ali kefayati

如果您想正确学习Fortran，最好的方法是从头开始编写代码。合成词比较简单。网上有几个教程。如果您以后遇到问题，您可以在此网站上再次发布您的问题:) – solalito

我不需要您关于学习Fortran的意见！如果你能解决问题，请帮忙。 –

首先尝试回答你的第一个问题。如果没有你原来的问题重复大码块我怀疑你需要做到以下几点：

更换

parameter(nvar=1)

与

parameter(nvar=2)

并且用类似

取代现有的 derivs程序

subroutine derivs(x,y,dydx) 
    implicit none 
    double precision x 
    double precision, dimension(:) y,dydx 
    dydx(1) = -x 
    dydx(2) = -1 
    return 
end

您还需要将main中的ystart,y和dydx的声明更改为double precision, dimension(2) :: ystart, y, dydx，然后确保这些设置正确。这可能足以给你正确的答案。

对于第二个问题，获取中间值x和y的一种方法不是从开始到结束进行集成，而是逐步集成。例如像

do i=1,nstep 
    call odeint(y,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,rkqs) 
    print*,"At x=",x2," y= ",y 
    !Update start and end points 
    x1=x2 
    x2=x1+stepSize 
enddo

但是，如果你的目标是不学FORTRAN（如你的评论所说），但正好解决了这些方程，你可能想看看在odeint模块从scipy在Python它提供了所有这些功能。

来源

2016-06-07 13:11:48

Fortran中五阶Runge-Kutta方法的自适应渐进方法

回答

相关问题