Адаптивный пошаговый метод для метода Рунге-Кутты 5-го порядка в Фортране
Я хочу решить систему уравнений, используя метод Рунге-Кутты 5-го порядка с адаптивным методом шага. Я нашел полезный код, написанный Танер Акгун. Вот код:
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
К сожалению, я не знаком с Фортраном вообще. Я собираюсь решить следующий набор уравнений, используя этот код.
- ду / дх =-x
- ду / дх =-1
Внутри кода написано, что переменная nvar - это число уравнений, и в этом коде она установлена на 1. Если я хочу изменить ее на значение, отличное от 1, как мне изменить код?
Кроме того, я хочу сохранить все значения x и y в выходном файле. Как я мог это сделать?
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
Вам также нужно будет изменить декларацию ystart, y а также dydx в main быть 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
Однако, если ваша цель не состоит в том, чтобы изучать фортран (как вы предложили в комментарии), а просто решить эти уравнения, вы можете обратиться к odeint модуль из scipy в Python, который предоставляет все эти функции.