Adaptive Stepsize Method for 5th Order Runge-Kutta Method in Fortran

Question

I want to solve a set of equations using 5th order Runge-Kutta method with adaptive stepsize method. I have found a useful code written by Taner Akgun. Here is the code:

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

Unfortunately, I'm not familiar with Fortran at all. I'm going to solve the following set of equations using this code.

dy/dx=-x
dy/dx=-1

Inside the code, it says the nvar variable is the number of equations and in this code, it is set to 1. If I want to change it to other than 1, how should I change the code?

Also, I want to save all the x's and y's values in the output file. How could I do it?

If you want to properly learn Fortran, the best way to do so would be to write the code from scratch yourself. The synthax is _relatively_ simple. There are several tutorials available online. If you are having problems later on, you could post your questions on this website again :) — solalito, Jun 07 '16 at 12:26
I don't need your opinion about learning Fortran! If you could help with the problem, please just help. — Ali Kefayati, Jun 07 '16 at 12:57

d_1999 · Answer 1 · 2016-06-07T13:25:30.440

Firstly to try to answer your first question. Without repeating the large code block from your original question I suspect you will need to do the following:

Replace

parameter(nvar=1)

with

parameter(nvar=2)

and replace the existing derivs routine with something like

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

You will also need to change the declaration of ystart, y and dydx in main to be double precision, dimension(2) :: ystart, y, dydx and then ensure these are set properly. This may be sufficient to give you the correct answer.

For your second question, one way to get the intermediate x and y values is rather than integrate from the start right to the end, instead integrate in steps. For example something like

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

However, if your goal is not to learn fortran (as you've suggested in a comment) but just to solve these equations you may want to look to the odeint module from scipy in python which provides all of this functionality.

Adaptive Stepsize Method for 5th Order Runge-Kutta Method in Fortran

1 Answers1