Méthode Adaptive Stepsize pour la méthode Runge-Kutta du 5ème ordre dans Fortran

Question

Je souhaite résoudre un ensemble d'équations en utilisant la méthode de Runge-Kutta du 5ème ordre avec la méthode d'échelonnement adaptatif. J'ai trouvé un code utile écrit par Taner Akgun. Voici le 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 

Malheureusement, je ne connais pas du tout Fortran. Je vais résoudre l'ensemble d'équations suivant en utilisant ce code.

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

A l'intérieur du code, il est dit la variable nvar est le nombre d'équations et dans ce code, il est réglé sur 1. Si je veux le changer à autre chose que 1, comment dois-je changer le code?

En outre, je veux enregistrer toutes les valeurs x et y dans le fichier de sortie. Comment pourrais-je le faire?

Answer 1

Tout d'abord pour essayer de répondre à votre première question. Sans répéter le grand bloc de code de votre question initiale, je soupçonne que vous devez faire ce qui suit:

Remplacer

parameter(nvar=1) 

avec

parameter(nvar=2) 

et remplacer la routine derivs existante avec quelque chose comme

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

Vous aurez également besoin pour modifier la déclaration de ystart, y et dydx dans main pour être double precision, dimension(2) :: ystart, y, dydx puis assurez-vous que ceux-ci sont définis correctement. Cela peut être suffisant pour vous donner la bonne réponse. Pour votre deuxième question, une façon d'obtenir les valeurs intermédiaires x et y est plutôt que de les intégrer du début à la fin, mais de les intégrer par étapes. Par exemple quelque chose comme

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 

Cependant, si votre objectif est de ne pas apprendre Fortran (comme vous l'avez dit dans un commentaire), mais juste pour résoudre ces équations que vous pouvez regarder le module odeint de scipy en python qui fournit toutes ces fonctionnalités.