Fortran中5阶龙格-库塔方法的自适应步长方法

Question

我想用自适应步长法求解一组方程。我找到了Taner Akgun编写的有用代码。这是代码：

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。我将使用此代码解决以下一组方程。

1. dy/dx=-x
2. dy/dx=-1

在代码中，它

说 nvar 变量是方程的数量，在此代码中，它设置为 1。如果我想将其更改为 1 以外的其他内容，我应该如何更改代码？

另外，我想将所有 x 和 y 的值保存在输出文件中。我该怎么做？

Answer 1

首先尝试回答您的第一个问题。在不重复您原始问题中的大代码块的情况下，我怀疑您需要执行以下操作：

取代

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（正如您在评论中建议的那样），而只是为了解决这些方程，您可能需要查看 python 中scipy的 odeint 模块，它提供了所有这些功能。