「プログラミングと数値計算」カテゴリーアーカイブ

プログラミングと数値計算, 弾道計算

弾道計算(BB弾)のコード

2015年7月28日 sikino 17件のコメント

2016/03/13更新
2016/03/18更新
2016/03/20更新
2016/07/23更新
2016/07/25更新
弾道計算(BB弾)のfortran90での計算コードです。

2016/03/18 以降のコードに関して
コードは改変、複製、配布、営利、非営利目的に使用していただいて構いません。
ただし著作権を放棄している訳ではありません。
また、このコードを利用した場合、クレジットの表記をお願いします。

Linuxで、OSがUbuntuやLinuxMintの場合はfortranコンパイラは

sudo apt-get install gfortran

でインストールできます。
intel®fortranコンパイラとgfortranコンパイラで正常に動くことは確かめています。

もろもろのファイルをmain.f90に記述した後、
gfortranコンパイラなら

gfortran main.f90
./a.out

intel®fortranコンパイラなら

ifort main.f90
./a.out

で実行、ファイル”orbit.txt”が生成されます。

弾道計算に関するその他ページ
弾道計算(BB弾)の理論
 BB弾の回転量について(実験との比較)
弾道計算(BB弾)の結果
 弾道計算の結果2, 比較と詳細データ
 弾道計算(BB弾)のコード(fortran90)←今ここ
弾道計算のコード(Excel)
バレル内部でのBB弾の方程式
 水中下でのBB弾の弾道計算

できること

・BB弾の弾道計算(重力、慣性抵抗、粘性抵抗、ホップ有、回転の減衰を考慮)ができます。
式で表すと
BB弾の運動方程式(式1)
$
\displaystyle m\frac{d^2 \vec{r}}{dt^2}=
-mg\vec{k}+\left\{-6\pi \eta R |\vec{V}|-\frac{1}{2}C_d \rho \pi R^2 |\vec{V}|^2\right\}\frac{\vec{V}}{|\vec{V}|}
-\frac{4}{3}\pi R^3 2\rho |\vec{\omega}| |\vec{V}|\frac{\vec{V}\times\vec{\omega}}{|\vec{V}\times\vec{\omega}|}
$
とBB弾の回転角速度の減衰を記述する方程式(式2)
$
\begin{align}
\frac{d\omega_z}{dt}&=N_z/I \\
N_z&=\frac{\rho C_f}{2}R^3\int_0^{2\pi}d\varphi\int_0^{\pi}d\theta |u\sin\varphi-R\omega\sin\theta|(u\sin\varphi-R\omega\sin\theta)\sin^2\theta
\end{align}
$
の2つを考慮したものになります。

・ゼロインを決めた時の最適な軌道
・上下方向が最小の振れ幅になる軌道
も出力することができます。

※BB弾の回転の減衰を考えるか考えないかは、inputパラメータ”omgdecay”により制御できます。
回転の減衰に疑問がある場合はこれを利用してください。

注意事項

・とんでもないパラメータでは計算できない場合がありますので、そこら辺はご了承ください。
・僕のプログラミング知識ではこれが限界です。
僕はこのプログラムを使用したことによって生じた、いかなる場合でも一切の責任を負いません。

主なルーチンと説明

▼説明

モジュール

GBL
グローバル変数用のモジュールです。基本的には計算途中で何があっても変化しないパラメータを入れて至る所で呼び出して使っています。
rkmod
ルンゲクッタ法です。式1の時間発展を計算する際に用います。

サブルーチン

detwy
与えられたゼロイン位置(zeroinx)に対して最適な回転角速度(exwy)と最適な下向き方向の初速(exvz)を与えます。最適な角速度は関数wycondの根として与えられ、デフォルトではfalse position法を使って値がepsを下回るまで収束させています。
- detvz
  与えられた初期状態の回転角速度(omg(1:3))に対して最適な下向き方向の初速(exvz)を与えます。この時、エネルギーを基準として初期の速さを求めているため下向きの初速が増えるにしたがって水平方向の速度は減少しています。最適な下向き速度は関数vzcondの根として与えられ、デフォルトでは二分法を使って値がepsを下回るまで収束させています。二分法を使う理由は、ゼロイン距離が短い場合、関数vzcondが激しい振る舞いを見せるため、ほかの方法ではうまく収束しないためです。
- seeklm
  与えられた初期位置,初期速度、初期角速度、風速のもとで、z軸方向の極値とその時の値xを求めます。nodesは着弾するまでに射出高さを横切った回数、polesは着弾するまでの極値の数を意味します。考えられるnodesとpolesは最大で2つであるため、返す値は極値(xa,za)と極値(xb,zb)です。極値が無い場合、(xa,za)と(xb,zb)にはどちらも(0,0)が入ります。極値が1つの場合、(xa,za)にその値が入り、(xb,zb)には値(0,0)が入ります。
- seekzeroin
  与えられた初期位置,初期速度、初期角速度、風速のもとで、射出高さと同じ高さを通過した時の値xを求めます。nodesは着弾するまでに射出高さを横切った回数を意味します。考えられるnodesは最大で2つであるため、返す値はその時の射出方向の位置xaと位置xbです。横切ることが無い場合、xaとxbにはどちらも0が入ります。1回だけ横切る場合、xaにその値が入り、zbには値0が入ります。
rootfinding
与えられた関数funcのゼロ点を範囲ad~bd内で手法methodにより最大Nr個求めます。ad>bdでなければなりません。Nr>countの時、返すのはad~bd内のゼロ点の個数countとその値rd(1:count)です。Nr≤countの時、返すのはadから数えてNr個のゼロ点とその値rd(1:Nr)です。デフォルトではrootfindingのfuncに入る関数はvzcondのみです。
rootfindingw
与えられた関数funcのゼロ点を範囲ad~bd内で手法methodにより最大Nr個求めます。ad>bdでなければなりません。Nr>countの時、返すのはad~bd内のゼロ点の個数countとその値rd(1:count)です。Nr≤countの時、返すのはadから数えてNr個のゼロ点とその値rd(1:Nr)です。デフォルトではrootfindingwのfuncに入る関数はwycondのみです。

関数

vzcond
最適な下向き方向の初速を求めるための関数です。この関数はサブルーチンseeklmによって求めた極値(xa,za)と(xb,zb)の情報を使い、値(zb-z0)+(za-z0)を出力します。
wycond
与えられたゼロイン位置zeroinxに合わせるために最適な角速度を求めるための関数です。この関数はサブルーチンseekzeroinによって求めた極値xbの情報を使い、値xb-zeroinxを出力します。

計算理論

式1の時間発展はルンゲ-クッタ-フェールベルグ法による適応刻み幅制御を行っています。デフォルトでは8桁の精度で(※1)一致するようになっています。モジュールGBL内のtol=1d-8を変更することで精度を制御することができます。
式2の右辺の積分は、デフォルトでは15次のガウス=クロンロッド求積法を2回用いています。今回の問題に対してはあまり良い方法とは言えません。ですが、5-6桁はロンバーグ積分法による結果と一致するので、実用上問題ないと判断し、計算速度を優先しました。
式2の回転減衰の効果はルンゲクッタ法を行う際の連立方程式に組み込んでいません。これだけはオイラー法で行っています。なので、1ステップ目の刻み幅をどうするかで、本来はずれないはずですが、変更するとそれが伝搬してゼロイン位置がずれます。デフォルトでは全ての1ステップ目の刻み幅は統一してあるので上記問題は起こりません。本当に厳密にやるならば、コードを変更して減衰を連立方程式内に組み込んでください。
2016/03/13更新
コードを変更して減衰を連立方程式内に組み込みました。計算は若干遅くなりますが、上記問題はもう起こり得ません。

※1 相対誤差もしくは絶対誤差で、です。この理由は、相対誤差に統一してしまうと値が非常に小さいとき、機械誤差によっていつまでたっても収束しないためです。なので、元になる値が1より小さいときは絶対誤差を、1より大きいときは相対誤差を取るようにしています。

2016/03/13更新
2016/03/18更新
2016/03/20更新
2016/07/23更新
2016/07/25更新
2017/10/07更新

inputファイルとコード

inputというファイル名のファイルを作り、中に

&input
m=0.25d-3, ! Weight of bullet [kg]
a=3d-3, ! Radious of bullet [m]
energy=0.90d0, ! Energy of bullet [J]
g=9.80665d0, ! Gravity [m s^-2]
temperature=20d0, ! Temperature [degree of Celsius]
pressure=101325d0, ! Atmospheric pressure [Pa]
moisture=0.6d0, ! Value of moisture 0~1
omgdecay="yes", ! Introduce omega decay. "yes" or "no"
search_zeroin="yes", ! Search good vz and omgy. "yes" or "no"
zeroinx=50d0, ! +-- if yes, set zeroin length [m]
search_updown="yes",! Search good vz and omgy. "yes" or "no"
updownz=0.20d0, ! +-- if yes, set up-down henght [m]
theta=-1.7d0,
rx=0.d0,
ry=0.d0,
rz=1.0d0,
vy=0.d0,
ux=0.d0,
uy=0.d0,
uz=0.d0,
omgx=0.d0, !*2pi
omgy=-340d0, !*2pi
omgz=0.d0, !*2pi
outputfile="orbit", ! Output file name prefix
stept=0.01d0, ! Output interval [s]
ik=0, ! ik=0 exact, ik=1 approximation
&end

と記述してください。このinputファイルの例では、

重さ(m)0.25g
エネルギー(energy)0.90J
重力加速度(g)9.80665m/(s^2)
温度(temperature)20度
大気圧(pressure)101325Pa
湿度(moisture)60%
回転角速度の減衰(omgdecay)有り
ゼロインによる最適な弾道の探索(search_zeroin)しない
↑がyesなら、ゼロインの位置(zeroinx)50m
上下幅による最適な弾道の探索(search)する
↑がyesなら、上下幅(updownz)0.1m
銃口の位置(rz)1m
y方向の速度(vy)0m/s
射出角度(theta)-0.5度
風速(ux=uy=uz=0)0m/s
ホップ(omgy)-210回転/s
出力ファイル名の接頭語(outputfile)orbit
出力の時間間隔(stept)0.01秒ごとに出力
回転の減衰を表す積分の近似方法(ik=0 : 積分の直接計算、ik=1 : 積分を近似)

となっています。
回転の減衰を考えたくない時は、omgdecay=”no”としてください。

そして下のコードを展開し、コピー&ペーストして一つのファイル(例えばmain.f90という名前)に入れてください。
(カーソルをmodule GBL…のmに合わせた後、一番下までスクロールしてshift+クリックで楽に選択できます。)
個人的にファイルをダウンロードして…は嫌いなのでこういう形で公開しています。

2016/03/13 更新
2016/03/18 更新
2016/03/20 更新
2016/07/08 更新
2016/07/23 更新
2016/07/25 更新
2016/08/06 更新
ソースコードは1523行あります。

ここをクリックすると、ソースコードが展開されます。

module GBL
!sikinote http://slpr.sakura.ne.jp/qp/
!Author : sikino
!Date : 2016/03/18 (yyyy/mm/dd)
! : 2016/03/20
! : 2016/07/23
! : 2016/07/25
! : 2016/08/06
implicit none
!Constant of system
double precision::m,g,a,eta,rho,temperature,pressure,moisture
character(48)::omgdecay

!Example parameter
! 1, tol=1d-10, eps=1d-9
! => very exact, slow computing speed
! 2, tol=1d-6, eps=1d-6
! => enough computation practically.

!Accuracy of rungekutta(tol must be constant)
! 1d-2 < tol : Rouph condition, computing speed is fast.
! 1d-10 < tol < 1d-2 : Intermediate
! 1d-12 < tol < 1d-10 : Exact condition, computing speed is slow.
! tol < 1d-12 : Not reccomended, Maybe cannot computing due to machine epsilon.
double precision,parameter::tol=1d-10

!Conversence epsilon at rootfinding method
! 1d0 < eps : Rouph condition, computing speed is fast.
! 1d-10 < eps < 1d0 : Intermediate
! eps < 1d-10 : Exact condition, computing speed is slow.
! or cannot converge less than 1d-10.
double precision,parameter::eps=1d-7
end module GBL

!================================
function rho_humid(T,P,M)
! T : temperature [degree]
! P : pressure [Pa] ( 1 atm = 101325 Pa )
! M : moisture [no-dimension]
implicit none
double precision,intent(in)::T,P,M
double precision::rho_humid,et,rhoair

et=6.1078d0*10d0**(7.5d0*T/(T+237.3d0)) ! et[hPa]
et=100d0*et
et=M*et
!write(6,*)"et[Pa]",et

rhoair=0.0034856447d0*P/(T+273.15d0-0.670d0)
rho_humid=rhoair*(1d0-0.378d0*et/P)

return
end function rho_humid

function eta_air(T)
! T : temperature [degree]
implicit none
double precision::T,eta_air

eta_air=1.487d-6*((T+273.15d0)**1.5d0)/((T+273.15d0)+117d0)

return
end function eta_air

module RKmod
!sikinote http://slpr.sakura.ne.jp/qp/
!developer => sikino
!date => 2015/07/13
implicit none
!For Butcher table
integer,private::s
double precision,private,allocatable::a(:,:),b1(:),b2(:),c(:),Rc(:)
!-----------------------------
contains
subroutine rk_preparation(method)
!set Butcher table
implicit none
character(*)::method

if(trim(method).eq."rk4")then
s=4
allocate(a(1:s,1:s),b1(1:s),c(1:s))
a=0d0; b1=0d0; c=0d0

c(1:4)=(/0d0, 0.5d0, 0.5d0, 1d0/)
a(1,1:4)=(/0d0,0d0,0d0,0d0/)
a(2,1:4)=(/0.5d0,0d0,0d0,0d0/)
a(3,1:4)=(/0d0,0.5d0,0d0,0d0/)
a(4,1:4)=(/0d0,0d0,1d0,0d0/)
b1(1:4)=(/1d0/6d0, 1d0/3d0, 1d0/3d0, 1d0/6d0/)

elseif(trim(method).eq."rkf45")then
s=6
allocate(a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s))
a=0d0; b1=0d0; b2=0d0; c=0d0; Rc=0d0

c(1:6)=(/0d0, 0.25d0, 3d0/8d0, 12d0/13d0, 1d0, 0.5d0/)

a(1:6,1:6)=0d0
a(1,1:6)=(/0d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(2,1:6)=(/0.25d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(3,1:6)=(/0.09375d0, 0.28125d0, 0d0, 0d0, 0d0, 0d0/)
a(4,1:6)=(/1932d0/2197d0, -7200d0/2197d0, 7296d0/2197d0, 0d0, 0d0, 0d0/)
a(5,1:6)=(/439d0/216d0, -8d0, 3680d0/513d0, -845d0/4104d0, 0d0, 0d0/)
a(6,1:6)=(/-8d0/27d0, 2d0, -3544d0/2565d0, 1859d0/4104d0, -11d0/40d0, 0d0/)

b2(1:6)=(/16d0/135d0, 0d0, 6656d0/12825d0, 28561d0/56430d0, -9d0/50d0, 2d0/55d0/)
b1(1:6)=(/25d0/216d0, 0d0, 1408d0/2565d0, 2197d0/4104d0, -0.2d0, 0d0/)

Rc(1:6)=(/1d0/360d0,0d0,-128d0/4275d0,-2197d0/75240d0,1d0/50d0,2d0/55d0/)
else
write(6,*)"program stop at rk_preparation"
stop
end if

return
end subroutine rk_preparation
!------------------------------
subroutine rk_deallocation(method)
implicit none
character(*)::method

if(trim(method).eq."rk4")then
deallocate(a,b1,c)
elseif(trim(method).eq."rkf45")then
deallocate(a,b1,b2,c,Rc)
else
write(6,*)"program stop at rk_deallocation"
stop
end if

return
end subroutine rk_deallocation
!-------------------
subroutine rkf451(func,N,x,h,y,xbound,info,tol,ik)

!------------
!info = -2 (Failed. calclation range has already over.)
! = -1 (Failed.
! h becomes too small. change tol or check condition of func.)
! = 0 (Success. running now)
! = 1 (Success. x reach xbound normally)
!------------
implicit none
interface
function func(iN,ix,iy,is,iik)
implicit none
integer,intent(in)::iN,is,iik
double precision,intent(in)::ix,iy(1:iN)
double precision::func
end function func
end interface
integer,intent(in)::N,ik
double precision,intent(in)::xbound
double precision,intent(in)::tol
double precision,intent(inout)::x,h,y(1:N)
integer,intent(inout)::info
double precision,parameter::hmin=1d-14,hmax=1.d0

integer::i,j,FLAG
double precision::R,delta,tx,tmp(1:N),K(1:s,1:N),Sy,err

if(abs(h).ge.hmax)then
if(h.le.0d0)then
h=-hmax
else
h=hmax
endif
endif

FLAG=1
if(abs(x-xbound).le.1d-15)then
info=1
FLAG=0
else
if(abs(h).le.1d-15)then
write(6,*)"maybe overflow or underflow, please change tol."
write(6,'(A)')"====Err info===="
write(6,'(A,e17.9e3)')"x --> ",x
write(6,'(A,e17.9e3)')"h --> ",h
do i=1,N
write(6,'(A,i0,A,e17.9e3)')"y(",i,") --> ",y(i)
enddo
write(6,'(A)')"================"
info=-1
FLAG=0
stop
endif
if(h.le.0d0.and.xbound-x.ge.0d0)then
info=-2
FLAG=0
elseif(h.gt.0d0.and.xbound-x.le.0d0)then
info=-2
FLAG=0
endif
endif

do while(FLAG.eq.1)
tx=x
do j=1,s
tx=x+c(j)*h
tmp(1:N)=y(1:N)
do i=1,j-1
tmp(1:N)=tmp(1:N)+K(i,1:N)*a(j,i)
enddo
do i=1,N
K(j,i)=h*func(N,tx,tmp,i,ik)
enddo
enddo

!step 4
R=0d0
do i=1,N
R=R+(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i)+Rc(5)*K(5,i)+Rc(6)*K(6,i))**2d0
enddo
R=abs(dsqrt(R)/h)

Sy=0d0
do i=1,N
Sy=Sy+(y(i)*y(i))
enddo
Sy=dsqrt(Sy)
if(Sy.ge.1d0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err)then
x=x+h
do i=1,s
y(1:N)=y(1:N)+b1(i)*K(i,1:N)
enddo
FLAG=0
endif

!step 6
! Avoid zero deviding.
if(R.ge.1d-20)then
delta=(err/(2d0*R))**0.25d0
else
delta=4d0
endif

!step 7
if(delta.le.0.1d0)then
!function changes dramatically.
h=0.1d0*h
elseif(delta.ge.4d0)then
!function changes loosely.
h=4d0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
if(h.le.0d0)then
h=-hmax
else
h=hmax
endif
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if
enddo

return
end subroutine rkf451
subroutine rkf451_e(func,x,y,xbound,info,tol,ik)
!sikinote
! propagate from y(x) to y(xbound) without interval
!
! info = -1 : h < hmin. Maybe path the singular point.
! = 1 : x reach xbound.
!
implicit none
interface
function func(iN,ix,iy,is,iik)
implicit none
integer,intent(in)::iN,is,iik
double precision,intent(in)::ix,iy(1:iN)
double precision::func
end function func
end interface
integer::N,ik
double precision,intent(in)::xbound,tol
double precision,intent(inout)::x,y(:)
integer,intent(inout)::info
double precision,parameter::hmin=1d-14,hmax=1.d0

integer::i,j,FLAG,key,disc
double precision::R,delta,tx,Sy,err,h,h0
double precision,allocatable::tmp(:),K(:,:)

disc=0
key=0
h0=999d0
N=size(y,1)
allocate(tmp(1:N),K(1:s,1:N))
tmp=0d0; K=0d0

h=xbound-x
if(abs(h).ge.hmax)then
if(h.le.0d0)then
h=-hmax
else
h=hmax
endif
endif

FLAG=1
if(abs(x-xbound).le.hmin*0.1d0)then
info=1
FLAG=0
endif

do while(FLAG.eq.1)
tx=x
do j=1,s
tx=x+c(j)*h
tmp(1:N)=y(1:N)
do i=1,j-1
tmp(1:N)=tmp(1:N)+K(i,1:N)*a(j,i)
enddo
do i=1,N
K(j,i)=h*func(N,tx,tmp,i,ik)
enddo
enddo

!step 4
R=0d0
do i=1,N
R=R+(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i)+Rc(5)*K(5,i)+Rc(6)*K(6,i))**2d0
enddo
R=abs(dsqrt(R)/h)

Sy=0d0
do i=1,N
Sy=Sy+(y(i)*y(i))
enddo
Sy=dsqrt(Sy)
if(Sy.ge.1d0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
do i=1,s
y(1:N)=y(1:N)+b1(i)*K(i,1:N)
enddo
key=0
endif

!step 6
! Avoid zero deviding.
if(R.ge.1d-20)then
delta=(err/(2d0*R))**0.25d0
else
delta=4d0
endif

!step 7
if(delta.le.0.1d0)then
!function changes dramatically.
h=0.1d0*h
elseif(delta.ge.4d0)then
!function changes loosely.
h=4d0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1d0,h)*hmin
key=1
disc=1
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if
enddo

if(disc.eq.1)then
info=-1
endif

deallocate(tmp,K)
return
end subroutine rkf451_e
end module RKmod

!--------------------
subroutine seeklm(xa,xb,za,zb,nodes,poles,r0,v0,omg0,u0,ik)
!nodes --> number of nodes
!poles --> number of local maximam/minimam
!
!nodes=0
!
!
! x--------------------------------
! \_______
! \______
! \
! orbit of bullet
!
!nodes=1
! za_____________xa
! __/\_
! ________/ \_
! / \
! x------------------\------------
! \
! \
! \
! \orbit of bullet
!
!nodes=2
!
! zb________________ xb
! /\
! xa / \
! x-----------/--------\------------
! \ / \
! \ / \
! za______\/ \
! \orbit of bullet
!
use RKmod
use GBL
implicit none
integer,intent(in)::ik
double precision,intent(in)::omg0(1:3),r0(1:3),v0(1:3),u0(1:3)
double precision,intent(out)::xa,xb,za,zb
integer,intent(out)::nodes,poles

integer::i,j,info
integer,parameter::N=12
double precision::vz0,vz1,vz2,t0,t1,t2,ta,tb,tbef,rz0,rz1
double precision::x(1:N),tt,th,tx(1:N),bt,bh,bx(1:N)
double precision::t,h,xbound,rkfd
external::rkfd

x(1:3)=r0(1:3)
x(4:6)=v0(1:3)
x(7:9)=omg0(1:3)
x(10:12)=u0(1:3)

xbound=10d0
i=1; t=0d0; h=1d-1; info=0
info=0; nodes=0; poles=0

xa=0d0; xb=0d0; za=0d0; zb=0d0

call rk_preparation("rkf45")
do while(info.eq.0)
tt=t; th=h; tx=x

rz0=x(3); vz0=x(6); t0=t
call rkf451(rkfd,size(x),t,h,x,xbound,info,tol,ik)
rz1=x(3); vz1=x(6); t1=t

if(vz0*vz1.lt.0.d0)then
bt=tt; bh=th; bx=tx

t0=bt; t2=0.5d0*(2d0*bt+th); t1=bt+th
do j=1,60
tt=bt; th=bh; tx=bx

th=t2-bt
call rkf451(rkfd,size(tx),tt,th,tx,xbound,info,tol,ik)
vz2=tx(6)
if(vz2*vz0.ge.0.d0)then
t2=t2+0.5d0*(t1-t2)
else
tbef=t2
t2=t2-0.5d0*(t1-t2)
t1=tbef
endif

if(abs(t1-t2)/t1.lt.1d-10)then
if(poles.eq.0)then
ta=t2; xa=tx(1); za=tx(3)
else
tb=t2; xb=tx(1); zb=tx(3)
endif
exit
endif
enddo
poles=poles+1
endif
if((rz0-r0(3))*(rz1-r0(3)).lt.0d0)then
nodes=nodes+1
endif

if(x(3).lt.0d0)then
info=1
endif
i=i+1
enddo
call rk_deallocation("rkf45")

return
end subroutine seeklm

subroutine seekzeroin(xa,xb,nodes,r0,v0,omg0,u0,ik)
!nodes --> number of nodes
!
!nodes=0
!
!
! x--------------------------------
! \_______
! \______
! \
! orbit of bullet
!
!nodes=1
!
! __/\_
! ________/ \_
! / \ xa
! x------------------\------------
! \
! \
! \
! \orbit of bullet
!
!nodes=2
!
!
! /\
! xa / \ xb
! x-----------/--------\------------
! \ / \
! \ / \
! \/ \
! \orbit of bullet
!
use RKmod
use GBL
implicit none
integer,intent(in)::ik
double precision,intent(in)::omg0(1:3),r0(1:3),v0(1:3),u0(1:3)
double precision,intent(out)::xa,xb
integer,intent(out)::nodes

integer::i,j,info
integer,parameter::N=12
double precision::rz0,rz1,rz2,t0,t1,t2,ta,tb,tbef,rzini
double precision::x(1:N),tt,th,tx(1:N),bt,bh,bx(1:N)
double precision::t,h,xbound,rkfd
external::rkfd

rzini=r0(3)
x(1:3)=r0(1:3)
x(4:6)=v0(1:3)
x(7:9)=omg0(1:3)
x(10:12)=u0(1:3)

xbound=10d0
i=1; t=0d0; h=1d-1; info=0
info=0; nodes=0

xa=0d0; xb=0d0

call rk_preparation("rkf45")

!-----1step-----
call rkf451(rkfd,size(x),t,h,x,xbound,info,tol,ik)
!----------------
do while(info.eq.0)
tt=t; th=h; tx=x

rz0=x(3); t0=t
call rkf451(rkfd,size(x),t,h,x,xbound,info,tol,ik)
rz1=x(3); t1=t

if(rz0-rzini.lt.0.d0.neqv.rz1-rzini.lt.0.d0)then
bt=tt; bh=th; bx=tx

t0=bt; t2=0.5d0*(2d0*bt+th); t1=bt+th
do j=1,60
tt=bt; th=bh; tx=bx

th=t2-bt
call rkf451(rkfd,size(tx),tt,th,tx,xbound,info,tol,ik)
rz2=tx(3)
if(rz2-rzini.lt.0.d0.eqv.rz0-rzini.lt.0.d0)then
t2=t2+0.5d0*(t1-t2)
else
tbef=t2
t2=t2-0.5d0*(t1-t2)
t1=tbef
endif

if(abs(t1-t2)/t1.lt.1d-10)then
if(nodes.eq.0)then
ta=t2; xa=tx(1)
else
tb=t2; xb=tx(1)
endif
nodes=nodes+1
exit
endif
enddo
endif

if(x(3).lt.0d0)then
info=1
endif
i=i+1
enddo
call rk_deallocation("rkf45")

return
end subroutine seekzeroin

subroutine detwy(energy,r,vy,omgx,omgz,u,zeroinx,exwy,exvz,ik)
implicit none
integer,intent(in)::ik
double precision,intent(in)::energy,zeroinx
double precision,intent(in)::r(1:3),omgx,omgz,u(1:3),vy
double precision,intent(out)::exwy,exvz

integer::Nr,count,Nsplit
double precision,parameter::pi=dacos(-1d0)
double precision,allocatable::rd(:)
double precision::ad,bd,wycond
external::wycond

ad=-800d0*2d0*pi; bd=2d0*pi; Nsplit=25; Nr=1
allocate(rd(1:Nr))
rd(1:Nr)=0d0
count=0
call rootfindingw(wycond,energy,r,vy,omgx,omgz,u,zeroinx,ad,bd,Nsplit,Nr,rd,exvz,count,"bisection",1)
exwy=rd(1)

if(ik.eq.0)then
ad=exwy-5d0*2d0*pi; bd=exwy+5d0*2d0*pi; Nsplit=1; Nr=1
rd(1:Nr)=0d0
count=0
call rootfindingw(wycond,energy,r,vy,omgx,omgz,u,zeroinx,ad,bd,Nsplit,Nr,rd,exvz,count,"false_position",0)
exwy=rd(1)
endif

deallocate(rd)
return
end subroutine detwy
!--------------------------
function wycond(energy,r,vy,omgx,omgy,omgz,u,zeroinx,exvz,ik)
use GBL, only:m
implicit none
integer,intent(in)::ik
double precision,intent(in)::energy,r(1:3),u(1:3),vy,omgx,omgy,omgz,zeroinx
double precision,intent(out)::exvz

integer::nodes
double precision::zeroinx1,zeroinx2,wycond,v(1:3),omg(1:3)

omg(1)=omgx; omg(2)=omgy; omg(3)=omgz

call detvz(energy,r,vy,omg,u,exvz,ik)
v(1)=dsqrt((2d0*energy/m)-exvz**2-vy**2)
v(2)=vy
v(3)=exvz
call seekzeroin(zeroinx1,zeroinx2,nodes,r,v,omg,u,ik)
if(zeroinx2.ge.1d-10)then
wycond=zeroinx2-zeroinx
else
wycond=2d0*zeroinx
endif

return
end function wycond
!------------------------------
subroutine detwy2(energy,r,vy,omgx,omgz,u,udh,exwy,exvz,ik)
implicit none
integer,intent(in)::ik
double precision,intent(in)::energy,udh,r(1:3),omgx,omgz,u(1:3),vy
double precision,intent(out)::exwy,exvz

integer::Nr,count,Nsplit
double precision,parameter::pi=dacos(-1d0)
double precision,allocatable::rd(:)
double precision::ad,bd,wycond2
external::wycond2

ad=-1000d0*2d0*pi; bd=2d0*pi; Nsplit=25; Nr=1
allocate(rd(1:Nr))
rd(1:Nr)=0d0
count=0
call rootfindingw(wycond2,energy,r,vy,omgx,omgz,u,udh,ad,bd,Nsplit,Nr,rd,exvz,count,"bisection",1)
exwy=rd(1)

if(ik.eq.0)then
ad=exwy-5d0*2d0*pi; bd=exwy+5d0*2d0*pi; Nsplit=1; Nr=1
count=0; rd(1:Nr)=0d0
call rootfindingw(wycond2,energy,r,vy,omgx,omgz,u,udh,ad,bd,Nsplit,Nr,rd,exvz,count,"false_position",0)
exwy=rd(1)
endif

deallocate(rd)
return
end subroutine detwy2
!--------------------------
function wycond2(energy,r,vy,omgx,omgy,omgz,u,udh,exvz,ik)
use GBL, only:m
implicit none
integer,intent(in)::ik
double precision,intent(in)::energy,r(1:3),u(1:3),vy,omgx,omgy,omgz,udh
double precision,intent(out)::exvz

integer::nodes,poles
double precision::xa,xb,za,zb,z0,wycond2,v(1:3),omg(1:3)

omg(1)=omgx; omg(2)=omgy; omg(3)=omgz

call detvz(energy,r,vy,omg,u,exvz,ik)
v(1)=dsqrt((2d0*energy/m)-exvz**2-vy**2)
v(2)=vy
v(3)=exvz

z0=r(3)
call seeklm(xa,xb,za,zb,nodes,poles,r,v,omg,u,ik)
if(poles.eq.0.or.za.lt.0d0)then
wycond2=10d0
elseif(poles.eq.1)then
wycond2=za-z0
elseif(poles.eq.2)then
wycond2=zb-za
else
write(6,*)"unknown,poles-->",poles
stop
endif
wycond2=wycond2-udh

return
end function wycond2
!-------------------------
subroutine detvz(energy,r,vy,omg,u,exvz,ik)
implicit none
integer,intent(in)::ik
double precision,intent(in)::energy
double precision,intent(in)::r(1:3),omg(1:3),u(1:3),vy
integer::Nr,count,Nsplit
double precision,intent(out)::exvz

double precision,allocatable::rd(:)
double precision::ad,bd,vzcond
external::vzcond

ad=-30d0; bd=1d0; Nsplit=80; Nr=1
allocate(rd(1:Nr))
rd(1:Nr)=0d0
count=0
call vzroot(vzcond,energy,r,vy,omg,u,ad,bd,Nsplit,Nr,rd,count,ik)
exvz=rd(1)

return
end subroutine detvz
!------------------------
function vzcond(energy,r,vy,vz,omg,u,ik)
use GBL, only:m
implicit none
integer,intent(in)::ik
double precision,intent(in)::energy,r(1:3),vy,vz,omg(1:3),u(1:3)
double precision::vzcond,v(1:3)

integer::nodes,poles
double precision::xa,xb,za,zb,z0

v(1)=dsqrt((2d0*energy/m)-vz**2-vy**2)
v(2)=vy
v(3)=vz

z0=r(3)
call seeklm(xa,xb,za,zb,nodes,poles,r,v,omg,u,ik)
if(poles.eq.0)then
vzcond=-z0-1d0
elseif(poles.eq.1)then
vzcond=za-z0
elseif(poles.eq.2)then
vzcond=za+zb-2d0*z0
else
write(6,*)"unknown,poles-->",poles
stop
endif

return
end function vzcond
!--------------------
subroutine vzroot(func,energy,r,vy,omg,u,ad,bd,Nsplit,Nr,rd,count,ik)
!Developer : sikino
use GBL
implicit none
interface
function func(ienergy,ir,ivy,ivz,iomg,iu,iik)
implicit none
integer,intent(in)::iik
double precision,intent(in)::ienergy,ir(1:3),ivy,ivz,iomg(1:3),iu(1:3)
double precision::func
end function func
end interface
integer,intent(in)::Nr,Nsplit,ik
double precision,intent(in)::ad,bd,energy,r(1:3),vy,omg(1:3),u(1:3)
integer,intent(out)::count
double precision,intent(out)::rd(1:Nr)

integer::N,i,j
double precision::h,hd,x0,x1,x2,y0,y1,y2,Dy0,tx1,ty1

!ad < bd
if(ad.ge.bd)then
write(6,'(A,2e15.6e3)')"must be ad < bd, your ad,bd --> ",ad,bd
write(6,'(A)')"program stop at rootfinding"
stop
end if

N=Nsplit
hd=abs(ad-bd)/dble(N)
count=0

x0=ad
y0=func(energy,r,vy,x0,omg,u,ik)
do i=1,N
x1=x0+hd
y1=func(energy,r,vy,x1,omg,u,ik)
if(y0.lt.0.d0.neqv.y1.lt.0.d0)then
tx1=x1
ty1=y1

!bisection rule
do j=1,60
x2=0.5d0*(x0+x1)
y2=func(energy,r,vy,x2,omg,u,ik)
if(y2.lt.0.d0.eqv.y0.lt.0.d0)then
x0=x2
else
x1=x2
endif
if(abs((x1-x0)/x2).lt.eps)then
if(y2.le.1d0)then
count=count+1
rd(count)=x2
exit
endif
endif
enddo
if(j.eq.61)write(6,*)"+--cannot convergence at root-finding method bs--+"

if(count.ge.Nr)exit
x1=tx1
y1=ty1
endif
x0=x1
y0=y1
enddo

return
end subroutine vzroot

subroutine rootfindingw(func,energy,r,vy, &
omgx,omgz,u,zeroinx,ad,bd,Nsplit,Nr,rd,exvz,count,method,ik)
!Date : 2015/07/28
!Developer : sikino
use GBL
implicit none
interface
function func(ienergy,ir,ivy,iomgx,iomgy,iomgz,iu,izin,iexvz,iik)
implicit none
integer,intent(in)::iik
double precision,intent(in)::ienergy,ir(1:3),ivy,iomgx,iomgy,iomgz,izin,iu(1:3)
double precision,intent(out)::iexvz
double precision::func
end function func
end interface
integer,intent(in)::Nr,Nsplit,ik
double precision,intent(in)::ad,bd,energy,r(1:3),vy,u(1:3)
double precision,intent(in)::omgx,omgz,zeroinx
integer,intent(out)::count
double precision,intent(out)::rd(1:Nr),exvz
double precision,parameter::pi=dacos(-1d0)
character(*),intent(in)::method

integer::N,i,j

double precision::h,hd,x0,x1,x2,y0,y1,y2,Dy0,tx1,ty1

!ad < bd
if(ad.ge.bd)then
write(6,'(A,2e15.6e3)')"must be ad < bd, your ad,bd --> ",ad,bd
write(6,'(A)')"program stop at rootfinding"
stop
end if

!Announsment
write(6,'(A,A)')" ---- ",trim(method)

N=Nsplit
hd=abs(ad-bd)/dble(N)
count=0

x0=ad
y0=func(energy,r,vy,omgx,x0,omgz,u,zeroinx,exvz,ik)
do i=1,N
x1=x0+hd
y1=func(energy,r,vy,omgx,x1,omgz,u,zeroinx,exvz,ik)
write(6,'(A,i0,A,i0,A,f14.8,A,f14.8)')" Sequence ",i," of ",N," : omgy/2pi ",x1/2d0/pi, " : zeroin ",y1
if(y0.lt.0.d0.neqv.y1.lt.0.d0)then
write(6,'(A)')" +--- Root find, go to conversion phase --+ "
tx1=x1
ty1=y1
if(trim(method).eq."bisection")then
!bisection rule
do j=1,60
x2=0.5d0*(x0+x1)
y2=func(energy,r,vy,omgx,x2,omgz,u,zeroinx,exvz,ik)
if(y2.lt.0.d0.eqv.y0.lt.0.d0)then
x0=x2
else
x1=x2
endif
write(6,'(A,f16.8,A)')"Conversion => ",eps*abs(x2/(x1-x0))*100d0,"%"
if(abs((x1-x0)/x2).lt.eps)then
!If y2 is large, y2 is singular point.
if(y2.le.1d0)then
count=count+1
rd(count)=x2
endif
exit
endif
enddo
if(j.eq.60)write(6,*)"+--cannot convergence at root-finding method bs--+"
elseif(trim(method).eq."newton_raphson")then
!Newton-Raphson method
do j=1,10
h=1.d-4
y0=func(energy,r,vy,omgx,x0,omgz,u,zeroinx,exvz,ik)
Dy0=(func(energy,r,vy,omgx,x0+h,omgz,u,zeroinx,exvz,ik)-y0)/h
x2=x0-y0/Dy0
if(abs((x0-x2)/x2).lt.eps)then
if(y0.le.1d0)then
count=count+1
rd(count)=x2
endif
exit
endif
x0=x2
enddo
if(j.eq.10)write(6,*)"+--cannot convergence at root-finding method nr--+"
elseif(trim(method).eq."false_position")then
!false position method
if(abs(y1).gt.abs(y0))then
do j=1,30
x2=x0-y0*(x1-x0)/(y1-y0)
write(6,'(A,f16.8,A)')"Conversion => ",eps*abs(x2/(x2-x0))*100d0,"%"
if(abs((x2-x0)/x2).lt.eps)then
if(y0.le.1d0)then
count=count+1
rd(count)=x2
end if
exit
endif
x0=x2
y0=func(energy,r,vy,omgx,x0,omgz,u,zeroinx,exvz,ik)
enddo
else
do j=1,30
x2=x0-y0*(x1-x0)/(y1-y0)
write(6,'(A,f15.8,A)')"Conversion => ",eps*abs(x2/(x1-x2))*100d0,"%"
if(abs((x2-x1)/x2).lt.eps)then
if(y1.le.1d0)then
count=count+1
rd(count)=x2
endif
exit
endif
x1=x2
y1=func(energy,r,vy,omgx,x1,omgz,u,zeroinx,exvz,ik)
enddo
endif
if(j.eq.30)write(6,*)"+--cannot convergence at root-finding method fp--+"
else
write(6,'(A,A)')"unknown type of method, your method--> ",trim(method)
write(6,'(A,A)')"program stop"
stop
end if
if(count.ge.Nr)exit
x1=tx1
y1=ty1
endif
x0=x1
y0=y1
enddo

return
end subroutine rootfindingw
!---------------------------
function rkfd(N,t,x,s,ik)
use GBL
implicit none
integer,intent(in)::N,s,ik
double precision,intent(in)::t,x(1:N)
double precision::r(1:3),v(1:3),relv(1:3),omg(1:3),u(1:3),nw,nv,I
double precision::rkfd,gravity,vis1,vis2,mag,Nz,Nze
double precision,parameter::pi=datan(1d0)*4d0
external::gravity,vis1,vis2,mag,Nz,Nze

r(1:3)=x(1:3)
v(1:3)=x(4:6)
omg(1:3)=x(7:9)
u(1:3)=x(10:12)
relv(1:3)=v(1:3)-u(1:3)

rkfd=0.d0
if(s.le.3)then
! Differential equation of position
! d {x,y,z}/dt = v{x,y,z}
rkfd=v(s)
elseif(s.le.6)then
! Differential equation of velocity
! d v_{x,y,z}/dt = F{x,y,z}
rkfd=gravity(s-3,m,g,r)+vis1(s-3,relv,a,eta) &
+vis2(s-3,relv,a,eta,rho)+mag(s-3,omg,relv,a,rho)
rkfd=rkfd/m
elseif(s.le.9)then
! Differential equation for omega
! d omega{x,y,z}/dt = N_{x,y,z}/I
if(trim(omgdecay).eq."yes")then
I=0.4d0*m*a*a
nv=dsqrt(relv(1)**2+relv(2)**2+relv(3)**2)
nw=dsqrt(omg(1)**2+omg(2)**2+omg(3)**2)
if(nw.le.1d-13)then
rkfd=0d0
else
if(ik.eq.0)then
rkfd=(Nz(nv,a,eta,rho,nw)/I)*x(s)/nw
else
rkfd=(Nze(nv,a,eta,rho,nw)/I)*x(s)/nw
end if
endif
else
rkfd=0d0
endif
else
! Differential equation for wind
! d u{x,y,z}/dt = Fu_{x,y,z}
rkfd=0d0
endif

return
end function rkfd
!-----------------------------
function gravity(dir,m,g,r)
implicit none
integer,intent(in)::dir
double precision::gravity
double precision,intent(in)::m,g,r(1:3)

if(dir.eq.3)then
gravity=-m*g
else
gravity=0.d0
endif

return
end function gravity
!-------------------------------
function c1(a,eta)
implicit none
double precision::c1
double precision,intent(in)::eta,a
double precision,parameter::pi=dacos(-1.d0)

c1=6.d0*pi*eta*a

return
end function c1
!------------------------------
function c2(nv,a,eta,rho)
implicit none
double precision,parameter::pi=dacos(-1.d0)
double precision,intent(in)::nv,eta,a,rho
double precision::Cd,Reynolds,c2
external::Cd,Reynolds

c2=0.5d0*Cd(Reynolds(nv,a,eta,rho))*rho*pi*a*a

return
end function c2
!-----------------------------
function vis1(dir,v,a,eta)
implicit none
integer,intent(in)::dir
double precision,intent(in)::v(1:3),a,eta
double precision::vis1,norm,nv,c1
external::c1

nv=dsqrt(v(1)**2.d0+v(2)**2.d0+v(3)**2.d0)
norm=c1(a,eta)*nv

vis1=-norm*v(dir)/nv

return
end function vis1
!-------------------------
function vis2(dir,v,a,eta,rho)
implicit none
integer,intent(in)::dir
double precision,intent(in)::v(1:3),a,eta,rho
double precision::vis2,norm,nv,c2
external::c2

nv=dsqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
norm=c2(nv,a,eta,rho)*nv*nv

vis2=-norm*v(dir)/nv

return
end function vis2
!--------------------------
function mag(dir,omg,v,a,rho)
implicit none
integer,intent(in)::dir
double precision,intent(in)::omg(1:3),v(1:3),a,rho
double precision,parameter::pi=dacos(-1.d0)
double precision::mag,L(1:3),nomg,nv,nL,Cl

nomg=dsqrt(omg(1)*omg(1)+omg(2)*omg(2)+omg(3)*omg(3))
if(nomg.le.1.d-14)then
mag=0.d0
return
endif

nv=dsqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
L(1)=v(2)*omg(3)-v(3)*omg(2)
L(2)=v(3)*omg(1)-v(1)*omg(3)
L(3)=v(1)*omg(2)-v(2)*omg(1)
nL=dsqrt(L(1)*L(1)+L(2)*L(2)+L(3)*L(3))

if(nL.le.1.d-14)then
mag=0.d0
return
endif

! for BB bullet
Cl=0.12d0
mag=-Cl*(4.d0/3.d0)*pi*(a**3d0)*2.d0*rho*nomg*nv*L(dir)/nL

return
end function mag
!--------------------------
function Reynolds(nv,a,eta,rho)
!Reynolds number
! nv : norm of velocity of object
! a : radius of object
!eta : viscosity (not Kinetic viscosity)
!rho : density of fluid
implicit none
double precision,intent(in)::nv,a,eta,rho
double precision::keta,Reynolds

!keta means Kinetic viscosity
keta=eta/rho

Reynolds=nv*2.d0*a/keta

return
end function Reynolds
!---------------------------
function Cd(Re)
!From
!http://www.chem.mtu.edu/~fmorriso/DataCorrelationForSphereDrag2013.pdf
!Drag coefficient Cd,
!Cd depend on Reynolds number,Re.
!Fource of Drug,D is written by
! 1
! D= ---Cd*rho*pi*a**2*|V|**2
! 2
! ^ This Cd!
implicit none
double precision,intent(in)::Re
double precision::Cd,c1,c2,c3,c4

c1=24.d0/Re
c2=Re/5.d0
c2=2.6d0*c2/(1.d0+c2**(1.52d0))
c3=Re/263000d0
c3=0.411d0*c3**(-7.94d0)/(1.d0+c3**(-8d0))
c4=Re**0.8d0/461000d0

Cd=c1+c2+c3+c4

return
end function Cd
!----------------------------
function Nz(nv,a,eta,rho,omg)
!Moment of omg direction
implicit none
double precision,intent(in)::nv,a,eta,rho,omg
double precision::Nz,Cf,Fintegral
external::Cf,Fintegral

Nz=0.5d0*rho*Cf(nv,a,eta,rho)*a**3
Nz=Nz*Fintegral(nv,a,omg)

return
end function Nz
function Nze(nv,a,eta,rho,omg)
implicit none
double precision,intent(in)::nv,a,eta,rho,omg
double precision::Nze,Cf,Fintegral
double precision,parameter::pi=atan(1d0)*4d0
external::Cf,Fintegral
double precision::pc,tc,vu,vd,t

pc=pi/5.32065d0 ! magic phi
tc=pi/3.60475d0 ! magic theta

vu= nv*sin(pc)-a*omg*sin(tc)
vd=-nv*sin(pc)-a*omg*sin(tc)
Nze=-0.5d0*rho*Cf(nv,a,eta,rho)
Nze=Nze*(4d0*pi*a**2)*a*0.5d0
Nze=-Nze*(abs(vu)*vu+abs(vd)*vd)

return
end function Nze
!-------------------
function Cf(nv,a,eta,rho)
implicit none
! nv : norm of velocity of object
! a : radius of object
!eta : viscosity (not Kinetic viscosity)
!rho : density of fluid
double precision,intent(in)::nv,a,eta,rho
double precision::Reynolds,Cf
external::Reynolds

! Laminar flow
Cf=1.328d0/dsqrt(Reynolds(nv,a,eta,rho))

! Turbulent flow Re > 10^7
!Cf=0.455d0/(log10(Reynolds(nv,a,eta,rho))**(2.58d0))

return
end function Cf
!----------------------------------
function Fintegral(u,R,omega)
! 2016/07/23
! 2pi pi
! / /
! | dphi | |u*sin(phi)-R*omega*sin(theta)|*{u*sin(phi)-R*omega*sin(theta)}*sin^2(theta) d theta
! / /
! 0 0
implicit none
double precision,intent(in)::u,R,omega
double precision,parameter::pi=datan(1.d0)*4d0
double precision::Fintegral,s,Fphi,x(1:15),w(1:15)
external::Fphi
integer::i

s=0d0
call GaussKronrod15ab(0d0,pi,x,w)
do i=1,15
s=s+w(i)*Fphi(u,R,omega,x(i))
enddo
call GaussKronrod15ab(pi,2d0*pi,x,w)
do i=1,15
s=s+w(i)*Fphi(u,R,omega,x(i))
enddo
Fintegral=s

return
end function Fintegral
!------------------------------
subroutine GaussKronrod15ab(a,b,x,w)
!Gauss-Kronrod Quadrature Nodes and Weights
!http://www.advanpix.com/2011/11/07/gauss-kronrod-quadrature-nodes-weights/
implicit none
double precision,intent(in)::a,b
double precision,intent(out)::x(1:15),w(1:15)

integer::i
integer,parameter::N=15

x=0d0; w=0d0

x( 8) = 0d0
x( 9) = 2.077849550078984676006894037732449d-1
x(10) = 4.058451513773971669066064120769615d-1
x(11) = 5.860872354676911302941448382587296d-1
x(12) = 7.415311855993944398638647732807884d-1
x(13) = 8.648644233597690727897127886409262d-1
x(14) = 9.491079123427585245261896840478513d-1
x(15) = 9.914553711208126392068546975263285d-1
do i=1,7
x(i)=-x(N-i+1)
enddo

w( 8) = 2.094821410847278280129991748917143d-1
w( 9) = 2.044329400752988924141619992346491d-1
w(10) = 1.903505780647854099132564024210137d-1
w(11) = 1.690047266392679028265834265985503d-1
w(12) = 1.406532597155259187451895905102379d-1
w(13) = 1.047900103222501838398763225415180d-1
w(14) = 6.309209262997855329070066318920429d-2
w(15) = 2.293532201052922496373200805896959d-2
do i=1,7
w(i)=w(N-i+1)
enddo

x=0.5d0*((b-a)*x+(a+b))
w=0.5d0*(b-a)*w

return
end subroutine GaussKronrod15ab
!----------------------------------------
function Fphi(u,R,omega,phi)
! using analysis solution.
!
! pi
! /
! | |u*sin(phi)-R*omega*sin(theta)|*{u*sin(phi)-R*omega*sin(theta)}*sin^2(theta) d theta
! /
! 0
!
implicit none
double precision,intent(in)::u,R,omega,phi
double precision,parameter::pi=dacos(-1d0)
double precision::a,b,c,tp,Fphi,Rw

a=(u/(R*omega))*sin(phi)

if(a.le.0d0)then
Fphi=-(0.5d0*a*a*pi-8d0*a/3d0+3d0*pi/8d0)
elseif(a.ge.1d0)then
Fphi=0.5d0*a*a*pi-8d0*a/3d0+3d0*pi/8d0
elseif(a.gt.0d0.and.a.lt.1d0)then
tp=asin(a)
b=-a*a*pi*0.5d0-8d0*a/3d0-3d0*pi/8d0+(2d0*a*a+3d0/2d0)*tp
c=-(a*a+1d0)*dsin(2d0*tp)+dsin(4d0*tp)/8d0+6d0*a*dcos(tp)-2d0*a*dcos(3d0*tp)/3d0
Fphi=b+c
else
write(6,*)"undefined at Fphi, program stop",a
stop
endif

Rw=R*omega
Fphi=Fphi*Rw*abs(Rw)

return
end function Fphi

!---------------------------------
subroutine bullet_orbit(N,x0,Nt,time,filename,ik)
!sikinote http://slpr.sakura.ne.jp/qp/
!Author : sikino
!Date : 2016/03/18 (yyyy/mm/dd)
! : 2016/03/20
use GBL
use RKmod
implicit none
integer,intent(in)::ik,N,Nt
double precision,intent(in)::x0(1:N),time(0:Nt)
character(*),intent(in)::filename
double precision,parameter::pi=datan(1d0)*4d0
integer::info,i,j
double precision::x(1:N),Ene,h,t,tbound,rkfd
external::rkfd

x=x0
info=0

t=time(0)
open(22,file=trim(filename),status='replace')
write(22,'(A,e12.5e1,A)')"# m: ",m,"[kg]"
write(22,'(A,f12.5,A)')"# g: ",g,"[m s^{-2}]"
write(22,'(A,e12.5e1,A)')"# a: ",a,"[m]"
write(22,'(A,f12.5,A)')"# temperature: ",temperature,"[degree Celsius]"
write(22,'(A,f12.5,A)')"# moisture: ",moisture,"[no-dimension]"
write(22,'(A,f12.5,A)')"# pressure: ",pressure,"[Pa]"
write(22,'(A,e12.5e1,A)')"# eta:",eta,"[kg m^{-1} s^{-1}]"
write(22,'(A,f12.5,A)')"# rho:",rho,"[kg m^{-3}]"
write(22,'(A,e12.5e1)')"# tol:",tol
write(22,'(A,e12.5e1)')"# eps:",eps
write(22,'(A,14A13)')"#","t[s]","x[m]","y[m]","z[m]","vx[m/s]","vy[m/s]","vz[m/s]" &
,"wx[rot/s]","wy[rot/s]","wz[rot/s]","windx[m/s]","windy[m/s]","windz[m/s]","Energy[J]"
Ene=0.5d0*m*(x(4)**2+x(5)**2+x(6)**2)
write(22,'(14f13.6)')t,(x(i),i=1,6),(x(i)/(2d0*pi),i=7,9),(x(i),i=10,12),Ene

call rk_preparation("rkf45")
do j=0,Nt-1
t=time(j); tbound=time(j+1); h=tbound-t
call rkf451_e(rkfd,t,x,tbound,info,tol,ik)

Ene=0.5d0*m*(x(4)**2+x(5)**2+x(6)**2)
write(22,'(14f13.6)')t,(x(i),i=1,6),(x(i)/(2d0*pi),i=7,9),(x(i),i=10,12),Ene
!Stop calculation when bullet reach ground.
if(x(3).lt.0d0)exit
enddo
call rk_deallocation("rkf45")

close(22)
write(6,'(A,e10.3)')"Reach ground at ",t

return
end subroutine bullet_orbit

program main
use GBL
implicit none
integer::i,N,Nt,ik
double precision,allocatable::x(:),time(:),r(:),u(:)
double precision,parameter::pi=datan(1d0)*4d0
double precision::rx,ry,rz,vx,vy,vz,ux,uy,uz,omgx,omgy,omgz
double precision::energy,zeroinx,updownz,stept,exvz,exwy,theta
character(48)::outputfile,filename,search_zeroin,search_updown
double precision::tbound
double precision::rho_humid,eta_air
double precision::r0(1:3),v0(1:3),omg0(1:3),u0(1:3)
external::rho_humid,eta_air

double precision::xa,xb,za,zb
integer::nodes,poles
real::t1,t0

namelist /input/m,a,energy,g,temperature,pressure,moisture &
,omgdecay,search_zeroin,zeroinx,search_updown,updownz &
,theta,rx,ry,rz,vy,ux,uy,uz,omgx,omgy,omgz &
,stept,outputfile,ik

open(10,file="./input")
read(10,nml=input)
close(10)
write(*,nml=input)

N=12
eta=eta_air(temperature)
rho=rho_humid(temperature,pressure,moisture)

vx=dsqrt((2d0*energy/m)-vy**2)*dcos(pi*theta/180d0)
vz=dsqrt((2d0*energy/m)-vy**2)*dsin(pi*theta/180d0)
omgx=omgx*2*pi
omgy=omgy*2*pi
omgz=omgz*2*pi

allocate(x(1:N))
x(1)=rx; x(4)=vx; x(7)=omgx; x(10)=ux
x(2)=ry; x(5)=vy; x(8)=omgy; x(11)=uy
x(3)=rz; x(6)=vz; x(9)=omgz; x(12)=uz

tbound=100d0
Nt=nint(tbound/stept)
allocate(time(0:Nt)); time=0d0
do i=0,Nt
time(i)=dble(i)*stept
enddo

filename=trim(outputfile)//".txt"
call bullet_orbit(N,x,Nt,time,filename,ik)

call cpu_time(t0)
if(trim(search_zeroin).eq."yes")then
write(6,'(A)')"==============================="
write(6,'(A,f10.5)')" search vz and wy by zeroin --> ", zeroinx
allocate(r(1:3),u(1:3)); r(1:3)=x(1:3); u(1:3)=x(10:12)
call detwy(energy,r,x(5),x(7),x(9),u,zeroinx,exwy,exvz,ik)
write(6,*)exwy/(2d0*pi),exvz
x(4)=dsqrt((2d0*energy/m)-exvz**2-vy**2)
x(5)=vy
x(6)=exvz
x(8)=exwy
filename=trim(outputfile)//"_opt.txt"
call bullet_orbit(N,x,Nt,time,filename,ik)

deallocate(r,u)
endif

if(trim(search_updown).eq."yes")then
write(6,'(A)')"==============================="
write(6,'(A,f10.5)')" search vz and wy by up-down height --> ", updownz
allocate(r(1:3),u(1:3)); r(1:3)=x(1:3); u(1:3)=x(10:12)
call detwy2(energy,r,x(5),x(7),x(9),u,updownz,exwy,exvz,ik)
write(6,*)exwy/(2d0*pi),exvz
x(4)=dsqrt((2d0*energy/m)-exvz**2-vy**2)
x(5)=vy
x(6)=exvz
x(8)=exwy
filename=trim(outputfile)//"_h.txt"
call bullet_orbit(N,x,Nt,time,filename,ik)
deallocate(r,u)
endif

call cpu_time(t1)
write(6,'(f10.3,A)')(t1-t0),"[CPU sec]"

stop
end program main

使い方

・fortranコンパイラ(ifortやgfortran)
・gnuplot
を想定しています。

2017/02/04追)
Bernadotte66さんの報告より、windows上でgfortranコンパイラを用いる際、gfortranコンパイラがver5.1以前のものである場合エラーが出てしまうようです。ver6.2以上のコンパイラを推奨します。

Linuxでは、gfortranコンパイラver4.8.4で問題ないことは確かめています。

それでコンパイルして実行してください。
一連の流れとしては、以下の通りです。

$ gfortran main.f90
$ ./a.out
&INPUT
M= 2.5000000000000001E-004,
A= 3.0000000000000001E-003,
ENERGY= 0.90000000000000002 ,
G= 9.8066499999999994 ,
TEMPERATURE= 20.000000000000000 ,
PRESSURE= 101325.00000000000 ,
MOISTURE= 0.59999999999999998 ,
OMGDECAY="yes ",
SEARCH_ZEROIN="yes ",
ZEROINX= 50.000000000000000 ,
SEARCH_UPDOWN="yes ",
UPDOWNZ= 0.20000000000000001 ,
THETA=-0.50000000000000000 ,
RX= 0.0000000000000000 ,
RY= 0.0000000000000000 ,
RZ= 1.0000000000000000 ,
VY= 0.0000000000000000 ,
UX= 0.0000000000000000 ,
UY= 0.0000000000000000 ,
UZ= 0.0000000000000000 ,
OMGX= 0.0000000000000000 ,
OMGY= -210.00000000000000 ,
OMGZ= 0.0000000000000000 ,
STEPT= 1.0000000000000000E-002,
OUTPUTFILE="orbit ",
IK= 0,
/
Reach ground at 0.980E+00
===============================
search vz and wy by zeroin --> 50.00000
---- false_position
Sequence 1 of 25 : omgy/2pi -767.96000000 : zeroin 100.00000000
Sequence 2 of 25 : omgy/2pi -735.92000000 : zeroin 100.00000000
Sequence 3 of 25 : omgy/2pi -703.88000000 : zeroin 100.00000000
Sequence 4 of 25 : omgy/2pi -671.84000000 : zeroin 100.00000000
Sequence 5 of 25 : omgy/2pi -639.80000000 : zeroin 100.00000000
Sequence 6 of 25 : omgy/2pi -607.76000000 : zeroin 100.00000000
Sequence 7 of 25 : omgy/2pi -575.72000000 : zeroin 17.22935680
Sequence 8 of 25 : omgy/2pi -543.68000000 : zeroin 15.16813299
Sequence 9 of 25 : omgy/2pi -511.64000000 : zeroin 12.90696782
Sequence 10 of 25 : omgy/2pi -479.60000000 : zeroin 10.42356698
Sequence 11 of 25 : omgy/2pi -447.56000000 : zeroin 7.68923851
Sequence 12 of 25 : omgy/2pi -415.52000000 : zeroin 4.66717441
Sequence 13 of 25 : omgy/2pi -383.48000000 : zeroin 1.30910154
Sequence 14 of 25 : omgy/2pi -351.44000000 : zeroin -2.44916818
+--- Root find, go to conversion phase --+
Conversion => 0.00033361%
Conversion => 0.00902931%
Conversion => 0.24389215%
Conversion => 6.77148580%
Conversion => 104.78145034%
---- false_position
Sequence 1 of 1 : omgy/2pi -366.89192543 : zeroin -0.57095843
+--- Root find, go to conversion phase --+
Conversion => 0.00075051%
Conversion => 0.08150935%
Conversion => 8.77790138%
Conversion => 579.19923386%
-371.80132940800553 -2.5094934272766087
Reach ground at 0.170E+01
===============================
search vz and wy by up-down height --> 0.20000
---- false_position
Sequence 1 of 25 : omgy/2pi -959.96000000 : zeroin 19.21215678
Sequence 2 of 25 : omgy/2pi -919.92000000 : zeroin 17.82618788
Sequence 3 of 25 : omgy/2pi -879.88000000 : zeroin 16.40387340
Sequence 4 of 25 : omgy/2pi -839.84000000 : zeroin 14.95914779
Sequence 5 of 25 : omgy/2pi -799.80000000 : zeroin 9.80000000
Sequence 6 of 25 : omgy/2pi -759.76000000 : zeroin 9.80000000
Sequence 7 of 25 : omgy/2pi -719.72000000 : zeroin 9.80000000
Sequence 8 of 25 : omgy/2pi -679.68000000 : zeroin 9.80000000
Sequence 9 of 25 : omgy/2pi -639.64000000 : zeroin 9.80000000
Sequence 10 of 25 : omgy/2pi -599.60000000 : zeroin 9.80000000
Sequence 11 of 25 : omgy/2pi -559.56000000 : zeroin 1.45287636
Sequence 12 of 25 : omgy/2pi -519.52000000 : zeroin 1.10739833
Sequence 13 of 25 : omgy/2pi -479.48000000 : zeroin 0.80320434
Sequence 14 of 25 : omgy/2pi -439.44000000 : zeroin 0.54089117
Sequence 15 of 25 : omgy/2pi -399.40000000 : zeroin 0.32062896
Sequence 16 of 25 : omgy/2pi -359.36000000 : zeroin 0.14208536
Sequence 17 of 25 : omgy/2pi -319.32000000 : zeroin 0.00427472
Sequence 18 of 25 : omgy/2pi -279.28000000 : zeroin -0.09464709
+--- Root find, go to conversion phase --+
Conversion => 0.00183551%
Conversion => 0.00989716%
Conversion => 0.05297720%
Conversion => 0.28399231%
Conversion => 1.52644398%
Conversion => 8.26417579%
Conversion => 39.67268127%
Conversion => 179.58494188%
---- false_position
Sequence 1 of 1 : omgy/2pi -312.86039405 : zeroin -0.01405310
+--- Root find, go to conversion phase --+
Conversion => 0.00065868%
Conversion => 0.03115036%
Conversion => 1.49063057%
Conversion => 102.25144566%
-317.78761751909860 -1.5473682498931858
Reach ground at 0.149E+01
2.681[CPU sec]
$ gnuplot

G N U P L O T
Version 4.6 patchlevel 4 last modified 2013-10-02
Build System: Linux x86_64

Copyright (C) 1986-1993, 1998, 2004, 2007-2013
Thomas Williams, Colin Kelley and many others

gnuplot home: http://www.gnuplot.info
faq, bugs, etc: type "help FAQ"
immediate help: type "help" (plot window: hit 'h')

Terminal type set to 'wxt'
gnuplot> plot "orbit.txt" u 2:4 w l
gnuplot> replot "orbit_h.txt" u 2:4 w l
gnuplot> replot "orbit_opt.txt" u 2:4 w l

gnuplot上での出力が

のようになったら成功です。

出力ファイルは最大で3つ生成(inputパラメータ “search_??????” )に依存します。
orbit.txtは、inputによって与えられたパラメータでそのまま計算したもの
orbit_opt.txtはそのエネルギー、重さで上下振れ幅が最小でゼロイン位置がinputのzeroinxのもの
orbit_h.txtはそのエネルギー、重さで上下振れ幅が最小で上下振れ幅がinputのupdownzのもの
となります。

出力ファイルはそれぞれ14列あり、
時刻t[秒],位置x,y,z[m]:x,y,z方向の速度[m/s]:BB弾の回転数x,y,z[回転/s]:x,y,z方向の風速:エネルギー[J]
になっています。
出力ファイルの中身は下のように書いてあるので、見てください。

# m: 0.25000E-3[kg]
# g: 9.80665[m s^{-2}]
# a: 0.30000E-2[m]
# temperature: 20.00000[degree Celsius]
# moisture: 0.60000[no-dimension]
# pressure: 101325.00000[Pa]
# eta: 0.18197E-4[kg m^{-1} s^{-1}]
# rho: 1.20123[kg m^{-3}]
# tol: 0.10000E-7
# eps: 0.10000E-8
# t[s] x[m] y[m] z[m] vx[m/s] vy[m/s] vz[m/s] wx[rot/s] wy[rot/s] wz[rot/s] windx[m/s] windy[m/s] windz[m/s] Energy[J]
0.000000 0.000000 0.000000 1.000000 84.849583 0.000000 -0.740471 0.000000 -210.000000 0.000000 0.000000 0.000000 0.000000 0.900000
0.010000 0.838863 0.000000 0.992911 82.937695 0.000000 -0.678265 0.000000 -209.261019 0.000000 0.000000 0.000000 0.000000 0.859890
...

実銃のパラメータでの計算

実銃は特に想定していませんが、計算ができると言えばできます。
実銃と比較して考慮されていない点は3つあり、

滑らかな球として考えているので抗力係数$C_d$が異なる
銃弾の歳差運動が考慮されていない(球ではなく細長いので)
ライフリングによる回転が異なる
コリオリ力を考慮に入れていない

1番目はかなり重大です。
私のプログラムを使うと1000メートル地点でのエネルギーの値が2倍ほど大きく出ます。
ライフリングによる横回転だけなら問題なく考慮できます。バックスピンで入っていたものを変えればいいだけです。
ただし、風が無い、またはコリオリ力を考慮しない状況では計算上は横回転があっても無くても変わりません。
風とコリオリ力がある時にライフリングによる弾道のずれが影響します。
2番目はどんな力が働いたりするか私には見当がつきません。専門の方じゃないと分かるまでに時間がかかると思います。
コリオリ力はプログラムで入れればいいです。これはそこまで辛くないでしょう。
ライフリングの効果は入れられるには入れられますが、ちょっと変な結果になるので使わないほうが良いでしょう。
なのでこの効果は入れていません。
パッと思いつく限り、違う点は以上です。

よって、本プログラムによる実銃を想定した計算とは、
風が無く、赤道上で実銃と同じ口径の滑らかな球が実銃のエネルギーで射出されたときの弾道
となります。

ちなみに、.338ラプアマグナム弾の場合を想定したinputパラメータはこちら
角度を使って計算しているのでそこのとこのプログラムは各自直してください。

&input
! For .338 lapua magnum( but treated as sphire)
m=16.2d-3, ! Weight of bullet [kg]
a=4.29d-3, ! Radious of bullet [m]
energy=6562.16d0, ! Energy of bullet [J]
g=9.80665d0, ! Acceleration of Gravity [m s^-2]
temperature=20d0, ! Temperature [degree of Celsius]
pressure=101325d0, ! Atmospheric pressure [Pa]
moisture=0.6d0, ! Value of moisture 0~1
omgdecay="no", ! Introduce omega decay. "yes" or "no"
search_zeroin="no", ! Search good vz and omgy. "yes" or "no"
zeroinx=50d0, ! +-- if yes, set zeroin length [m]
search_updown="no",! Search good vz and omgy. "yes" or "no"
updownz=0.10d0, ! +-- if yes, set up-down henght [m]
theta=0.6d0, ! default, theta is neglected. change program main.
rx=0.d0,
ry=0.d0,
rz=1.0d0,
vy=0.d0,
vz=0.d0,
ux=0.d0,
uy=0.d0,
uz=0.d0,
omgx=0.d0, !*2pi [omgx*2*pi rotate/s]
omgy=0.d0, !*2pi [omgy*2*pi rotate/s]
omgz=0.d0, !*2pi
outputfile="orbit", ! Output file name prefix
stept=0.01d0, ! Output interval [s]
&end

上記パラメータで計算するとこんな軌道が得られます(点は0.05秒ごとの軌跡)

.338ラプアマグナム弾のデータは↓
.338 ラプアマグナム弾 / .338 Lapua Magnum -MEDIAGUN DATABASE
を参照しました。

改良版

これから最適な回転量、角度を探すプログラムに更新したいと思いますが、とりあえず
改良したプログラムを載せておきます。
input
main.f90

プログラミングと数値計算

対角化

2015年7月4日 Fernandeskun コメントする

fortran90でLapackを用います。

倍精度ルーチン

・複素一般行列の対角化, diag(N,A,Ev)
・複素エルミート行列の対角化, zdiag(N,A,Ev)
・実エルミート行列の対角化, herdiag(N,A,Ev), diag_slct(N,A,Nq,Ev)
・実対称三重対角行列の対角化, tridiag(N,A,Ev)

4倍精度ルーチン

・複素一般行列の対角化, qdiag(N,A,Ev)

複素一般行列　　　　　　　　　：行列要素は複素数で、値は適当に入っている。
複素エルミート行列　　　　　　：行列要素は複素数で、$A_{i,j}=A^*_{j,i}$の関係がある。
実エルミート行列（実対称行列）：行列要素は実数で、$A_{i,j}=A_{j,i}$の関係がある。
実三重対角行列　　　　　　　　：行列要素は実数で、$A_{i,i}\ne 0, ~~A_{i,i-1}=A_{i,i+1}$、それ以外の要素は0。

（エルミートな行列と言われたら、その固有値は必ず実数です。）

複素一般行列の場合

wolframのページで、
Eigenvectors[{-2,2,4},{-2,4,2},{-2,1,4}]
と入力することで固有値、固有ベクトルを出してくれます。これで確かめを行いました。

・この行列を入れた時、正しい固有値と固有ベクトルを出力してくれる事を確認しています。
・40×40程度の行列でも確認しています。
・複素数の固有値、固有ベクトルでも確認しています。

Lapackを使います。

ifort -mkl main.f90

か、g95, もしくはgfortranでlapackにリンク

gfortran -llapack main.f90

して、コンパイルしてください。

大きさN×Nの正方行列Aを入力すると, 固有値Ev(1:N)とそれぞれの固有値の要素に属する固有ベクトルを上書きして返します。(1列目に固有値Ev(1)に属する固有ベクトルが、2列目に固有値Ev(2)に属する固有ベクトルが･･･という感じに。)
また固有値は実部の大小を比較して小さい順に並べ替えられ、固有ベクトルもそれに対応した列に格納されます。
計算速度は$O(N^3)$です。

▼ここクリックで展開

program main
implicit none
integer::i,j,N
complex(kind(0d0)),allocatable::A(:,:),Ev(:)

N=3
allocate(A(1:N,1:N),Ev(1:N))
A(1:N,1:N)=dcmplx(0d0,0d0)
Ev(1:N)=dcmplx(0d0,0d0)

A(1,1:N)=(/(-2d0,0d0),(2d0,0d0),(4d0,0d0)/)
A(2,1:N)=(/(-2d0,0d0),(4d0,0d0),(2d0,0d0)/)
A(3,1:N)=(/(-2d0,0d0),(1d0,0d0),(4d0,0d0)/)

call diag(N,A,Ev)

do i=1,N
write(6,'(A,e12.5e2,A,e12.5e2,A)')"(",dble(Ev(i)),",",dimag(Ev(i)),")"
enddo

do j=1,N
do i=1,N
write(6,'(A,e12.5e2,A,e12.5e2,A,$)')"(",dble(A(j,i)),",",dimag(A(j,i)),") "
enddo
write(6,*)
enddo

stop
end program main
!---------------------

subroutine diag(N,A,Ev)
! sikinote
!date : 2015/07/07
! 2015/08/21
!developer : sikino & fernandeskun
implicit none
integer,intent(in)::N
complex(kind(0d0)),intent(inout)::A(1:N,1:N)
complex(kind(0d0)),intent(out)::Ev(1:N)

integer::ilo,ihi,info,lwork,turn(1:N),tmp,i
double precision::scale(1:N),rwork(1:N)
complex(kind(0d0))::tau(1:N-1),w(1:N),z(1:N,1:N),Q(1:N,1:N),vr(1:N,1:N),tw(1:3)
complex(kind(0d0)),allocatable::work(:)

tau(1:N-1)=dcmplx(0d0,0d0)
w(1:N)=dcmplx(0d0,0d0)
z(1:N,1:N)=dcmplx(0d0,0d0)
Q(1:N,1:N)=dcmplx(0d0,0d0)
vr(1:N,1:N)=dcmplx(0d0,0d0)
tw(1:3)=dcmplx(0d0,0d0)
Ev(1:N)=dcmplx(0d0,0d0)

!Equilibrate matrix A to equilibrated matrix A' to improve accuracy.
! i i io i o o o o
call zgebal('P', N, A, N, ilo, ihi, scale, info)
if(info.ne.0)then
write(6,'(A,i0)')" At zgebal error, info --> ",info
write(6,'(A)')" Program stop"
stop
endif

!Size Query
call zgehrd(N, ilo, ihi, A, N, tau, tw, -1, info)
lwork=nint(dble(tw(1)))
allocate(work(1:lwork)); work=dcmplx(0d0,0d0)

!Degenerate matrix A to upper Hessenberg matrix H.
! i i i io i o i i o
call zgehrd(N, ilo, ihi, A, N, tau, work, lwork, info)
if(info.ne.0)then
write(6,'(A,i0)')" At zgehrd error, info --> ",info
write(6,'(A)')" Program stop"
stop
endif
deallocate(work)

Q=a
!Size Query
call zunghr(N, ilo, ihi, Q, N, tau, tw, -1, info)
lwork=nint(dble(tw(1)))
allocate(work(1:lwork)); work=dcmplx(0d0,0d0)

!Make complex unitary matrix Q from upper Hessenberg matrix H.
! i i i io i i i i o
call zunghr(N, ilo, ihi, Q, N, tau, work, lwork, info)
if(info.ne.0)then
write(6,'(A,i0)')" At zunghr error, info --> ",info
write(6,'(A)')" Program stop"
stop
endif
deallocate(work)

z=Q
!Size Query
call zhseqr('S', 'V', N, ilo, ihi, A, N, Ev, z, N, tw, -1, info)
lwork=nint(dble(tw(1)))
allocate(work(1:lwork)); work=dcmplx(0d0,0d0)

!Get eigenvalue of upper Hessenberg matrix H and Get Schur vector.
! i i i io i o o i i i o
call zhseqr('S', 'V', N, ilo, ihi, A, N, Ev, z, N, work, lwork, info)
if(info.ne.0)then
write(6,'(A,i0)')" At zhseqr error, info --> ",info
write(6,'(A)')" Program stop"
stop
endif
deallocate(work)

!Get right eigenvector X from upper triangular matrix T.
allocate(work(1:2*N))
vr=z
! i i i o i i o i i i
call ztrevc('R', 'B', 0, N, A, N, 0, 1, vr, N, N, tmp, work, rwork, info)
if(info.ne.0)then
write(6,'(A,i0)')" At zhseqr error, info --> ",info
write(6,'(A)')" Program stop"
stop
endif
deallocate(work)

!Transrate right eigenvector X of Equilibrated matrix A' to right eigenvector of matrix A
! i i i i i o i o
call zgebak('P', 'R', N, ilo, ihi, scale, N, vr, N, info)
if(info.ne.0)then
write(6,'(A,i0)')" At zhseqr error, info --> ",info
write(6,'(A)')" Program stop"
stop
endif

A=vr

!swap Eigenvectol as same arrangement for Eigenvalue
call sortdp2(N,Ev,turn)

Q=A
do i=1,N
tmp=turn(i)
A(1:N,i)=Q(1:N,tmp)
enddo
return

!sort Eigenvalue of real part from small to big.
contains
subroutine sortdp2(N,data,turn)
implicit none
integer::i,ti,j,N,turn(1:N)
complex(kind(0d0))::data(1:N),tmp

do i=1,N
turn(i)=i
enddo

do i=1,N-1
do j=i+1,N
if(dble(data(i)) > dble(data(j)))then
tmp=data(i)
data(i)=data(j)
data(j)=tmp

ti=turn(i)
turn(i)=turn(j)
turn(j)=ti
end if
end do
end do

return
end subroutine sortdp2
end subroutine diag

複素エルミート行列の対角化

対称行列に限る場合は計算を早くすることができます。
ここでは複素エルミート行列を対角化して固有値を出力するケースを考えましょう。

そのまま使うのはやはり面倒です。なぜなら、このルーチンを用いるためにはワーク配列を何個か宣言して用いる必要があるからです。
このワーク配列は計算に直接関係ありません。
僕らが欲しいのは行列を入れたら固有値と固有ベクトルを返してくれさえすればいいものです。
下はサブルーチンzdiagは僕が作ったもので、複素エルミート行列$A(1:N,1:N)$を入れると固有値の配列$Ev(1:N)$と$A$に上書きして固有ベクトル$A(1:N,1:N)$を入れて返します。

エルミート行列なので、実は全ての行列要素を求める必要ありません。
下の例では、
エルミート行列$A(i,j)$は,$j\ge i$の上三角の要素だけ分かっていればいいです。
なので、コメントアウトしている行列でも同じ固有値と固有ベクトルが得られます。
以下のプログラムをlapackと共にコンパイルし、実行すると、

>./a.out
(-0.42419E+01)
( 0.21580E+01)
( 0.80839E+01)
(-0.88841E+00,-0.92008E-01) (-0.15371E+00,-0.13692E+00) (-0.39557E+00, 0.58357E-01)
( 0.17660E+00, 0.14837E-01) ( 0.55488E+00,-0.56209E+00) (-0.35177E+00, 0.47012E+00)
( 0.41334E+00, 0.00000E+00) (-0.57775E+00, 0.00000E+00) (-0.70382E+00, 0.00000E+00)
>

という結果が得られます。
1番目の固有値(-0.42419E+01)に対応する固有ベクトルは1列目に、
2番目の固有値( 0.21580E+01)に対応する固有ベクトルは2列目…というように代入されていきます。

プログラム例)
確かめ用のリンクEigenvectors[{-2,2+i,4},{2-i,4,1-2i},{4,1+2i,4}] -wolfram alpha

▼ここクリックでこの場に展開

行と列を入れ替えたければ、以下のルーチンを使用してください。

subroutine swap2d(N,A)
implicit none
integer::N,i,j
double precision::temp,A(1:N,1:N)

do i=1,N
do j=i+1,N
temp=A(i,j)
A(i,j)=A(j,i)
A(j,i)=temp
enddo
enddo

return
end subroutine swap2d

実エルミート行列の場合

以下のサブルーチンを使うと良いでしょう。

▼ここクリックでこの場に展開

$N\times N$の実対称行列で、固有値の小さい方から順に$Nq$個だけ欲しい場合

▼ここクリックでこの場に展開

実エルミート行列の、
完全対角化(herdiag)と、
選択的に固有値を計算する(diag_slct)
の比較(実時間計測)を行いました。

横軸が選択的に何個固有値を計算するか、を表し、縦軸は実際に掛かった時間をあらわします。
色は行列サイズによって変化し、黄色,紫、緑の順に行列サイズが$2^{16}\times 2^{16}, 2^{14}\times 2^{14}, 2^{12}\times 2^{12}$を表します。
破線はherdiagによる対角化時間、実線はdiag_slctによる計算時間を表します。
また、図中の値2,4,8は計算に用いたCPUの数です。CPUは8コア16スレッドです。

ここから分かるのは、欲しい固有値の数が行列サイズ$N\times N$の$N/4$より小さい場合は選択的に計算した方が早い、ということです。
選択的に計算する方は計算時に要求されるメモリも少なくて済みます。

(2016/08/27)

4倍精度の一般行列の対角化

4倍精度用の一般行列(複素非エルミートを含む)を対角化し、その固有値と固有ベクトルを出力するfortranコードが以下のページ
４倍精度対角化 -ただの備忘録
で公開されていました。
ここで公開されているサブルーチン”qeispack.f90”
を使うことで対角化できます。
gfortranとifortで動くことは確かめました。
コンパイルは例えば

ifort qeispack.f90 main.f90

のようにしてください。
本ページと同じように書き加えて実装すると以下のようになります。

▼ここクリックでこの場に展開

program main
implicit none
integer::i,j,N
real*16::qS,theta
complex*32,allocatable::qA(:,:),qEv(:)

N=3
allocate(qA(1:N,1:N),qEv(1:N))
qA(1:N,1:N)=cmplx(0q0,0q0,kind=16)
qEv(1:N)=cmplx(0q0,0q0,kind=16)

qA(1,1:N)=(/(-2q0,0q0),(2q0,0q0),(4q0,0q0)/)
qA(2,1:N)=(/(-2q0,0q0),(4q0,0q0),(2q0,0q0)/)
qA(3,1:N)=(/(-2q0,0q0),(1q0,0q0),(4q0,1q0)/)

call qdiag(N,qA,qEv)

do i=1,N
write(6,'(A,e45.33e2,A,e45.33e2,A)')"(",real(qEv(i)),",",imag(qEv(i)),")"
enddo

do j=1,N
do i=1,N
write(6,'(A,e45.33e2,A,e45.33e2,A,$)')"(",real(qA(j,i)),",",imag(qA(j,i)),") "
enddo
write(6,*)
enddo

stop
end program main
!---------------------
subroutine qdiag(N,A,Ev)
implicit none
integer,intent(in)::N
complex*32,intent(inout)::A(1:N,1:N)
complex*32,intent(out)::Ev(1:N)
!
! Date : 2016/08/27
! Author : sikino
!
! For the eigenvalues and eigenvectors of a complex general matrix.
! Use qeispack.f90 at "https://github.com/hanada-yasutaka/QEISPACK"
!
integer::matz,ierr,i,j,turn(1:N),tmp
real*16,dimension(1:N)::wi,wr
real*16,dimension(1:N,1:N)::ai,ar,xi,xr
complex*32::Q(1:N,1:N)

ar=real(A)
ai=imag(A)
!
! matz = 0 for eigenvalues only,
! matz = 1 for eigenvalues and eigenvectors.
!
matz = 1
call cg ( N, ar, ai, wr, wi, matz, xr, xi, ierr )

if ( ierr .ne. 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'qdiag - Warning!'
write ( *, '(a,i8)' ) ' The error return flag IERR = ', ierr
return
end if

A=cmplx(xr,xi,kind=16)
Ev=cmplx(wr,wi,kind=16)

!swap Eigenvectol as same arrangement for Eigenvalue
call sortqp2(N,Ev,turn)
Q=A
do i=1,N
tmp=turn(i)
A(1:N,i)=Q(1:N,tmp)
enddo
return

!sort Eigenvalue of real part from small to big.
contains
subroutine sortqp2(N,data,turn)
implicit none
integer::i,ti,j,N,turn(1:N)
complex*32::data(1:N),tmp

do i=1,N
turn(i)=i
enddo

do i=1,N-1
do j=i+1,N
if(real(data(i)) > real(data(j)))then
tmp=data(i)
data(i)=data(j)
data(j)=tmp

ti=turn(i)
turn(i)=turn(j)
turn(j)=ti
end if
end do
end do

return
end subroutine sortqp2
end subroutine qdiag

実対称三重対角行列の場合

実対称三角行列の場合も置いておきます。

▼ここクリックでこの場に展開

プログラミングと数値計算, ルンゲ=クッタ法

ルンゲ=クッタ法の説明と刻み幅制御

2015年4月19日 sikino 5件のコメント

ルンゲ=クッタ法(Runge-Kutta method、RK法)とは?
僕の知る限りの知識で紹介します。

特に良く使われる陽的ルンゲ=クッタ法は、
・実装が簡単
・良いアルゴリズムではない
という手法です。

良いアルゴリズムである陰的ルンゲ=クッタ法は、
陰的ルンゲ=クッタ法
をご覧ください。

埋め込まれた陽的ルンゲ=クッタ法
ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御(適応刻み幅制御)
ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御のプログラム(実数/複素数で１/２階微分方程式を解くプログラム)
不連続な点を含む場合
刻み幅制御のベンチマーク(振り子)
（追記）ルンゲ=クッタ=ドルマンド=プリンス法
陽的ルンゲ=クッタ法の導出
参考文献

理論はいいから4次ルンゲ=クッタ法の数値計算ではどうやるかだけ教えてくれ！という方は4次ルンゲ=クッタ法へどうぞ。

ルンゲ=クッタ法の系統

ルンゲ=クッタ法は微分方程式の数値計算解を得るための手法のことを指します。

通常の参考書で微分方程式を解くために良く紹介されているのは、オイラー法と中点法、４次ルンゲ=クッタ法でしょうか。
オイラー法も中点法も4次ルンゲ=クッタ法も、”陽的ルンゲクッタ法” と呼ばれる枠組みの1つです。

オイラー法は正確には “陽的1段1次ルンゲ=クッタ法” と呼ばれ、
中点法は “陽的2段2次ルンゲ=クッタ法”、
4次ルンゲクッタ法(RK4)は、”陽的4段4次ルンゲ=クッタ法” と呼ばれています。

“段”と”次”とはなんなんでしょう？それは、
計算の大変さ(段)と、計算の正確さ(次)
です。”段“の値が小さければ小さいほど計算時間が少なくて済みますし、”次“の値が高ければ高いほど計算が正確です。

オイラー法は1という計算コストで正確さ1が得られますし、RK4は4という計算コストで正確さ4が得られます。

4次ルンゲ=クッタ法が使われる理由

理由は実装が簡単でそれなりの精度を持つから。です。

陽的ルンゲ=クッタ法において、pという計算コスト(p段)で、pより大きな正確さq(q次)を得ることはできません。
Derivation of Runge–Kutta methodsによれば、
$q$次の正確さ(q次のオーダー)を得たい場合、最低限必要な段数$p_{\mbox{min}}(q)$は

という関係にあります。

ここで注目するべきは4次の時までは計算コストに比例して計算精度が上がっていきます。
しかし、4次以上では、計算コストの増加と計算の正確さが見合わなくなっていきます。故に
計算効率が良いのは4次だろう、
と予想できます。
また、重要な理由として、4次ルンゲ=クッタ法に現れる係数が5次以降と比べて圧倒的にシンプルであることが挙げられます。4次では$0,1/2,1/6$程度の係数だけが使われ、プログラム作成時の入力ミスがほぼ生まれません。しかし、5次では$28561/56430, -7200/2197$といった係数が数多く出てきます。
これらの理由から,4次ルンゲ=クッタ法(RK4)が数値計算科学の世界でよく使われるのです。

陽的ルンゲ=クッタ法に限って言えばプログラムの実装が非常に簡単であることが挙げられます。陰的ルンゲ=クッタ法と呼ばれるアルゴリズムもあり、これは陽的ルンゲ=クッタよりも優れていますが、計算量が多くなり、若干複雑なアルゴリズムになります。陰的ルンゲ=クッタ法を詳しく知りたい方は陰的ルンゲ=クッタ法をご覧ください。

Q. オイラー法もものすごく細かい分点を取れば、その計算効率はRK4と同じなんじゃないの？
A. 刻み幅の乗数で効いてくるのでそうではありません。高次を使っても計算が信頼できるのであれば、大きなステップで進めるほうが早いです。例え、目標精度への計算時間が同じだとしても、計算機の有効桁数によって否定されてしまいます。
RK4で典型的にとられる時間ステップの間隔は、おおよそ$10^{-2}\sim 10^{-4}$程度であり、RK4のエラーのオーダーは$O(h^5)$です。
そして、科学計算で使う際の有効桁数は倍精度型で16桁です。
1ステップ当りの誤差は$h$の減少に伴い、解が$h^4$の早さで収束していく、と言えます。
だから16桁の計算では$h=10^{-1}\to h=10^{-4}$に変化させる時、誤差は$O(h^5)=c 10^{-5}\to c 10^{-20}$と変化します。
おおよそ$c\approx 1$と見積もれば、(有効桁数16桁を多少超えてしまいますが、)有効桁数いっぱいまで正しい値が出るであろうと期待できます。

これに対し、オイラー法で同じことをするには$h$を$10^{-16}$にしなくてはいけません。
$t$の値が$10^{-16}$変わった時に、桁落ちの問題を回避できるほど関数$x$の値に変化が生まれるか？
が問題になり、まぁそんな急激な変化は生まれないでしょう、と予想できます。これでは桁落ちの問題を回避するほどの変化は到底望めません。

よって計算の効率と有効桁数の限界から、RK4なのです。

また、あまりに高次の方法を使うとルンゲ現象に代表される不安定性といったことが起こるかもしれません。
高次は高精度という意味とイコールではないことに注意しましょう。この高次≠高精度については等間隔の分点における数値積分の時に書いたので気になる方はどうぞ。。

誤差について

4次ルンゲ=クッタ法の、1ステップ当りの誤差は$h^5$に比例,すなわち$O(h^5)$です。
しかし、通常は細かいルンゲ=クッタ法を何度も繰り返して計算します。
区間$[a,b]$を刻み幅$h$の4次ルンゲ=クッタ法で$N$回のステップを繰り返し計算する場合、誤差は
$
\displaystyle N\cdot O(h^5)=\frac{b-a}{h}\cdot O(h^5)=O(h^4)$

となります。よって、$N$回繰り返すような計算では、オーダーが１つ落ちる事に注意しましょう。

[adsense1]

陽的ルンゲ=クッタ法の段数と次数について

さて、ここまで“段”は計算コスト、で“次”は計算の正確さ、という曖昧な表現でしたが、その表現をちゃんと知りましょう。
段と次を知るためにはルンゲ=クッタ法の計算方法を知る必要があります。
具体例を載せます。
$
\displaystyle \frac{dx}{dt}=f(t,x)
$
の、$t_{n}\rightarrow t_{n}+h\ \ (=t_{n+1})$における初期値問題に対する、
陽的1段1次ルンゲ=クッタ法(オイラー法)の計算スキームは、
$
x_{n+1}=x_{n}+h\cdot f(t_{n},x_{n})
$
です。

陽的4段4次ルンゲ=クッタ法(RK4)の計算スキームは、
$
\begin{align}
k_1&=f(t_n, x_n) \\
k_2&=f(t_n+h/2, x_n+h k_1/2) \\
k_3&=f(t_n+h/2, x_n+h k_2/2) \\
k_4&=f(t_n+h, x_n+h k_3) \\
x_{n+1}&=x_{n}+{(k_1+2k_2+2k_3+k_4)}h/6
\end{align}
$
として与えられます。

一般的に、陽的s段のルンゲ=クッタ法とは
$
\begin{align}
g_i&=x_n+h\sum_{j=1}a_{i,j}k_j\ \ \ (j\lt i, \ i=1,2,…,s) \\
k_i&=f(t_n+c_ih, g_i) \\
x_{n+1}&=x_n+h\sum_{i=1}^s b_ik_i
\end{align}
$
として書けます。
ここで行列形式で与えられる係数$a_{i,j}, b_{i},c_{i}$によって、そのs段ルンゲ=クッタ法が持つ次数が決められます。段数はここから由来します。

点$(t_n, x(t_n))$周りで関数をテーラー展開し、その関数が点$(t_n+h\ \ (=t_{n+1}), x(t_{n+1}))$で作る点を近似解とするのがルンゲ=クッタ法です。
故に、$x(t_{n+1})$は、
$
\begin{align}
x(t_{n+1})=x(t_n)+\left.\frac{h}{1!}\frac{dx}{dt}\right|_{t=t_n}+\left.\frac{h^2}{2!}\frac{d^2x}{dt^2}\right|_{t=t_n}+\left.\frac{h^3}{3!}\frac{d^3x}{dt^3}\right|_{t=t_n}+\left.\frac{h^4}{4!}\frac{d^4x}{dt^4}\right|_{t=t_n}+…
\end{align}
$
と書けます。
ここで、テイラー展開としてどの程度一致させて$x(t_n+h)$を決定するか？を表すのが次数に当たります。

言葉で書くなら、

局所的な精度は，高次項が解のテイラー(Taylor)展開とどの程度マッチしているかによって測定される．オイラー法は一次精度なので，誤差は一次高く$h^2$というベキ乗から始まる． ~NDSolveの”ExplicitRungeKutta”メソッドより

ともあります。

Butcher tableによるルンゲ=クッタ法の記述

行列形式で与えられるルンゲ=クッタ法での係数$a_{i,j}, b_{i},c_{i}$は何なのか？
具体的に記述してみましょう。
オイラー法(1段1次)はもっとも単純で、係数は
$
\begin{align}
a_{1,1}&=0　\\
b_{1}&=1 \\
c_{1}&=0
\end{align}
$
です。これを一般的な表記法の式に当てはめれば、
$
\begin{align}
g_1&=x_n+h a_{1,1}k_1 \\
k_1&=f(t_n+c_1h, g_1) \\
x_{n+1}&=x_n+h b_1k_1
\end{align}
$
となります。

中点法は、
$
\begin{align}
a_{1,1}&=0 \\
a_{1,2}&=0 \\
a_{2,1}&=1/2 \\
a_{2,2}&=0 \\
b_{1}&=0 \\
b_{2}&=1 \\
c_{1}&=0 \\
c_{2}&=1/2
\end{align}
$
という組で与えられます。

この係数行列の組はまとめてButcher tableと呼ばれる表記をするのが便利です。

これは、$a_{i,j}, b_{i},c_{i}$を

としてまとめて書く表記法です。

再び、オイラー法はButcher tableで書くと

とまとめて書くことができます。
中点法は

RK4は

です。

高次のルンゲ=クッタ法(10,12,14次)

4次、5次…とずっとあるわけです。
こんなページがありました[3]。
High-Order Explicit Runge-Kutta Methods
この上のページには

17段10次(8次が埋め込まれてる)
25段12次(10次が埋め込まれてる)
35段14次(12次が埋め込まれてる)

といったButcher tableにおける係数の値が書かれています。埋め込まれてる、の意味は次の節で説明します。
ただし、上のページのbutcher tableは

となっているので注意が必要です。

埋め込まれた陽的ルンゲ=クッタ法

“埋め込まれた“という表現が出てきたのでその説明を行いましょう。
日本語では『埋め込まれた陽的ルンゲ=クッタ法』、英語では『embedded explicit runge-kutta method』と呼ばれるものがあります。
これは、p段q次陽的ルンゲ=クッタ法を作ったら、別の次数の陽的ルンゲ=クッタ法も、係数行列$a_{i,j}, c_{i}$を使って作れるじゃありませんか!
というものです。

Butcher tableは、この場合extended Butcher tableと呼ばれ、こういう形式で書かれます。

この埋め込まれたルンゲ=クッタ法のいいところは、

計算誤差の評価ができる
刻み幅を自動的に制御できる、適応刻み幅制御。(応用として。)

という点です。ルンゲ=クッタ法によって得られた解が真の解とどのくらい違っているのか？が評価できるんです。

例えば、4次のルンゲ=クッタ法を使って得られた解$x^{(4)}(t)$と5次のルンゲ=クッタ法を使って得られた解$x^{(5)}(t)$があったとします。
もしも、$x^{(4)}(t)$と解$x^{(5)}(t)$の解の差を調べ、その差が無かったらその数値計算解は真の解に限りなく近い、と判断することができ、差が大きかったらその解は真の解から離れていて、数値計算の精度が足らない、と判断することができます。どちらも1つだけの解では出来なかったことです。

精度が足らない場合、刻み幅を小さくすれば精度が上がります。また逆に、精度が十分に足りている場合、刻み幅を大きくし、計算時間を減らすことができます。
これが適応刻み幅制御なのです。

違った次数のルンゲ=クッタ法を、まるで別々に計算してもいいのですが、そうすると計算時間が単純に考えておおよそ2倍になります。
解を評価するために2倍の計算時間が必要というのは良くない計算効率です。
そこで考えられたのが埋め込まれたルンゲ=クッタ法なのです。

具体例を挙げましょう。
一番簡単な埋め込まれたルンゲ=クッタ法は、ホイン法と呼ばれています。

1行目は2次のオーダーを持ち、2行目は1次のオーダーを持ちます。

また、4次と5次を持つ埋め込まれたルンゲ=クッタ法は、ルンゲ=クッタ=フェールベルグ(Runge-Kutta-Fehlberg)法と呼ばれています。
その埋め込まれたルンゲ=クッタ法は、

と書かれます。1行目は5次のオーダー、2行目は4次のオーダーを持ちます。

ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御(適応刻み幅制御)

さて、次数の違う2つのルンゲ=クッタ法を用いて、適応刻み幅制御を行いたいと考えます。
刻み幅を制御するにあたって、適当に精度良かったから2倍にしてもまだ大丈夫だろ、とか差が大きすぎるから刻み幅半分にしよう、ということをやってはいけません。
適当にやったら計算時間が余計にかかり、精度が良くない変な結果が得られます。

[5~9]によれば、ルンゲ=クッタ=フェールベルグ法において区間$i$での最適な刻み幅$h’$は区間$i$の誤差評価の結果を使って、
$
\displaystyle h’=\delta h=\left(\frac{\varepsilon h}{2|x^{(5)}_{i}-x^{(4)}_{i}|}\right)^{1/4} h
$
と予想できます。ここで$\varepsilon$はエラーを制御する精度の目安で、おおよそ計算時に要求する相対誤差です。もちろん、この$h’$は区間$i$の最適な幅ですが、関数に劇的な変化は無いだろうとして、次の区間の計算の一番初めに用いる刻み幅を推定するのです。
なので、$i+1$番目の計算区間では、計算するときはこの$h’$の値を使えばいいんです。
(ちなみに、m次ルンゲ=クッタ法の場合では
$
\displaystyle h’=\delta h=\left(\frac{\varepsilon h}{2|x^{(m+1)}_{i}-x^{(m)}_{i}|}\right)^{1/m} h
$
と推測されます。)

詳しい理由は分かりませんが、5次オーダーではなく、4次です。5次のオーダーは誤差評価のためだけに用いられているようです。
ルンゲ=クッタ=フェールベルグ法の計算スキームは[7]に詳しく書かれています。
日本語訳して、その計算スキームを書けば下のようになります。

ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御(プログラム)

講義のレポート等の宿題で使うのは僕の意向と異なるので使用はご控えください。
研究目的、趣味、確かめの場合はミスがあるかもしれないことを念頭に置いたうえならば使用と改変をご自由にしてください。
このプログラムを使用して生じた責任は取りません。

fortran90によるプログラムです。ほぼ上の説明をそのままプログラミングしたものです。

実数、1階微分方程式の場合

微分方程式
$
\displaystyle \frac{d}{dx}y=y\cos x,~~y(x=0)=1
$
を倍精度実数、刻み幅制御で$x=10$まで解く事を考えます。
解析解は
$
\displaystyle y(x)=exp(sin(x))
$
です。コードは以下の通りです。

▼ここクリックでこの場に展開

program main
implicit none
integer::i,N,info
double precision::x,h,tol,xbound
double precision,allocatable::y(:)
external::grk

N=1 ! Number of differential equations
allocate(y(1:N))

x=0d0; xbound=10d0
y(1)=1d0 !initial condition

h=xbound-x
i=0; info=0; tol=1d-8;
do while(info.le.0)
call drkf45(grk,x,h,N,y,xbound,info,tol)
write(10,'(3e25.10e3)')x,y(1),h
i=i+1
enddo
write(6,*)"Number of referenced times -->",i

stop
end program main

subroutine grk(N,x,y,f)
implicit none
integer,intent(in)::N
double precision,intent(in)::x,y(1:N)
double precision,intent(out)::f(1:N)

! Solve
! d y(1) / dt = y(1) * cos(t)

f(1)=y(1)*cos(x)

return
end subroutine grk

!===============================

subroutine drkf45(grk,x,h,N,y,xbound,info,tol)
! if h < hmin, propagate forcibly with warning.
!
!-----------------
!info = -9 (maybe path the discontinue points)
! = 0 (Running now)
! = 1 (x reach xbound)
!-----------------
!
implicit none
integer,intent(in)::N
double precision,intent(in)::xbound,tol
double precision,intent(inout)::x,h,y(1:N)
integer,intent(inout)::info

integer::i,j,FLAG,key
double precision::R,delta,tx,Sy,err
double precision,allocatable::ty(:),K(:,:),tf(:)
double precision,parameter::hmin=1d-14,hmax=0.5d0
integer,parameter::s=6
double precision::a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s)
external::grk

c(1:6)=(/0d0, 0.25d0, 0.375d0,&
0.9230769230769230769230769230769230769231d0, 1d0, 0.5d0/)
a(1:6,1:6)=0d0
a(1,1:6)=(/0d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(2,1:6)=(/0.25d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(3,1:6)=(/0.09375d0, 0.28125d0, 0d0, 0d0, 0d0, 0d0/)
a(4,1:6)=(/0.8793809740555302685480200273099681383705d0, &
-3.277196176604460628129267182521620391443d0, &
3.320892125625853436504324078288575329995d0, 0d0, 0d0, 0d0/)
a(5,1:6)=(/2.032407407407407407407407407407407407407d0,-8d0, &
7.173489278752436647173489278752436647173d0, &
-0.2058966861598440545808966861598440545809d0, 0d0, 0d0/)
a(6,1:6)=(/-0.2962962962962962962962962962962962962963d0,2d0, &
-1.381676413255360623781676413255360623782d0, &
0.4529727095516569200779727095516569200780d0,-0.275d0,0d0/)
b2(1:6)=(/0.1185185185185185185185185185185185185185d0, 0.d0,&
0.5189863547758284600389863547758284600390d0, &
0.5061314903420166578061314903420166578061d0, &
-0.18d0, 0.03636363636363636363636363636363636363636d0/)
b1(1:6)=(/0.1157407407407407407407407407407407407407d0, 0d0,&
0.5489278752436647173489278752436647173489d0, &
0.5353313840155945419103313840155945419103d0, -0.2d0, 0d0/)
Rc(1:6)=(/0.002777777777777777777777777777777777777778d0,0d0, &
-0.02994152046783625730994152046783625730994d0, &
-0.02919989367357788410419989367357788410420d0, 0.02d0, &
0.03636363636363636363636363636363636363636d0/)

key=0
allocate(ty(1:N),tf(1:N),K(1:s,1:N))
ty=0d0; tf=0d0; K=0d0

if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
endif

if(h.ge.abs(xbound-x))h=xbound-x

FLAG=1
if(abs(x-xbound).le.hmin)then
info=1
FLAG=0
endif

do while(FLAG.eq.1)
tx=x
do j=1,s
tx=x+c(j)*h
ty(1:N)=y(1:N)
do i=1,j-1
ty(1:N)=ty(1:N)+K(i,1:N)*a(j,i)
enddo
call grk(N,tx,ty,tf)
K(j,1:N)=h*tf(1:N)
enddo

!step 4
R=0d0
do i=1,N
R=R+(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i) &
+Rc(5)*K(5,i)+Rc(6)*K(6,i))**2d0
enddo
R=abs(dsqrt(R)/h/dble(N))

Sy=0d0
do i=1,N
Sy=Sy+(y(i)*y(i))
enddo
Sy=dsqrt(Sy)
if(Sy.ge.1d0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
y(1:N)=y(1:N)+b1(1)*K(1,1:N)+b1(3)*K(3,1:N) &
+b1(4)*K(4,1:N)+b1(5)*K(5,1:N)
FLAG=0
endif

!step 6
! Avoid zero deviding.
if(R.ge.1d-20)then
delta=(err/(2d0*R))**0.25d0
else
delta=4d0
endif

!step 7
if(delta.le.0.1d0)then
!function changes dramatically.
h=0.1d0*h
elseif(delta.ge.4d0)then
!function changes loosely.
h=4d0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1d0,h)*hmin
key=1
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if

if(h.le.0d0.and.xbound-x.ge.0d0)then
info=1
FLAG=0
elseif(h.ge.0d0.and.xbound-x.le.0d0)then
info=1
FLAG=0
endif
enddo

if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')&
"Strange point between ",x-h," and ",x
info=-9
endif

deallocate(ty,tf,K)
return
end subroutine drkf45

実数、2階微分方程式の場合

微分方程式
$
\displaystyle \frac{d^2}{dx^2}y=-\frac{1}{4}y,~~y(x=0)=1,~ y'(x=0)=0
$
を倍精度実数、刻み幅制御で$x=20$まで解く事を考えます。コードは以下の通りです。

▼ここクリックでこの場に展開

program main
implicit none
integer::i,N,info
double precision::x,h,tol,xbound
double precision,allocatable::y(:)
external::grk

N=2
allocate(y(1:N))

x=0d0; xbound=20d0
!initial condition
y(1)=1d0 ! y (0) = 1d0
y(2)=0d0 ! y'(0) = 0d0

h=xbound-x
i=0; info=0; tol=1d-8;
do while(info.le.0)
call drkf45(grk,x,h,N,y,xbound,info,tol)
write(10,'(4e25.10e3)')x,y(1),y(2),h
i=i+1
enddo
write(6,*)"Number of referenced times -->",i

stop
end program main

subroutine grk(N,x,y,f)
implicit none
integer,intent(in)::N
double precision,intent(in)::x,y(1:N)
double precision,intent(out)::f(1:N)

! Solve
! d^2 y(1) / dt^2 = - 0.5 * y(1)

f(1)=y(2)
f(2)=-0.25d0*y(1)

return
end subroutine grk

!===============================

subroutine drkf45(grk,x,h,N,y,xbound,info,tol)
! if h < hmin, propagate forcibly with warning.
!
!-----------------
!info = -9 (maybe path the discontinue points)
! = 0 (Running now)
! = 1 (x reach xbound)
!-----------------
!
implicit none
integer,intent(in)::N
double precision,intent(in)::xbound,tol
double precision,intent(inout)::x,h,y(1:N)
integer,intent(inout)::info

integer::i,j,FLAG,key
double precision::R,delta,tx,Sy,err
double precision,allocatable::ty(:),K(:,:),tf(:)
double precision,parameter::hmin=1d-14,hmax=0.5d0
integer,parameter::s=6
double precision::a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s)
external::grk

c(1:6)=(/0d0, 0.25d0, 0.375d0,&
0.9230769230769230769230769230769230769231d0, 1d0, 0.5d0/)
a(1:6,1:6)=0d0
a(1,1:6)=(/0d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(2,1:6)=(/0.25d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(3,1:6)=(/0.09375d0, 0.28125d0, 0d0, 0d0, 0d0, 0d0/)
a(4,1:6)=(/0.8793809740555302685480200273099681383705d0, &
-3.277196176604460628129267182521620391443d0, &
3.320892125625853436504324078288575329995d0, 0d0, 0d0, 0d0/)
a(5,1:6)=(/2.032407407407407407407407407407407407407d0,-8d0, &
7.173489278752436647173489278752436647173d0, &
-0.2058966861598440545808966861598440545809d0, 0d0, 0d0/)
a(6,1:6)=(/-0.2962962962962962962962962962962962962963d0,2d0, &
-1.381676413255360623781676413255360623782d0, &
0.4529727095516569200779727095516569200780d0,-0.275d0,0d0/)
b2(1:6)=(/0.1185185185185185185185185185185185185185d0, 0.d0,&
0.5189863547758284600389863547758284600390d0, &
0.5061314903420166578061314903420166578061d0, &
-0.18d0, 0.03636363636363636363636363636363636363636d0/)
b1(1:6)=(/0.1157407407407407407407407407407407407407d0, 0d0,&
0.5489278752436647173489278752436647173489d0, &
0.5353313840155945419103313840155945419103d0, -0.2d0, 0d0/)
Rc(1:6)=(/0.002777777777777777777777777777777777777778d0,0d0, &
-0.02994152046783625730994152046783625730994d0, &
-0.02919989367357788410419989367357788410420d0, 0.02d0, &
0.03636363636363636363636363636363636363636d0/)

key=0
allocate(ty(1:N),tf(1:N),K(1:s,1:N))
ty=0d0; tf=0d0; K=0d0

if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
endif

if(h.ge.abs(xbound-x))h=xbound-x

FLAG=1
if(abs(x-xbound).le.hmin)then
info=1
FLAG=0
endif

do while(FLAG.eq.1)
tx=x
do j=1,s
tx=x+c(j)*h
ty(1:N)=y(1:N)
do i=1,j-1
ty(1:N)=ty(1:N)+K(i,1:N)*a(j,i)
enddo
call grk(N,tx,ty,tf)
K(j,1:N)=h*tf(1:N)
enddo

!step 4
R=0d0
do i=1,N
R=R+(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i) &
+Rc(5)*K(5,i)+Rc(6)*K(6,i))**2d0
enddo
R=abs(dsqrt(R)/h/dble(N))

Sy=0d0
do i=1,N
Sy=Sy+(y(i)*y(i))
enddo
Sy=dsqrt(Sy)
if(Sy.ge.1d0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
y(1:N)=y(1:N)+b1(1)*K(1,1:N)+b1(3)*K(3,1:N) &
+b1(4)*K(4,1:N)+b1(5)*K(5,1:N)
FLAG=0
endif

!step 6
! Avoid zero deviding.
if(R.ge.1d-20)then
delta=(err/(2d0*R))**0.25d0
else
delta=4d0
endif

!step 7
if(delta.le.0.1d0)then
!function changes dramatically.
h=0.1d0*h
elseif(delta.ge.4d0)then
!function changes loosely.
h=4d0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1d0,h)*hmin
key=1
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if

if(h.le.0d0.and.xbound-x.ge.0d0)then
info=1
FLAG=0
elseif(h.ge.0d0.and.xbound-x.le.0d0)then
info=1
FLAG=0
endif
enddo

if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')&
"Strange point between ",x-h," and ",x
info=-9
endif

deallocate(ty,tf,K)
return
end subroutine drkf45

※ここで使われているサブルーチンdrkf45は実数一階微分のプログラム内にあるルーチンと一字一句同一です。

複素数、1階微分方程式の場合

微分方程式
$
\displaystyle \frac{d}{dx}y=y\cos x,~~ y(x=0)=1+i\frac{1}{2}
$
を倍精度実数、刻み幅制御で$x=10$まで解く事を考えます。コードは以下の通りです。

▼ここクリックでこの場に展開

program main
implicit none
integer::i,N,info
double precision::x,h,tol,xbound
complex(kind(0d0)),allocatable::y(:)
external::grk

N=1
allocate(y(1:N))

x=0d0; xbound=10d0
y(1)=dcmplx(1d0,0.5d0) !initial condition

h=xbound-x
i=0; info=0; tol=1d-8;
do while(info.le.0)
call crkf45(grk,x,h,N,y,xbound,info,tol)
write(10,'(4e25.10e3)')x,y(1),h
i=i+1
enddo
write(6,*)"Number of referenced times -->",i

stop
end program main

subroutine grk(N,x,y,f)
implicit none
integer,intent(in)::N
double precision,intent(in)::x
complex(kind(0d0)),intent(in)::y(1:N)
complex(kind(0d0)),intent(out)::f(1:N)

f(1)=y(1)*cos(x)

return
end subroutine grk

!===============================

subroutine crkf45(grk,x,h,N,y,xbound,info,tol)
! if h < hmin, propagate forcibly with warning.
!
!-----------------
!info = -9 (maybe path the discontinue points)
! = 0 (Running now)
! = 1 (x reach xbound)
!-----------------
!
implicit none
integer,intent(in)::N
double precision,intent(in)::xbound,tol
double precision,intent(inout)::x,h
complex(kind(0d0)),intent(inout)::y(1:N)
integer,intent(inout)::info
external::grk

integer::i,j,FLAG,key
double precision::R,delta,tx,Sy,err
complex(kind(0d0)),allocatable::ty(:),tf(:),K(:,:)
double precision,parameter::hmin=1d-14,hmax=0.5d0
integer,parameter::s=6
double precision::a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s)

c(1:6)=(/0d0, 0.25d0, 0.375d0,&
0.9230769230769230769230769230769230769231d0, 1d0, 0.5d0/)
a(1:6,1:6)=0d0
a(1,1:6)=(/0d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(2,1:6)=(/0.25d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(3,1:6)=(/0.09375d0, 0.28125d0, 0d0, 0d0, 0d0, 0d0/)
a(4,1:6)=(/0.8793809740555302685480200273099681383705d0, &
-3.277196176604460628129267182521620391443d0, &
3.320892125625853436504324078288575329995d0, 0d0, 0d0, 0d0/)
a(5,1:6)=(/2.032407407407407407407407407407407407407d0,-8d0, &
7.173489278752436647173489278752436647173d0, &
-0.2058966861598440545808966861598440545809d0, 0d0, 0d0/)
a(6,1:6)=(/-0.2962962962962962962962962962962962962963d0,2d0, &
-1.381676413255360623781676413255360623782d0, &
0.4529727095516569200779727095516569200780d0,-0.275d0,0d0/)
b2(1:6)=(/0.1185185185185185185185185185185185185185d0, 0.d0,&
0.5189863547758284600389863547758284600390d0, &
0.5061314903420166578061314903420166578061d0, &
-0.18d0, 0.03636363636363636363636363636363636363636d0/)
b1(1:6)=(/0.1157407407407407407407407407407407407407d0, 0d0,&
0.5489278752436647173489278752436647173489d0, &
0.5353313840155945419103313840155945419103d0, -0.2d0, 0d0/)
Rc(1:6)=(/0.002777777777777777777777777777777777777778d0,0d0, &
-0.02994152046783625730994152046783625730994d0, &
-0.02919989367357788410419989367357788410420d0, 0.02d0, &
0.03636363636363636363636363636363636363636d0/)

key=0
allocate(ty(1:N),tf(1:N),K(1:s,1:N))
ty=dcmplx(0d0,0d0); K=dcmplx(0d0,0d0)
tf=dcmplx(0d0,0d0)

if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
endif

if(h.ge.abs(xbound-x))h=xbound-x

FLAG=1
if(abs(x-xbound).le.hmin)then
info=1
FLAG=0
endif

do while(FLAG.eq.1)
tx=x
do j=1,s
tx=x+c(j)*h
ty(1:N)=y(1:N)
do i=1,j-1
ty(1:N)=ty(1:N)+K(i,1:N)*a(j,i)
enddo
!do i=1,N
! K(j,i)=h*grk(N,tx,ty,i)
!enddo
call grk(N,tx,ty,tf)
K(j,1:N)=h*tf(1:N)
enddo

!step 4
R=0d0
do i=1,N
R=R+abs(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i) &
+Rc(5)*K(5,i)+Rc(6)*K(6,i))**2d0
enddo
R=abs(dsqrt(R)/h/dble(N))

Sy=0d0
do i=1,N
Sy=Sy+abs(y(i)*y(i))
enddo
Sy=dsqrt(Sy)
if(Sy.ge.1d0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
y(1:N)=y(1:N)+b1(1)*K(1,1:N)+b1(3)*K(3,1:N) &
+b1(4)*K(4,1:N)+b1(5)*K(5,1:N)
FLAG=0
endif

!step 6
! Avoid zero deviding.
if(R.ge.1d-20)then
delta=(err/(2d0*R))**0.25d0
else
delta=4d0
endif

!step 7
if(delta.le.0.1d0)then
!function changes dramatically.
h=0.1d0*h
elseif(delta.ge.4d0)then
!function changes loosely.
h=4d0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1d0,h)*hmin
key=1
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if

if(h.le.0d0.and.xbound-x.ge.0d0)then
info=1
FLAG=0
elseif(h.ge.0d0.and.xbound-x.le.0d0)then
info=1
FLAG=0
endif
enddo

if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')&
"Strange point between ",x-h," and ",x
info=-9
endif

deallocate(ty,tf,K)
return
end subroutine crkf45

複素数、2階微分方程式

微分方程式
$
\displaystyle \frac{d^2}{dx^2}y=-\frac{1}{4}y,~~ y(x=0)=1+i\frac{1}{2},~y'(x=0)=0+i0
$
を倍精度実数、刻み幅制御で$x=20$まで解く事を考えます。コードは以下の通りです。

▼ここクリックでこの場に展開

program main
implicit none
integer::i,N,info
double precision::x,h,tol,xbound
complex(kind(0d0)),allocatable::y(:)
external::grk

N=2
allocate(y(1:N))

x=0d0; xbound=20d0
y(1)=dcmplx(1d0,0.5d0) !initial condition
y(2)=dcmplx(0d0,0d0) !initial condition

h=xbound-x
i=0; info=0; tol=1d-8;
do while(info.le.0)
call crkf45(grk,x,h,N,y,xbound,info,tol)
write(10,'(6e25.10e3)')x,dble(y(1)),dimag(y(1)) &
,dble(y(2)),dimag(y(2)),h
i=i+1
enddo
write(6,*)"Number of referenced times -->",i

stop
end program main

subroutine grk(N,x,y,f)
implicit none
integer,intent(in)::N
double precision,intent(in)::x
complex(kind(0d0)),intent(in)::y(1:N)
complex(kind(0d0)),intent(out)::f(1:N)

f(1)=y(2)
f(2)=-0.25d0*y(1)

return
end subroutine grk

!===============================

subroutine crkf45(grk,x,h,N,y,xbound,info,tol)
! if h < hmin, propagate forcibly with warning.
!
!-----------------
!info = -9 (maybe path the discontinue points)
! = 0 (Running now)
! = 1 (x reach xbound)
!-----------------
!
implicit none
integer,intent(in)::N
double precision,intent(in)::xbound,tol
double precision,intent(inout)::x,h
complex(kind(0d0)),intent(inout)::y(1:N)
integer,intent(inout)::info
external::grk

integer::i,j,FLAG,key
double precision::R,delta,tx,Sy,err
complex(kind(0d0)),allocatable::ty(:),tf(:),K(:,:)
double precision,parameter::hmin=1d-14,hmax=0.5d0
integer,parameter::s=6
double precision::a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s)

c(1:6)=(/0d0, 0.25d0, 0.375d0,&
0.9230769230769230769230769230769230769231d0, 1d0, 0.5d0/)
a(1:6,1:6)=0d0
a(1,1:6)=(/0d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(2,1:6)=(/0.25d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(3,1:6)=(/0.09375d0, 0.28125d0, 0d0, 0d0, 0d0, 0d0/)
a(4,1:6)=(/0.8793809740555302685480200273099681383705d0, &
-3.277196176604460628129267182521620391443d0, &
3.320892125625853436504324078288575329995d0, 0d0, 0d0, 0d0/)
a(5,1:6)=(/2.032407407407407407407407407407407407407d0,-8d0, &
7.173489278752436647173489278752436647173d0, &
-0.2058966861598440545808966861598440545809d0, 0d0, 0d0/)
a(6,1:6)=(/-0.2962962962962962962962962962962962962963d0,2d0, &
-1.381676413255360623781676413255360623782d0, &
0.4529727095516569200779727095516569200780d0,-0.275d0,0d0/)
b2(1:6)=(/0.1185185185185185185185185185185185185185d0, 0.d0,&
0.5189863547758284600389863547758284600390d0, &
0.5061314903420166578061314903420166578061d0, &
-0.18d0, 0.03636363636363636363636363636363636363636d0/)
b1(1:6)=(/0.1157407407407407407407407407407407407407d0, 0d0,&
0.5489278752436647173489278752436647173489d0, &
0.5353313840155945419103313840155945419103d0, -0.2d0, 0d0/)
Rc(1:6)=(/0.002777777777777777777777777777777777777778d0,0d0, &
-0.02994152046783625730994152046783625730994d0, &
-0.02919989367357788410419989367357788410420d0, 0.02d0, &
0.03636363636363636363636363636363636363636d0/)

key=0
allocate(ty(1:N),tf(1:N),K(1:s,1:N))
ty=dcmplx(0d0,0d0); K=dcmplx(0d0,0d0)
tf=dcmplx(0d0,0d0)

if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
endif

if(h.ge.abs(xbound-x))h=xbound-x

FLAG=1
if(abs(x-xbound).le.hmin)then
info=1
FLAG=0
endif

do while(FLAG.eq.1)
tx=x
do j=1,s
tx=x+c(j)*h
ty(1:N)=y(1:N)
do i=1,j-1
ty(1:N)=ty(1:N)+K(i,1:N)*a(j,i)
enddo
call grk(N,tx,ty,tf)
K(j,1:N)=h*tf(1:N)
enddo

!step 4
R=0d0
do i=1,N
R=R+abs(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i) &
+Rc(5)*K(5,i)+Rc(6)*K(6,i))**2d0
enddo
R=abs(dsqrt(R)/h/dble(N))

Sy=0d0
do i=1,N
Sy=Sy+abs(y(i)*y(i))
enddo
Sy=dsqrt(Sy)
if(Sy.ge.1d0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
y(1:N)=y(1:N)+b1(1)*K(1,1:N)+b1(3)*K(3,1:N) &
+b1(4)*K(4,1:N)+b1(5)*K(5,1:N)
FLAG=0
endif

!step 6
! Avoid zero deviding.
if(R.ge.1d-20)then
delta=(err/(2d0*R))**0.25d0
else
delta=4d0
endif

!step 7
if(delta.le.0.1d0)then
!function changes dramatically.
h=0.1d0*h
elseif(delta.ge.4d0)then
!function changes loosely.
h=4d0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1d0,h)*hmin
key=1
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if

if(h.le.0d0.and.xbound-x.ge.0d0)then
info=1
FLAG=0
elseif(h.ge.0d0.and.xbound-x.le.0d0)then
info=1
FLAG=0
endif
enddo

if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')&
"Strange point between ",x-h," and ",x
info=-9
endif

deallocate(ty,tf,K)
return
end subroutine crkf45

※ここで使われているサブルーチンcrkf45は複素数一階微分のプログラム内にあるルーチンと一字一句同一です。

等間隔の出力の場合は、以下の通りで実行できます。
サブルーチンはdrkf45は変わっていません。

▼ここクリックでこの場に展開

program main
implicit none
integer::i,N,info,Nx
double precision::h,tol
double precision::xa,xb,tx
double precision,allocatable::y(:),x(:)
external::grk

N=1
allocate(y(1:N))

Nx=101
allocate(x(1:Nx))
xa=0d0
xb=10d0
do i=1,Nx
x(i)=(i-1)*(xb-xa)/dble(Nx-1)+xa
enddo

!initial conditions
y(1)=1d0 ! x (0)

tol=1d-8
write(10,'(2e25.10e3)')x(1),y(1)
do i=2,Nx
info=-1
h=x(i)-x(i-1)
tx=x(i-1)
do while(info.le.0)
call drkf45(grk,tx,h,N,y,x(i),info,tol)
enddo
write(10,'(2e25.10e3)')x(i),y(1)
enddo

stop
end program main

subroutine grk(N,x,y,f)
implicit none
integer,intent(in)::N
double precision,intent(in)::x,y(1:N)
double precision,intent(out)::f(1:N)

! Solve
! d y(1) / dt = y(1) * cos(t)

f(1)=y(1)*cos(x)

return
end subroutine grk

!===============================

subroutine drkf45(grk,x,h,N,y,xbound,info,tol)
! if h < hmin, propagate forcibly with warning.
!
!-----------------
!info = -9 (maybe path the discontinue points)
! = 0 (Running now)
! = 1 (x reach xbound)
!-----------------
!
implicit none
integer,intent(in)::N
double precision,intent(in)::xbound,tol
double precision,intent(inout)::x,h,y(1:N)
integer,intent(inout)::info

integer::i,j,FLAG,key
double precision::R,delta,tx,Sy,err
double precision,allocatable::ty(:),K(:,:),tf(:)
double precision,parameter::hmin=1d-14,hmax=0.5d0
integer,parameter::s=6
double precision::a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s)
external::grk

c(1:6)=(/0d0, 0.25d0, 0.375d0,&
0.9230769230769230769230769230769230769231d0, 1d0, 0.5d0/)
a(1:6,1:6)=0d0
a(1,1:6)=(/0d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(2,1:6)=(/0.25d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(3,1:6)=(/0.09375d0, 0.28125d0, 0d0, 0d0, 0d0, 0d0/)
a(4,1:6)=(/0.8793809740555302685480200273099681383705d0, &
-3.277196176604460628129267182521620391443d0, &
3.320892125625853436504324078288575329995d0, 0d0, 0d0, 0d0/)
a(5,1:6)=(/2.032407407407407407407407407407407407407d0,-8d0, &
7.173489278752436647173489278752436647173d0, &
-0.2058966861598440545808966861598440545809d0, 0d0, 0d0/)
a(6,1:6)=(/-0.2962962962962962962962962962962962962963d0,2d0, &
-1.381676413255360623781676413255360623782d0, &
0.4529727095516569200779727095516569200780d0,-0.275d0,0d0/)
b2(1:6)=(/0.1185185185185185185185185185185185185185d0, 0.d0,&
0.5189863547758284600389863547758284600390d0, &
0.5061314903420166578061314903420166578061d0, &
-0.18d0, 0.03636363636363636363636363636363636363636d0/)
b1(1:6)=(/0.1157407407407407407407407407407407407407d0, 0d0,&
0.5489278752436647173489278752436647173489d0, &
0.5353313840155945419103313840155945419103d0, -0.2d0, 0d0/)
Rc(1:6)=(/0.002777777777777777777777777777777777777778d0,0d0, &
-0.02994152046783625730994152046783625730994d0, &
-0.02919989367357788410419989367357788410420d0, 0.02d0, &
0.03636363636363636363636363636363636363636d0/)

key=0
allocate(ty(1:N),tf(1:N),K(1:s,1:N))
ty=0d0; tf=0d0; K=0d0

if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
endif

if(h.ge.abs(xbound-x))h=xbound-x

FLAG=1
if(abs(x-xbound).le.hmin)then
info=1
FLAG=0
endif

do while(FLAG.eq.1)
tx=x
do j=1,s
tx=x+c(j)*h
ty(1:N)=y(1:N)
do i=1,j-1
ty(1:N)=ty(1:N)+K(i,1:N)*a(j,i)
enddo
call grk(N,tx,ty,tf)
K(j,1:N)=h*tf(1:N)
enddo

!step 4
R=0d0
do i=1,N
R=R+(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i) &
+Rc(5)*K(5,i)+Rc(6)*K(6,i))**2d0
enddo
R=abs(dsqrt(R)/h/dble(N))

Sy=0d0
do i=1,N
Sy=Sy+(y(i)*y(i))
enddo
Sy=dsqrt(Sy)
if(Sy.ge.1d0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
y(1:N)=y(1:N)+b1(1)*K(1,1:N)+b1(3)*K(3,1:N) &
+b1(4)*K(4,1:N)+b1(5)*K(5,1:N)
FLAG=0
endif

!step 6
! Avoid zero deviding.
if(R.ge.1d-20)then
delta=(err/(2d0*R))**0.25d0
else
delta=4d0
endif

!step 7
if(delta.le.0.1d0)then
!function changes dramatically.
h=0.1d0*h
elseif(delta.ge.4d0)then
!function changes loosely.
h=4d0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1d0,h)*hmin
key=1
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if

if(h.le.0d0.and.xbound-x.ge.0d0)then
info=1
FLAG=0
elseif(h.ge.0d0.and.xbound-x.le.0d0)then
info=1
FLAG=0
endif
enddo

if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')&
"Strange point between ",x-h," and ",x
info=-9
endif

deallocate(ty,tf,K)
return
end subroutine drkf45

[adsense2]

不連続な点を含む場合

不連続な点を含む場合、境界条件を指定しないと解くことはできません。

さて、ここで微分方程式
$
\begin{eqnarray}
\frac{dy}{dx}=
\left\{
\begin{aligned}
0\;\;(x\le 0)\\
1\;\;(0\lt x)
\end{aligned}
\right.
\end{eqnarray}
$
を初期条件$y(-0.095)=0$の下で考えます(意図的に境界条件は考えません)。
これを4次ルンゲ=クッタと適応刻み幅ルンゲ=クッタで解いてみましょう。
プログラム上ではそのまま解くことが出来ます。
実際に解かせてみますと、

となり、2つの結果(赤：4次ルンゲ=クッタ、緑：刻み幅制御ルンゲ=クッタ)は異なってしまいます。これは1階微分の不連続性のため発生します。
不連続点$x=0$で関数$y(x)$に境界条件を指定しない限り、どちらも正しい解なのです。

さて、なぜこんなことが発生するのでしょうか？以下のように問題を表すことにします。

不連続な点を含む1階の微分方程式を考えます。
ここで不連続、という意味は関数$y(x)$の一階微分が、点$x’$で
$
\displaystyle \left. \frac{dy}{dx}\right|^{x=x’+0}_{x=x’-0}=a,\;\;\;(a\ne 0)
$
であるような点を指しているとします。

上記の例題を考えてみましょう。
上記の例題では$a=1$です。微分方程式を解析的に解いてみますと、
$
\begin{eqnarray}
y(x)=
\left\{
\begin{aligned}
C_0\;\;(x\le 0)\\
x+C_1\;\;(0\lt x)
\end{aligned}
\right.
\end{eqnarray}
$
となります。ここで$C_0, C_1$は定数です。
$C_0, C_1$は$y(x)$が解きたい問題の境界条件によって決まります。

例えば、$y(x)$は全領域に対して繋がっている、という条件を課しましょう。この場合、不連続点$x’$で
$
\displaystyle \left. y(x)\right|^{x=x’+0}_{x=x’-0}=0
$
という境界条件を満たさなければなりません。この条件を課すと、$C_1=C_0$となり、初めて関数$y(x)$を一意に決めることが出来ます。

1階微分方程式を解く場合、適応刻み幅制御では関数$y(x)$は計算領域内で繋がっている事が課されています。しかし、4次ルンゲ=クッタではその条件は課されません。$C_1$の値は初期条件に依存し、一意に関数が決まりません。
どちらが悪いという話ではありません。

通常は適応刻み幅でも、4次ルンゲ=クッタでも$y(x)$にどこか連続ではない変な点がある場合、その点で区間を別々に分けて解きます。その後、境界条件に従って値を調節して全体の関数を構成します。

ベンチマーク用

微分方程式の解法がどれくらい正しそうかのベンチマーク問題として振り子(角度が大きい時)を考えましょう(振り子の詳しい解説はこちら)。
以下の$\omega=1$としたときの運動方程式
$
\displaystyle \frac{d^2\theta}{dt^2}=-\sin\theta
$
を
$
t=0,\; \theta=0; \frac{d\theta}{dt}=1.9\cdots (k=0.95)
$
の初期条件の下解いた場合、
1周期$T$は
$
T=4K(0.95)=10.360044923498004876778\cdots
$
となります。
この値はwolfram alphaから求めました。
4EllipticK[0.950.95] wolfram alpha

刻み幅制御を行い、45000周期目の値を考えます。45000周期目は時刻
$
T_{45000}=466202.0215574102\cdots
$
です。刻み幅制御による精度を$10^{-12}$に設定し、数値計算を行わせます。

すると実行結果として”fort.10″に

0.4662020113E+006 -0.2103356901E-001 0.1899883579E+001 0.4224109363E-002
0.4662020155E+006 -0.1300808922E-001 0.1899955473E+001 0.4223843994E-002
0.4662020198E+006 -0.4982881533E-002 0.1899993468E+001 0.4223658015E-002
0.4662020240E+006 0.3042061693E-002 0.1899997567E+001 0.4223520921E-002

というデータが出力されます。
1列目が時刻$t$、2列目が$\theta(t)$,3列目が$\frac{d\theta(t)}{dt}$,4列目が刻み幅$h$です。
1回のステップでの要求精度12桁に対し、最終的な結果は8桁まで正しい値を出しています。
この時、計算回数は77,852,488回$~10^{8}$回行われているので、最終的な結果として4桁は少なくとも正しいと考えられます。
4桁合っていればいい状況で8桁もあっているのは、被積分関数が変な振る舞いをしないからでしょう。

また、60000周期で確認してみると($T_{60000}=621602.695409880292\cdots$)

0.6216026888E+006 -0.1531918479E-001 0.1899938246E+001 0.4223959920E-002
0.6216026930E+006 -0.7293808996E-002 0.1899986003E+001 0.4223717084E-002
0.6216026973E+006 0.7312355417E-003 0.1899999862E+001 0.4223582630E-002
0.6216027015E+006 0.8756011575E-002 0.1899979827E+001 0.4223462029E-002

です。
1回のステップでの要求精度12桁に対し、最終的な結果は8桁まで正しい値を出しています。
この時、計算回数は103,803,513回$~10^{8}$回行われているので、最終的な結果として4桁は少なくとも正しいと考えられます。
･･･まだまだ大丈夫そうですね。

少し特殊な初期条件($k=1$)でやってみましょう。
この$k=1$は、振り子の質点がちょうど真上に来て静止する非常に不安定な状態です。
何秒間静止していられるか試してみましょう。刻み幅の制御等は上記条件と同じです。
横軸に時間$t$、縦軸に$\theta(t)$を取った時のグラフです。

すぐに破綻しました。正しい値は$t=20$位までですかね。これは、機械的な誤差があることによって不安定な平衡点からほんの少しだけ外れてしまったからです。だからカオスとかを考えるときとかは注意しなければなりません。

ルンゲ=クッタ=ドルマンド=プリンス法

フェールベルグ法は初期の頃に考えられた埋め込まれた方法です。
研究が進み、実用上では別の係数の組み合わせが良いことが分かってきました。
その一つが1980年に発見されたドルマンド=プリンス(Dormand-Prince)法です。

これは、7段4(5)次の方法です。
フェールベルグ法は6段4(5)次の方法ですので、次数は同じです。

良く調べていませんが、この違いは、4次の結果を基準にして求めたブッチャー係数(フェールベルグ法)か、5次の結果を基準に求めたブッチャー係数(ドルマンド=プリンス法)か？の違いのようです。

単純に考えて、同じ次数なのにドルマンド=プリンス法の方が段数が増えていて効率が悪いです。
しかし、本来は7段なのですが、7段目に呼び出した結果を取って置けば、次のステップの1段目に同じ値が使えるように設計されているので、プログラム上は6段と(ほぼ)同じ関数の呼び出し回数になります。

プログラムはこんな感じになるでしょう。

適当な刻み幅で出力

▼ここクリックでこの場に展開

program main
implicit none
integer::i,N,info
double precision::x,h,tol,xbound
double precision,allocatable::y(:),work0(:)
double precision,external::grk

x=0d0; xbound=10d0

N=1
allocate(y(1:N),work0(1:N))

!initial conditions
y(1)=1d0 ! x (0)

h=xbound-x
info=-1
tol=1d-8
i=0
write(20,'(3e30.16e3)')x,y(1),h
do while(info.le.0)
call dDP45(grk,x,h,N,y,xbound,info,tol,work0)
write(20,'(3e30.16e3)')x,y(1),h
i=i+1
enddo
write(6,*)"the Number, routine called-->",i

stop
end program main

function grk(N,x,y,s)
implicit none
integer,intent(in)::N,s
double precision,intent(in)::x
double precision,intent(in)::y(1:N)
double precision::grk

grk=0d0
if(s.eq.1)then
grk=y(1)*cos(x)
else
write(6,*)"***Error grk"; stop
endif

return
end function grk

!===============================

subroutine dDP45(grk,x,h,N,y,xbound,info,tol,work0)
implicit none
integer,intent(in)::N
double precision,intent(in)::xbound,tol
double precision,intent(inout)::x,h,y(1:N),work0(1:N)
integer,intent(inout)::info
double precision,external::grk

! Runge-Kutta-Dormand-Prince method
!
! if h < hmin, propagate forcibly with warning.
!
!-----------------
!work0(1:N) : Just prepare, used and substituted by
! the coefficients K(7,1:N) to reduce computation costs
!info = -1 (Used for first loop;
! you should set info=0 when you try to use this routine)
! = 0 (Running now)
! = 1 (x reach xbound)
! = -2 (Path the discontinue points)
!-----------------
!

integer::i,j,FLAG,key
double precision::R,delta,tx,Sy,err
double precision,allocatable::tmp(:),K(:,:)
double precision,parameter::hmin=1d-12,hmax=0.5d0
integer,parameter::s=7
double precision::a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s),work(1:N)

c(1:7)=(/0d0, 0.2d0, 0.3d0,0.8d0,&
0.888888888888888888888888888888888888888889d0,&
1.d0, 1.d0/)
a(1:7,1:7)=0d0
a(1,1:7)=(/0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(2,1:7)=(/0.2d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(3,1:7)=(/0.075d0, 0.225d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(4,1:7)=(/0.977777777777777777777777777777777777778d0, &
-3.733333333333333333333333333333333333333d0, &
3.555555555555555555555555555555555555556d0, &
0.d0, 0d0, 0d0, 0d0/)
a(5,1:7)=(/2.95259868922420362749580856576741350403902d0, &
-11.59579332418838591678097850937357110196616d0, &
9.8228928516994360615759792714525224813290657d0, &
-0.29080932784636488340192043895747599451303155d0, &
0.d0, 0d0, 0d0/)
a(6,1:7)=(/2.846275252525252525252525252525252525252525d0, &
-10.757575757575757575757575757575757575757576d0, &
8.906422717743472460453592529064227177434725d0, &
! ^ Here, mistake compare with
! http://depa.fquim.unam.mx/amyd/archivero/DormandPrince_19856.pdf
0.278409090909090909090909090909090909090909d0, &
-0.273531303602058319039451114922813036020583d0, &
0.d0,0d0/)
a(7,1:7)=(/0.09114583333333333333333333333333333333333333d0, &
0.d0, &
0.4492362982929020664869721473495058400718778077d0, &
0.6510416666666666666666666666666666666666666667d0, &
-0.3223761792452830188679245283018867924528301887d0, &
0.13095238095238095238095238095238095238095238095d0, &
0.d0/)
b2(1:7)=(/0.0899131944444444444444444444444444444444444444d0, &
0.d0, &
0.45348906858340820604971548367774782869122491764d0, &
0.6140625d0, &
-0.271512382075471698113207547169811320754716981d0, &
0.089047619047619047619047619047619047619047619d0, &
0.025d0/)

b1(1:7)=a(7,1:7)

Rc(1:7)=(/0.001232638888888888888888888888888888888888888889d0,&
0.d0, &
-0.004252770290506139562743336328241988619347109913d0, &
0.036979166666666666666666666666666666666666666d0, &
-0.050863797169811320754716981132075471698113207547d0, &
0.0419047619047619047619047619047619047619047619d0, &
-0.025d0/)

key=0
allocate(tmp(1:N),K(1:s,1:N))
tmp=0d0; K=0d0

if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
endif

if(h.ge.abs(xbound-x))h=xbound-x

FLAG=1
if(abs(x-xbound).le.hmin)then
info=1
FLAG=0
endif

do while(FLAG.eq.1)
tx=x

if(info.ne.-1)then
K(1,1:N)=h*work0(1:N)
else
do i=1,N
work(i)=grk(N,x,y,i)
K(1,i)=h*work(i)
enddo
work0(1:N)=work(1:N)
info=0
endif

do j=2,s
tx=x+c(j)*h
tmp(1:N)=y(1:N)
do i=1,j-1
tmp(1:N)=tmp(1:N)+K(i,1:N)*a(j,i)
enddo
do i=1,N
work(i)=grk(N,tx,tmp,i)
K(j,i)=h*work(i)
enddo
enddo

!step 4
R=0d0
do i=1,N
R=R+(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i) &
+Rc(5)*K(5,i)+Rc(6)*K(6,i)+Rc(7)*K(7,i))**2
enddo
R=abs(sqrt(R)/h/dble(N))

Sy=0d0
do i=1,N
Sy=Sy+(y(i)*y(i))
enddo
Sy=dsqrt(Sy)
if(Sy.ge.1d0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
y(1:N)=y(1:N)+b1(1)*K(1,1:N)+b1(3)*K(3,1:N) &
+b1(4)*K(4,1:N)+b1(5)*K(5,1:N)+b1(6)*K(6,1:N)
FLAG=0
work0(1:N)=work(1:N)
endif

!step 6
! Avoid zero deviding.
if(R.ge.1d-20)then
delta=(err/(2d0*R))**0.2d0
else
delta=4d0
endif

!step 7
if(delta.le.0.1d0)then
!function changes dramatically.
h=0.1d0*h
elseif(delta.ge.4d0)then
!function changes loosely.
h=4d0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1d0,h)*hmin
key=1
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if

if(h.le.0d0.and.xbound-x.ge.0d0)then
info=1
FLAG=0
elseif(h.ge.0d0.and.xbound-x.le.0d0)then
info=1
FLAG=0
endif
enddo

if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')"Strange point between ",x-h," and ",x
info=-2
endif

deallocate(tmp,K)
return
end subroutine dDP45

等間隔(サブルーチンは上のものと同じなので省略)

program main
implicit none
integer::i,N,info,Nx
double precision::h,tol
double precision::xa,xb,tx
double precision,allocatable::y(:),x(:),work0(:)
double precision,external::grk

N=1
allocate(y(1:N),work0(1:N))

Nx=101
allocate(x(1:Nx))
xa=0d0
xb=10d0
do i=1,Nx
x(i)=(i-1)*(xb-xa)/dble(Nx-1)+xa
enddo

!initial conditions
y(1)=1d0 ! x (0)

tol=1d-8
write(10,'(2e25.10e3)')x(1),y(1)
do i=2,Nx
info=-1
h=x(i)-x(i-1)
tx=x(i-1)
do while(info.le.0)
call dDP45(grk,tx,h,N,y,x(i),info,tol,work0)
enddo
write(10,'(2e25.10e3)')x(i),y(1)
enddo

stop
end program main

function grk(N,x,y,s)
implicit none
integer,intent(in)::N,s
double precision,intent(in)::x
double precision,intent(in)::y(1:N)
double precision::grk

grk=0d0
if(s.eq.1)then
grk=y(1)<em>cos(x)
else
write(6,</em>)"***Error grk"; stop
endif

return
end function grk

4倍精度ルーチン

4倍精度のサブルーチンです。
計算速度は倍精度の30~50倍かかるので、必要なとき以外使わないようにしましょう。

▼ここクリックでこの場に展開

subroutine qDP45(grk,x,h,N,y,xbound,info,tol,work0)
implicit none
integer,intent(in)::N
real*16,intent(in)::xbound,tol
real*16,intent(inout)::x,h,y(1:N),work0(1:N)
integer,intent(inout)::info
real*16,external::grk

! Runge-Kutta-Dormand-Prince method
!
! if h < hmin, propagate forcibly with warning.
!
!-----------------
!work0(1:N) : Just prepare, used and substituted by
! the coefficients K(7,1:7) to reduce computation cost
!info = -1 (Used for first loop;
! you should set info=0 when you try to use this routine)
! = 0 (Running now)
! = 1 (x reach xbound)
! = -9 (Path the discontinue points)
!-----------------
!

integer::i,j,FLAG,key
real*16::R,delta,tx,Sy,err
real*16,allocatable::tmp(:),K(:,:)
real*16,parameter::hmin=1q-20,hmax=0.5q0
integer,parameter::s=7
real*16::a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s),work(1:N)

c(1:7)=(/0q0, 0.2q0, 0.3q0,0.8q0,&
0.888888888888888888888888888888888888888889q0,&
1.q0, 1.q0/)
a(1:7,1:7)=0q0
a(1,1:7)=(/0q0, 0q0, 0q0, 0q0, 0q0, 0q0, 0q0/)
a(2,1:7)=(/0.2q0, 0q0, 0q0, 0q0, 0q0, 0q0, 0q0/)
a(3,1:7)=(/0.075q0, 0.225q0, 0q0, 0q0, 0q0, 0q0, 0q0/)
a(4,1:7)=(/0.977777777777777777777777777777777777778q0, &
-3.733333333333333333333333333333333333333q0, &
3.555555555555555555555555555555555555556q0, &
0.q0, 0q0, 0q0, 0q0/)
a(5,1:7)=(/2.95259868922420362749580856576741350403902q0, &
-11.59579332418838591678097850937357110196616q0, &
9.8228928516994360615759792714525224813290657q0, &
-0.29080932784636488340192043895747599451303155q0, &
0.q0, 0q0, 0q0/)
a(6,1:7)=(/2.846275252525252525252525252525252525252525q0, &
-10.757575757575757575757575757575757575757576q0, &
8.906422717743472460453592529064227177434725q0, &
! ^ Here, mistake compare with
! http://depa.fquim.unam.mx/amyd/archivero/DormandPrince_19856.pdf
0.278409090909090909090909090909090909090909q0, &
-0.273531303602058319039451114922813036020583q0, &
0.q0,0q0/)
a(7,1:7)=(/0.09114583333333333333333333333333333333333333q0, &
0.q0, &
0.4492362982929020664869721473495058400718778077q0, &
0.6510416666666666666666666666666666666666666667q0, &
-0.3223761792452830188679245283018867924528301887q0, &
0.13095238095238095238095238095238095238095238095q0, &
0.q0/)
b2(1:7)=(/0.0899131944444444444444444444444444444444444444q0, &
0.q0, &
0.45348906858340820604971548367774782869122491764q0, &
0.6140625q0, &
-0.271512382075471698113207547169811320754716981q0, &
0.089047619047619047619047619047619047619047619q0, &
0.025q0/)

b1(1:7)=a(7,1:7)

Rc(1:7)=(/0.001232638888888888888888888888888888888888888889q0,&
0.q0, &
-0.004252770290506139562743336328241988619347109913q0, &
0.036979166666666666666666666666666666666666666q0, &
-0.050863797169811320754716981132075471698113207547q0, &
0.0419047619047619047619047619047619047619047619q0, &
-0.025q0/)

key=0
allocate(tmp(1:N),K(1:s,1:N))
tmp=0q0; K=0q0

if(abs(h).ge.hmax)then
h=sign(1q0,h)*hmax
endif

if(h.ge.abs(xbound-x))h=xbound-x

FLAG=1
if(abs(x-xbound).le.hmin)then
info=1
FLAG=0
endif

do while(FLAG.eq.1)
tx=x

if(info.ne.-1)then
K(1,1:N)=h*work0(1:N)
else
do i=1,N
work(i)=grk(N,x,y,i)
K(1,i)=h*work(i)
enddo
work0(1:N)=work(1:N)
info=0
endif

do j=2,s
tx=x+c(j)*h
tmp(1:N)=y(1:N)
do i=1,j-1
tmp(1:N)=tmp(1:N)+K(i,1:N)*a(j,i)
enddo
do i=1,N
work(i)=grk(N,tx,tmp,i)
K(j,i)=h*work(i)
enddo
enddo

!step 4
R=0q0
do i=1,N
R=R+(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i) &
+Rc(5)*K(5,i)+Rc(6)*K(6,i)+Rc(7)*K(7,i))**2
enddo
R=abs(sqrt(R)/h/dble(N))

Sy=0q0
do i=1,N
Sy=Sy+(y(i)*y(i))
enddo
Sy=sqrt(Sy)
if(Sy.ge.1q0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
y(1:N)=y(1:N)+b1(1)*K(1,1:N)+b1(3)*K(3,1:N) &
+b1(4)*K(4,1:N)+b1(5)*K(5,1:N)+b1(6)*K(6,1:N)
FLAG=0
work0(1:N)=work(1:N)
endif

!step 6
! Avoid zero deviding.
if(R.ge.1q-50)then
delta=(err/(2q0*R))**0.2q0
else
delta=4q0
endif

!step 7
if(delta.le.0.1q0)then
!function changes dramatically.
h=0.1q0*h
elseif(delta.ge.4q0)then
!function changes loosely.
h=4q0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1q0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1q0,h)*hmin
key=1
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if

if(h.le.0q0.and.xbound-x.ge.0q0)then
info=1
FLAG=0
elseif(h.ge.0q0.and.xbound-x.le.0q0)then
info=1
FLAG=0
endif
enddo

if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')"Strange point between ",x-h," and ",x
info=-9
endif

deallocate(tmp,K)
return
end subroutine qDP45

subroutine qrkf45(grk,x,h,N,y,xbound,info,tol)
! if h < hmin, propagate forcibly with warning.
!
!-----------------
!info = -9 (maybe path the discontinue points)
! = 0 (Running now)
! = 1 (x reach xbound)
!-----------------
!
implicit none
integer,intent(in)::N
real*16,intent(in)::xbound,tol
real*16,intent(inout)::x,h,y(1:N)
integer,intent(inout)::info
real*16,external::grk

integer::i,j,FLAG,key
real*16::R,delta,tx,Sy,err
real*16,allocatable::tmp(:),K(:,:)
real*16,parameter::hmin=1q-20,hmax=0.5q0
integer,parameter::s=6
real*16::a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s)

c(1:6)=(/0q0, 0.25q0, 0.375q0,&
0.9230769230769230769230769230769230769231q0, 1q0, 0.5q0/)
a(1:6,1:6)=0q0
a(1,1:6)=(/0q0, 0q0, 0q0, 0q0, 0q0, 0q0/)
a(2,1:6)=(/0.25q0, 0q0, 0q0, 0q0, 0q0, 0q0/)
a(3,1:6)=(/0.09375q0, 0.28125q0, 0q0, 0q0, 0q0, 0q0/)
a(4,1:6)=(/0.8793809740555302685480200273099681383705q0, &
-3.277196176604460628129267182521620391443q0, &
3.320892125625853436504324078288575329995q0, 0q0, 0q0, 0q0/)
a(5,1:6)=(/2.032407407407407407407407407407407407407q0,-8q0, &
7.173489278752436647173489278752436647173q0, &
-0.2058966861598440545808966861598440545809q0, 0q0, 0q0/)
a(6,1:6)=(/-0.2962962962962962962962962962962962962963q0,2q0, &
-1.381676413255360623781676413255360623782q0, &
0.4529727095516569200779727095516569200780q0,-0.275q0,0q0/)
b2(1:6)=(/0.1185185185185185185185185185185185185185q0, 0.q0,&
0.5189863547758284600389863547758284600390q0, &
0.5061314903420166578061314903420166578061q0, &
-0.18q0, 0.03636363636363636363636363636363636363636q0/)
b1(1:6)=(/0.1157407407407407407407407407407407407407q0, 0q0,&
0.5489278752436647173489278752436647173489q0, &
0.5353313840155945419103313840155945419103q0, -0.2q0, 0q0/)
Rc(1:6)=(/0.002777777777777777777777777777777777777778q0,0q0, &
-0.02994152046783625730994152046783625730994q0, &
-0.02919989367357788410419989367357788410420q0, 0.02q0, &
0.03636363636363636363636363636363636363636q0/)

if(info.eq.-1)info=0

key=0
allocate(tmp(1:N),K(1:s,1:N))
tmp=0q0; K=0q0

if(abs(h).ge.hmax)then
h=sign(1q0,h)*hmax
endif

if(h.ge.abs(xbound-x))h=xbound-x

FLAG=1
if(abs(x-xbound).le.hmin)then
info=1
FLAG=0
endif

do while(FLAG.eq.1)

tx=x
do j=1,s
tx=x+c(j)*h
tmp(1:N)=y(1:N)
do i=1,j-1
tmp(1:N)=tmp(1:N)+K(i,1:N)*a(j,i)
enddo
do i=1,N
K(j,i)=h*grk(N,tx,tmp,i)
enddo
enddo

!step 4
R=0q0
do i=1,N
R=R+(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i)+Rc(5)*K(5,i)+Rc(6)*K(6,i))**2q0
enddo
R=abs(sqrt(R)/h/(N*1q0))

Sy=0q0
do i=1,N
Sy=Sy+(y(i)*y(i))
enddo
Sy=sqrt(Sy)
if(Sy.ge.1q0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
y(1:N)=y(1:N)+b1(1)*K(1,1:N)+b1(3)*K(3,1:N)+b1(4)*K(4,1:N)+b1(5)*K(5,1:N)
FLAG=0
endif

!step 6
! Avoid zero deviding.
if(R.ge.1q-50)then
delta=(err/(2q0*R))**0.25q0
else
delta=4q0
endif

!step 7
if(delta.le.0.1q0)then
!function changes dramatically.
h=0.1q0*h
elseif(delta.ge.4q0)then
!function changes loosely.
h=4q0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1q0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1q0,h)*hmin
key=1
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if

if(h.le.0q0.and.xbound-x.ge.0q0)then
info=1
FLAG=0
elseif(h.ge.0q0.and.xbound-x.le.0q0)then
info=1
FLAG=0
endif
enddo

if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')"Strange point between ",x-h," and ",x
info=-9
endif

deallocate(tmp,K)
return
end subroutine qrkf45

陽的ルンゲ=クッタ法の導出

ルンゲ=クッタ法の導出は煩雑です。単に複雑なだけです。
導出過程について詳しく述べられているページは、早川様が記述された以下のpdfを参照すると良いと思います。
Runge-Kutta法についてのノート(早川尚男)
計算過程を含め記述されているので分かりやすいです。

参考文献

[1]Derivation of Runge–Kutta methods
[2]NDSolveの”ExplicitRungeKutta”メソッド
[3]High-Order Explicit Runge-Kutta Methods
[4]List of Runge–Kutta methods
[5]Runge-Kutta-Fehlberg Method (RKF45)
[6]Runge-Kutta-Fehlberg method
[7]Lecture:13Runge-Kutta-Fehlberg Method
[8]GPU acceleration of Runge Kutta-Fehlberg and its comparison with Dormand-Prince method
[9]William H. Pressら著『ニューメリカルレシピ・イン・シー日本語版―C言語による数値計算のレシピ– 』(1993)

ニューメリカルレシピ・イン・シー日本語版―C言語による数値計算のレシピ

↑この本が一番有益だと思います。Fortran版もありますが、英語しかありません。ちなみに、英語で,若干古いバージョンでいいのならば
公式のホームページ
Numerical Recipes in C (1992)
Numerical Recipes in Fortran 77 and Fortran 90 (1992, 1996)
で無料で公開されています。

プログラミングと数値計算

数値積分(等間隔)

2015年4月7日 sikino コメントする

数値計算での積分方法、特に等間隔の分点の場合であるニュートン・コーツ積分(とロンバーグ積分)に関する理論とのプログラムを載せます。

ニュートンコーツ型の公式
ルンゲ現象
ロンバーグ積分
サブルーチン”integral”のコードと使い方
- サブルーチン”integral”を用いた例題
サブルーチン”romberg”のコードと使い方
- サブルーチン”romberg”を用いた例題

等間隔の分点で積分を行う場合、よく使われる方法は台形則、もしくはシンプソン積分です。

台形則

台形則は、分点間を台形近似して面積を求める方法であり、下図のようなイメージで積分を近似する方法です。

積分$\displaystyle \int_a^b f(x) dx$を、
$
\displaystyle \int_a^b f(x) dx \sim \sum_{i=0}^{N-1} \frac{f(x_{i+1})+f(x_i)}{2}h
$
として近似します。この時、誤差のオーダーは$O(h^3)$となります。
言い換えると、
台形積分とは、隣り合う分点間を1次関数で近似して求める積分
言えます。

シンプソン積分

シンプソン積分は、隣り合う分点間を2次関数で近似して求める積分です。イメージは

な感じです。数式では
$
\displaystyle \int_a^b f(x) dx \sim \sum_{i=0}^{N/2} \frac{f(x_{2i})+4f(x_{2i+1})+f(x_{2i+2})}{3}h
$
となります。

高次≠高精度

これ以降もあります。分点間を3次、4次、…としていけば高次の積分をすることが可能になります。
この分点間の次数をどんどん上げていって積分をする方法をまとめてニュートンコーツの公式、と呼びます。
例えばの台形則はニュートンコーツの公式の1次に相当し、シンプソン則はニュートンコーツの公式の2次に相当しています。

じゃあ100次の公式作れば精度が凄い公式ができるんじゃないか？となりますが、これは違います。
高次であることと精度が高いことは違うのです。だから使うな、というわけではなくて知っておいてほしいことなんです。
ではなぜか？これはwikipediaのニュートンコーツの公式のページからとってきたものですが、この式に理由があります。

問題となるのが誤差項です。例えば台形近似を見ますと、誤差項は$h^3 f^{(2)}(\xi)$とあります。hは刻み幅なので、この項は刻み幅を小さくすれば小さくするほど精度がよくなることを言っています。
が、問題は$f^{(2)}(\xi)$の方です。この項は本当の関数が分からないと評価できません。
多項式だったら微分するごとに次数が1つづつ減少するのでこの項は高次になれば小さくなります。じゃあもしも微分するごとに値が増加していったら…？この場合、hを小さくしてもそれを上回るくらい$f^{(2)}(\xi)$が大きくなることがあります。これをルンゲ現象と呼びます。

ルンゲ現象

例えば関数
$
\displaystyle f(x)=\frac{1}{1+25x^2}
$
を考えます。wikipediaのページを参考にすれば、その導関数と値は、

となります。高次になればなるほど値が増加する関数なのです。
実際には$h^n f^{(m)}(\xi)$の積の値が小さいか大きいか？なので、一概には高次は悪い！とは言えません。高次は危ないくらいでしょう。

これを可能な限り回避するために、低次の積分を組み合わせて使う方法がとられます。
こうすることで高次の値を増加を抑えるのです。

そのほかの手法もいろいろあります。分点の位置に縛られなければガウス求積法が最高でしょう。素晴らしい方法です。ぜひ調べてみてください。

[adsense1]

ロンバーグ積分

さて、高次を追い求めるのもいいですが、一線を駕す…とまではいきませんが、ロンバーグ(Romberg)積分というものがあります。
このロンバーグ積分とは、台形則による計算を基本とし、無限小区間での積分結果を補外によって求めよう、というものです。
ここでいう補外の概念は以下のような感じです。

刻み幅1で台形則の結果から積分値S1が分かる
↓
刻み幅0.5で台形則による結果から積分値S2がわかる
↓
刻み幅0.25で台形則による結果から積分値S3がわかる
↓
刻み幅0.175で台形則による結果から積分値S4がわかる

とある程度計算していきます。そうすると刻み幅hを変化させるとそれに応じて積分値がS1→S2→S3→S4…と推移していくわけです。
この推移の仕方を計算すれば、刻み幅無限小の場合の積分値Sが分かる、こういう仕組みになっています。収束加速法と言ったり、リチャードソン補外、なんて言ったりします。
この補外による方法はかなり優れています。先ほどの高次がどうこう、といった問題が発生しないのです。台形則しか使ってないんですから。有限桁計算に適した方法である、とかwikipediaに載ってたりします。

ロンバーグ積分の詳しい計算方法や理論は
ニューメリカルレシピ・イン・シー日本語版―C言語による数値計算のレシピ

ニューメリカルレシピ・イン・シー日本語版―C言語による数値計算のレシピ

を見ると詳しく書かれています。

日本語では第6章数値積分と数値微分にロンバーグ積分の記述があります。
英語ではRomberg integralで調べると良いでしょう。参考までに、英語では
Romberg integration exampleや、Romberg Integrationなにかがいいと思います。

2016/03/01
上のRomberg Integrationを参考にして書いたプログラムはこちらです。これは、解析関数fをRomberg積分するプログラムです。
解析関数ではなく、サイズ$2^{n}+1$の配列に格納されている場合はこのページの下の方にあるプログラムを参照してください。

call romberg(a,b,s,pre)

で呼び出し、区間a~bの積分$s=\int_a^b f(x) dx$を精度preで求めるプログラムです。

▼クリックでこの場に展開されます

サブルーチン”integral”のコードと使い方

下のサブルーチン”integral”は1,2,3次元の短冊近似、台形則、シンプソン積分、シンプソン3/8積分、ブール則(ニュートン・コーツ型積分)をカバーします。
ただし、分点の個数に注意してください。分点の個数を気にせず使える近似方法は、短冊近似、台形近似のみです。

ロンバーグ積分では引数の関係上、ニュートン・コーツ型積分との間をうまく処理できなかったため、別のサブルーチンとして書きます。もしもロンバーグ積分を使いたい場合はもう少し先に進んでください。

まずはニュートン・コーツ型積分です。
下のサブルーチンを使ってください。
使い方は、

call integral(配列の大きさ,配列,刻み幅,積分結果,積分手法)

という感じです。
実際のプログラムでは、以下のように使用してください。配列yは複素数解列(ただし実軸上の積分)でもokです。

1次元

サブルーチンintegralの使い方
積分則	呼び方
短冊近似(長方形近似)	call integral(size(w,1),w,h,s,"box")	配列yの大きさに指定はない。
台形近似	call integral(size(y,1),y,h,s,"trapezoid")	配列yの大きさに指定はない。
シンプソン則	call integral(size(y,1),y,h,s,"simpson")	配列yの大きさが3,5,7,…個、$2n+1$個でないといけない
シンプソン3/8則	call integral(size(y,1),y,h,s,"simpson38")	配列yの大きさが4,7,10,…個、$3n+1$個でないといけない
ブール則	call integral(size(y,1),y,h,s,"boole")	配列yの大きさが5,9,13,…個、$4n+1$個でないといけない

2次元配列の場合は以下のように指定してください。

call integral(size(w,1),size(w,2),w,hx,hy,s,"simpson")

3次元の場合はこう指定してください。

call integral(size(w,1),size(w,2),size(w,3),w,hx,hy,hz,s,"simpson")

積分のためのモジュール

module integral_mod
!developer --> sikino
!date --> 2015/04/07
implicit none
interface integral
module procedure &
dintegral, &
dintegral2d, &
dintegral3d, &
cintegral, &
cintegral2d, &
cintegral3d
end interface integral
contains
subroutine dintegral(N,y,h,s,method)
integer,intent(in)::N
double precision,intent(in)::h,y(1:N)
character(*),intent(in)::method
double precision,intent(out)::s
integer::i
double precision::y0,y1,y2,y3

s=0.d0; y0=0.d0; y1=0.d0; y2=0.d0; y3=0.d0

if(trim(method).eq."box")then
s=h*sum(y(1:N-1))
elseif(trim(method).eq."trapezoid")then
y1=y(1)+y(N)
do i=2,N-1
y2=y2+y(i)
enddo
s=(y1+2.d0*y2)*h*0.5d0
elseif(trim(method).eq."simpson")then
if(mod(N,2).ne.1)then
write(6,*)"=====cannot calculation with simpson"
write(6,*)"=====program stop"
stop
endif

y1=y(1)+y(N)
do i=2,N-1,2
y2=y2+y(i)
enddo
do i=3,N-2,2
y3=y3+y(i)
enddo

s=(y1+4.d0*y2+2.d0*y3)*h/3.d0

elseif(trim(method).eq."simpson38")then

if(mod(N,3).ne.1)then
write(6,*)"=====cannot calculation with simpson38"
write(6,*)"=====program stop"
stop
endif

y0=y(1)+y(N)
do i=2,N-2,3
y1=y1+y(i)
enddo
do i=3,N-1,3
y2=y2+y(i)
enddo
do i=4,N-3,3
y3=y3+y(i)
enddo
s=(y0+3.d0*(y1+y2)+2.d0*y3)*3.d0*h/8.d0

elseif(trim(method).eq."boole")then
if(mod(N,4).ne.1)then
write(6,*)"=====cannot calculation with boole"
write(6,*)"=====program stop"
stop
endif

y0=y(1)+y(N)
do i=5,N-4,4
y1=y1+y(i)
enddo
do i=2,N-1,2
y2=y2+y(i)
enddo
do i=3,N-2,4
y3=y3+y(i)
enddo

s=(14.d0*y0+28.d0*y1+64.d0*y2+24.d0*y3)*h/45.d0
else
write(6,*)"=====cannot calculation in integral"
write(6,*)"=====program stop"
stop
end if

return
end subroutine dintegral

subroutine dintegral2d(Nx,Ny,z,hx,hy,s,method)
implicit none
integer,intent(in)::Nx,Ny
double precision,intent(in)::hx,hy,z(1:Nx,1:Ny)
character(*),intent(in)::method
double precision,intent(out)::s
integer::i
double precision::ty(1:Ny),r(1:Nx)

s=0.d0
ty(1:Ny)=0.d0
r(1:Nx)=0.d0
do i=1,Nx
ty(1:Ny)=z(i,1:Ny)
call integral(Ny,ty,hy,s,method)
r(i)=s
enddo
call integral(Nx,r,hx,s,method)

return
end subroutine dintegral2d

subroutine dintegral3d(Nx,Ny,Nz,w,hx,hy,hz,s,method)
implicit none
integer,intent(in)::Nx,Ny,Nz
double precision,intent(in)::hx,hy,hz,w(1:Nx,1:Ny,1:Nz)
character(*),intent(in)::method
double precision,intent(out)::s
integer::i
double precision::tyz(1:Ny,1:Nz),r(1:Nx)

s=0.d0
tyz(1:Ny,1:Nz)=0.d0
r(1:Nx)=0.d0
do i=1,Nx
tyz(1:Ny,1:Nz)=w(i,1:Ny,1:Nz)
call integral(Ny,Nz,tyz,hy,hz,s,method)
r(i)=s
enddo
call integral(Nx,r,hx,s,method)

return
end subroutine dintegral3d

subroutine cintegral(N,y,h,s,method)
integer,intent(in)::N
complex(kind(0d0)),intent(in)::y(1:N)
double precision,intent(in)::h
character(*),intent(in)::method
complex(kind(0d0)),intent(out)::s

double precision::res,ims

s=dcmplx(0d0,0d0); res=0.d0; ims=0.d0

call integral(N,dble(y),h,res,trim(method))
call integral(N,dimag(y),h,ims,trim(method))

s=dcmplx(res,ims)

return
end subroutine cintegral

subroutine cintegral2d(Nx,Ny,z,hx,hy,s,method)
integer,intent(in)::Nx,Ny
complex(kind(0d0)),intent(in)::z(1:Nx,1:Ny)
double precision,intent(in)::hx,hy
character(*),intent(in)::method
complex(kind(0d0)),intent(out)::s

double precision::res,ims

s=dcmplx(0d0,0d0); res=0.d0; ims=0.d0

call integral(Nx,Ny,dble(z),hx,hy,res,trim(method))
call integral(Nx,Ny,dimag(z),hx,hy,ims,trim(method))

s=dcmplx(res,ims)

return
end subroutine cintegral2d

subroutine cintegral3d(Nx,Ny,Nz,w,hx,hy,hz,s,method)
integer,intent(in)::Nx,Ny,Nz
complex(kind(0d0)),intent(in)::w(1:Nx,1:Ny,1:Nz)
double precision,intent(in)::hx,hy,hz
character(*),intent(in)::method
complex(kind(0d0)),intent(out)::s

double precision::res,ims

s=dcmplx(0d0,0d0); res=0.d0; ims=0.d0

call integral(Nx,Ny,Nz,dble(w),hx,hy,hz,res,trim(method))
call integral(Nx,Ny,Nz,dimag(w),hx,hy,hz,ims,trim(method))

s=dcmplx(res,ims)

return
end subroutine cintegral3d
end module integral_mod

サブルーチン”integral”を用いた例題

実際の使い方。参考にどうぞ。
下のコードは2次元の積分
$
\displaystyle \int\int xe^{-x^2}e^{-y}dxdy=\frac{e^{-x^2}}{2}e^{-y}+\mbox{const}
$
を数値積分します。積分範囲を$x=0\sim 3, \; y=0 \sim 5$にした場合、その解析解は
$
\begin{align}
\int_0^3 dx\int_0^5 dy xe^{-x^2}e^{-y}&=\frac{1}{2}(1-e^{-5}-e^{-9}+e^{-14})\\
&=0.4965697373627734784608751894356320535936864348993604\cdots
\end{align}
$
となります。実行すると

./a.out
simpson : 0.496569737366759E+000
romberg,3 : 0.496569737362773E+000

となります。コードは、

program main
use integral_mod
implicit none
integer::i,j,n,Nx,Ny
double precision::hx,hy,ans
double precision::xmin,xmax,ymin,ymax
double precision,allocatable::x(:),y(:)
double precision::s
double precision,allocatable::w(:,:)
double precision::f
external::f

n=10

xmin=0d0
xmax=3d0

ymin=0d0
ymax=5d0

Nx=2**(n)
Ny=2**(n)
allocate(x(0:Nx),y(0:Ny),w(0:Nx,0:Ny))
x=0d0; y=0d0; w=dcmplx(0d0,0d0)

hx=(xmax-xmin)/dble(Nx)
hy=(ymax-ymin)/dble(Ny)

do i=0,Nx
x(i)=xmin+hx*dble(i)
enddo
do i=0,Ny
y(i)=ymin+hy*dble(i)
enddo

do i=0,Nx
do j=0,Ny
w(i,j)=f(x(i),y(j))
end do
end do

call integral(size(w,1),size(w,2),w,hx,hy,s,"simpson")
write(6,'(A,e25.15e3)')"simpson : ",s

call integral(size(w,1),size(w,2),w,hx,hy,s,"romberg",3)
write(6,'(A,e25.15e3)')"romberg,3 : ",s

return
end program main

function f(x,y)
implicit none
double precision,intent(in)::x,y
double precision::f

f=x*exp(-x*x)*exp(-y)

return
end function f

3次元の計算の例題です。
下のコードは
$
\displaystyle \int\int\int x^4y^4z^4 dxdydz=\left(\frac{1}{5}\{0.7^5-(-1)^5\}\right)^3
$
を数値積分します。積分範囲を$x,y,z=-1\sim 0.7$にした場合、その値は
$
\displaystyle \int_{-1}^{0.7}\int_{-1}^{0.7}\int_{-1}^{0.7} x^4y^4z^4 dxdydz = 0.127496010896795\cdots
$
となります。実行すると、

$ gfortran integralmod.f90 main.f90
$ ./a.out
analysis : 0.127496010896795E-01
simpson : 0.10507613E-006 0.10507613E-006
romberg,3 : -0.17347235E-017 -0.17347235E-017

と得られます。中身はこうです。分点の数はx,y,z軸各々で違っていて構いません。

program main
use integral_mod
implicit none
integer::i,j,k,n,Nx,Ny,Nz
double precision::hx,hy,hz,ans
double precision,allocatable::x(:),y(:),z(:)
complex(kind(0d0))::s
complex(kind(0d0)),allocatable::w(:,:,:)

n=6
Nx=2**(n)
Ny=2**(n-1)
Nz=2**(n+1)

allocate(x(0:Nx),y(0:Ny),z(0:Nz),w(0:Nx,0:Ny,0:Nz))
x=0d0; y=0d0; z=0d0; w=dcmplx(0d0,0d0)

hx=1.7d0/dble(size(w,1)-1)
hy=1.7d0/dble(size(w,2)-1)
hz=1.7d0/dble(size(w,3)-1)

do i=0,Nx
x(i)=hx*dble(i)-1d0
enddo
do i=0,Ny
y(i)=hy*dble(i)-1d0
enddo
do i=0,Nz
z(i)=hz*dble(i)-1d0
enddo

do i=0,Nx
do j=0,Ny
do k=0,Nz
w(i,j,k)=dcmplx((x(i)*y(j)*z(k))**4,(x(i)*y(j)*z(k))**4)
end do
end do
end do

ans=(0.7d0**5.d0+1.d0)/5.d0
ans=ans**3.d0
write(6,'(A,e23.15e2)')"analysis : ",ans

call integral(size(w,1),size(w,2),size(w,3),w,hx,hy,hz,s,"simpson")
write(6,'(A,2e17.8e3)')"simpson : ",dble(s)-ans,dimag(s)-ans

call integral(size(w,1),size(w,2),size(w,3),w,hx,hy,hz,s,"romberg",3)
write(6,'(A,2e17.8e3)')"romberg,3 : ",dble(s)-ans,dimag(s)-ans

return
end program main

ロンバーグ積分のプログラム

ロンバーグ積分

call romberg(jx,x,y,s,pre)

配列yの大きさは2,3,5,10,…,$2^{jx}+1$個でないといけない

ロンバーグ積分では収束精度”pre”を指定すること。もしも収束精度に達しなかった場合、警告と共に、与えられた値での収束限界の積分結果を返す。積分精度を高めたければ刻み幅を小さくする操作(例えば分点数を増やす等)をしてください。”pre”は、積分結果の絶対値が1より大きくなる場合は相対誤差を、小さくなる場合は絶対誤差を取ります。この理由は以前の結果の補正を加え、収束をさせるためであり、非常に小さい誤差の場合いつまでたっても機械誤差のため収束しないからです。

module romberg_mod
!developer --> sikino
!date --> 2016/03/01
! 2016/03/03
!
! 1D case :
! romberg(jx,x,y,s,pre)
! | | | | +- (in) precision (e.g. 1d-8)
! | | | +--- (out) integration result
! | | +----- (in) y(1:2**jx+1) f(x)
! | +------- (in) x(1:2**jx+1) x
! +---------- (in) related to array size
!
! 2D case :
! romberg(jx,jy,x,y,z,s,pre)
! | | | | | | +-- (in) precision (e.g. 1d-8)
! | | | | | +---- (out) integration result
! | | | | +------ (in) z(1:2**jx+1,1:2**jy+1) f(x,y)
! | | | +-------- (in) y(1:2**jy+1) y
! | | +---------- (in) x(1:2**jx+1) x
! | +------------- (in) related to array size y
! +---------------- (in) related to array size x
!
! 3D case :
! romberg(jx,jy,jz,x,y,z,w,s,pre)
! | | | | | | | | +-- (in) precision (e.g. 1d-8)
! | | | | | | | +---- (out) integration result
! | | | | | | +------ (in) w(1:2**jx+1,1:2**jy+1,1:2**jz+1) f(x,y,z)
! | | | | | +-------- (in) z(1:2**jz+1) z
! | | | | +---------- (in) y(1:2**jy+1) y
! | | | +------------ (in) x(1:2**jx+1) x
! | | +--------------- (in) related to array size z
! | +------------------ (in) related to array size y
! +--------------------- (in) related to array size x
!
implicit none
interface romberg
module procedure &
dromberg, &
dromberg2d, &
dromberg3d
end interface romberg
contains
subroutine dromberg(jx,x,y,s,pre)
implicit none
integer,intent(in)::jx
double precision,intent(in)::x(1:2**jx+1),y(1:2**jx+1),pre
double precision,intent(out)::s

integer,parameter::jm=6 !--> precision: O(h^(2*jm))
integer,parameter::nm=2**jm
double precision,allocatable::tx(:),ty(:)
double precision::ts
integer::k

s=0d0
if(jx.ge.jm)then
allocate(tx(1:nm+1),ty(1:nm+1)); tx=0d0; ty=0d0
do k=0,2**(jx-jm)-1
ts=0d0
tx(1:nm+1)=x(k*nm+1:(k+1)*nm+1)
ty(1:nm+1)=y(k*nm+1:(k+1)*nm+1)

call romberg_sub(jm,tx,ty,ts,pre)
s=s+ts
enddo
deallocate(tx,ty)
else
call romberg_sub(jx,x,y,s,pre)
endif

return
end subroutine dromberg

subroutine romberg_sub(jx,x,y,s,pre)
! reference "http://www.math.usm.edu/lambers/mat460/fall09/lecture29.pdf"
implicit none
integer,intent(in)::jx
double precision,intent(in)::x(1:2**jx+1),y(1:2**jx+1),pre
double precision,intent(out)::s

integer::i,j,k,n,dn
double precision::h,ps,tmp
double precision::T(1:jx+1,1:jx+1)

n=2**jx+1

h=x(n)-x(1)
dn=(n-1)/2

T(1,1)=0.5d0*h*(y(1)+y(n))
s=T(1,1)
ps=s
h=0.5d0*h
do j=2,jx+1

! trapezoidal rule
tmp=0d0
do i=1,2**(j-2)
tmp=tmp+y(1+(2*i-1)*(dn))
enddo
T(j,1)=0.5d0*T(j-1,1)+h*tmp

do k=2,j
! Richardson extrapolation
T(j,k)=T(j,k-1)+(T(j,k-1)-T(j-1,k-1))/(dble(4**(k-1))-1d0)
enddo
s=T(j,j)

! precision check
if(abs(s).ge.1d0)then
if(abs((ps-s)/s).le.pre)exit
else
if(abs((ps-s)).le.pre)exit
endif
ps=s
h=0.5d0*h
dn=dn/2
enddo

if(j-1.eq.jx)then
write(6,'(A)')" -+-+- didn't converge at romberg integral -+-+- "
write(6,'(A)')" Please change stepsize h of array x "
endif

return
end subroutine romberg_sub

subroutine dromberg2d(jx,jy,x,y,z,s,pre)
implicit none
integer,intent(in)::jx,jy
double precision,intent(in)::pre,x(1:2**jx+1),y(1:2**jy+1)
double precision,intent(in)::z(1:2**jx+1,1:2**jy+1)
double precision,intent(out)::s

integer::i,nx,ny
double precision::ty(1:2**jy+1),r(1:2**jx+1)
nx=2**jx+1
ny=2**jy+1

s=0.d0
ty(1:ny)=0.d0
r(1:nx)=0.d0
do i=1,nx
ty(1:ny)=z(i,1:ny)
call romberg(jy,y,ty,s,pre)
r(i)=s
enddo
call romberg(jx,x,r,s,pre)

return
end subroutine dromberg2d

subroutine dromberg3d(jx,jy,jz,x,y,z,w,s,pre)
implicit none
integer,intent(in)::jx,jy,jz
double precision,intent(in)::pre,x(1:2**jx+1),y(1:2**jy+1),z(1:2**jz+1)
double precision,intent(in)::w(1:2**jx+1,1:2**jy+1,1:2**jz+1)
double precision,intent(out)::s
integer::i,nx,ny,nz
double precision::tyz(1:2**jy+1,1:2**jz+1),r(1:2**jx+1)

nx=2**jx+1; ny=2**jy+1; nz=2**jz+1
s=0.d0
tyz(1:ny,1:nz)=0.d0
r(1:nx)=0.d0
do i=1,nx
tyz(1:ny,1:nz)=w(i,1:ny,1:nz)
call romberg(jy,jz,y,z,tyz,s,pre)
r(i)=s
enddo
call romberg(jx,x,r,s,pre)

return
end subroutine dromberg3d
end module romberg_mod

ロンバーグ積分(romberg_mod)の例

必要なのは、上で紹介したモジュール “romberg_mod” と下のメインプログラムです。
以下のプログラムは1次元の定積分
$
\int_1^10 \frac{1}{x^2} dx = 0.9
$
を分点数$2^8+1(jx=8)$個でロンバーグ積分するものです。
精度は$O(h^{2\cdot min{jm,jx}})$です。
ここで$jm$はモジュールromberg_modの中にパラメータとして宣言されています。

コンパイルは

ifort romberg_mod.f90 main.f90

とでもすればいいでしょう。

program main
use romberg_mod
implicit none
integer::i,jx,n
double precision::h,a,b,s
double precision,allocatable::x(:),y(:)
double precision::f
external::f

a=1d0
b=10d0
jx=8
n=2**jx+1
allocate(x(1:n),y(1:n))
h=(b-a)/dble(n-1)
do i=1,n
x(i)=a+h*dble(i-1)
y(i)=f(x(i))
enddo

call romberg(jx,x,y,s,1d-8)
write(6,*)s,"romberg"

stop
end program main

function f(x)
implicit none
double precision::f,x

f=1d0/(x*x)

return
end function f

実行例(要求精度は$10^{-8}$)

>./a.out
0.900000000062642
>

2次元ではこう。

program main
use romberg_mod
implicit none
integer::i,j,Nx,Ny,jx,jy
double precision::hx,hy
double precision::xmin,xmax,ymin,ymax
double precision,allocatable::x(:),y(:)
double precision::s
double precision,allocatable::w(:,:)
double precision::f
external::f

jx=10
jy=8
xmin=0d0
xmax=3d0

ymin=0d0
ymax=5d0

Nx=2**(jx)
Ny=2**(jy)
allocate(x(0:Nx),y(0:Ny),w(0:Nx,0:Ny))
x=0d0; y=0d0; w=dcmplx(0d0,0d0)

hx=(xmax-xmin)/dble(Nx)
hy=(ymax-ymin)/dble(Ny)

do i=0,Nx
x(i)=xmin+hx*dble(i)
enddo
do i=0,Ny
y(i)=ymin+hy*dble(i)
enddo

do i=0,Nx
do j=0,Ny
w(i,j)=f(x(i),y(j))
end do
end do

call romberg(jx,jy,x,y,w,s,1d-8)
write(6,*)s

stop
end program main

function f(x,y)
implicit none
double precision,intent(in)::x,y
double precision::f

f=x*exp(-x*x)*exp(-y)

return
end function f

実行例

>./a.out
0.496569737374163
>

[adsense2]

プログラミングと数値計算

ファイル読み込み

2015年3月31日 sikino コメントする

ファイルの読み込みに関するサブルーチンを記述します。

do i=1,4
do j=1,83
write(30,*)i,j,i*j
end do
write(30,*)
end do

によってファイル”fort.30″が作られたとします。
今、fort.30を読み込んで
“i”の数”4”
と
“j”の数”83″
を取り出したいとします。

read文を使って読み込みますが、そのまま読み込むと空白部分を読み込んだり、読み込まないだったりします。
ここでは上のdoループによって作られたファイルの場合に使えるものを想定します。

ここでは、
大きな塊を表す数である”4″をblock,
塊の中の要素数を表す”83″をelement
と呼ぶことにします。

この問題を解く戦略は、ファイルを読み込む2種類の方法です。
1つは本当の行数(空行含む)を数え、もう1つは空行を飛ばして読み込む方法です。
本当の行数を数えるサブルーチンは下の”linecount”
であり、
空行を読み飛ばして行数を数えるサブルーチンは下の”linecount_eff”
です。
この二つと、一番下の行に追加される余分な１行を対処するために”breaklinecheck”というルーチンを使います。

これによってblockとelementを出力するサブルーチン”blockelement”を作っています。

下のプログラムを実行すると

$ gfortran main.f90
$ ./a.out
===Apply for fortran file will done===
Nblock :: 4
Nelement :: 83

という結果が得られるかと思います。

program main
!developer => sikinote
!date => 2015/03/31
implicit none
integer::NBlock,Nelement
character(48)::filename

filename='./fort.30'
call blockelement(filename,Nblock,Nelement)

write(6,'(A,i10)')" Nblock ::",Nblock
write(6,'(A,i10)')" Nelement ::",Nelement

stop
end program
!===================================

subroutine blockelement(filename,Nblock,Nelement)
!developer => sikino
!date => 2015/03/31
implicit none
character(*),intent(in)::filename
integer,intent(out)::Nblock,Nelement
integer::c1,c2

call linecount(c1,filename)
call linecount_eff(c2,filename)
call breaklinecheck(c1,c2)
Nblock=c1-c2+1
Nelement=c2/Nblock

return
end subroutine blockelement
!------------------------------
subroutine linecount(c,filename)
implicit none
integer,intent(out)::c
character(*),intent(in)::filename

integer::ier
character(len_trim(filename))::fname

fname=trim(filename)
c=0
open(100,file=fname,status='old',iostat=ier,err=990)
do while(.true.)
read(100,*,end=999)
c=c+1
enddo
999 continue
close(100)

return
!+-----------------------+
990 write(6,'(A)')"!!!!==error when open file",trim(fname),"info-->",ier
write(6,*)"======program stop at linecount"
stop

end subroutine linecount
!--------------------------------
subroutine linecount_eff(c,filename)
implicit none
integer,intent(out)::c
character(*),intent(in)::filename

integer::ier
character(100)::cc
character(len_trim(filename))::fname

fname=trim(filename)

c=0
open(100,file=fname,status='old',iostat=ier,err=990)
do while(.true.)
read(100,*,end=998)cc
if(len_trim(cc).gt.0)c=c+1
enddo

998 continue
close(100)
return

990 write(6,'(A)')"!!!!==error when open file",trim(fname),"info==>",ier
write(6,*)"======program stop at linecount_eff"
stop

end subroutine linecount_eff
!-------------------------------------
subroutine breaklinecheck(c1,c2)
implicit none
integer,intent(inout)::c1
integer,intent(in)::c2
integer::Nb

Nb=c1-c2+1
if(Nb.eq.2.or.mod(c2,Nb).ne.0)then
write(6,*)"===Apply for fortran file will done==="
c1=c1-1
Nb=c1-c2+1
if(mod(c2,Nb).ne.0)then
write(6,*)"line is different(may be last break)"
write(6,*)"program stop at subroutine __breaklinecheck__"
stop
end if
endif

return
end subroutine breaklinecheck

データを読み込むには？

さて、ブロックの数と要素の数が上のサブルーチンを使うことにより求められることがわかりました。
実際にデータを配列に代入するためにはどうすればいいんでしょう？
型に応じて使うサブルーチンを変えます。
その手続きは下のモジュールを記述することでokです。これを書いた上で、
メインプログラムを以下のように書きます。そうすれば配列x(,)とy(,)に値がはいります。

program main
use read1
implicit none
integer::NBlock,Nelement,i,j
character(48)::filename

double precision,allocatable::x(:,:),y(:,:)

filename='./fort.30'
call blockelement(filename,Nblock,Nelement)

write(6,'(A,i10)')" Nblock ::",Nblock
write(6,'(A,i10)')" Nelement ::",Nelement

allocate(x(1:Nblock,1:Nelement),y(1:Nblock,1:Nelement))
call read_filedata(size(y,1),size(y,2),x,y,filename)

do i=1,Nblock
do j=1,Nelement
write(11,*)x(i,j),y(i,j)
enddo
write(11,*)
enddo

stop
end program

総称名を用いる場合の手続き(read_filedata())

module read1
implicit none
interface read_filedata
module procedure &
! dx1 -> double precision array, x(:)
! cy2 -> complex array, y(:,:)
! xyy -> coloum of file, | x y y |
read_dx0_dy1_xy, &
read_dx0_cy1_xyy, &
read_dx1_dy1_xy, &
read_dx1_cy1_xyy, &
read_dx1_dy2_xy, &
read_dx1_cy2_xyy, &
read_dx2_dy2_xy, &
read_dx2_cy2_xyy
end interface read_filedata
contains
subroutine read_dx0_dy1_xy(Ne,y,place,col12)
integer,intent(in)::Ne,col12
character(*),intent(in)::place
double precision,intent(out)::y(1:Ne)

integer::i,ier
double precision::a,b
character(len_trim(place))::fn

y=0d0
fn=trim(place)

open(28,file=fn,status='old',iostat=ier,err=977)
if(col12.eq.1)then
do i=1,Ne
read(28,*)a,b
y(i)=a
enddo
elseif(col12.eq.2)then
do i=1,Ne
read(28,*)a,b
y(i)=b
enddo
else
go to 977
endif

close(28)
return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
write(6,*)"======program stop at read_dx0_dy1_11"
stop
end subroutine read_dx0_dy1_xy

subroutine read_dx0_cy1_xyy(Ne,y,place)
integer,intent(in)::Ne
character(*),intent(in)::place
complex(kind(0d0)),intent(out)::y(1:Ne)

integer::i,ier
double precision::a,b,c
character(len_trim(place))::fn

y=0d0
fn=trim(place)

open(28,file=fn,status='old',iostat=ier,err=977)
do i=1,Ne
read(28,*)a,b,c
y(i)=dcmplx(b,c)
enddo
close(28)
return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
write(6,*)"======program stop at read_dx0_cy1_xyy"
stop
end subroutine read_dx0_cy1_xyy

subroutine read_dx1_dy1_xy(Ne,x,y,place)
integer,intent(in)::Ne
character(*),intent(in)::place
double precision,intent(out)::x(1:Ne),y(1:Ne)

character(len_trim(place))::fn
double precision::a,b
integer::i,ier

x=0d0; y=0d0

fn=trim(place)
open(28,file=trim(fn),status='old',iostat=ier,err=977)
do i=1,Ne
read(28,*)a,b
x(i)=a
y(i)=b
enddo
close(28)
return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
write(6,*)"======program stop at read_dx1_dy1_xy"
stop
end subroutine read_dx1_dy1_xy

subroutine read_dx1_cy1_xyy(Ne,x,y,place)
integer,intent(in)::Ne
character(*),intent(in)::place
double precision,intent(out)::x(1:Ne)
complex(kind(0d0)),intent(out)::y(1:Ne)

integer::i,ier
double precision::a,b,c
character(len_trim(place))::fn

y=0d0
fn=trim(place)

open(28,file=fn,status='old',iostat=ier,err=977)
do i=1,Ne
read(28,*)a,b,c
x(i)=a
y(i)=dcmplx(b,c)
enddo
close(28)
return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
write(6,*)"======program stop at read_dx1_cy1_xyy"
stop
end subroutine read_dx1_cy1_xyy

subroutine read_dx1_dy2_xy(Nb,Ne,x,y,place)
integer,intent(in)::Nb,Ne
double precision,intent(out)::x(1:Ne),y(1:Nb,1:Ne)
character(*),intent(in)::place

integer::i,j,ier
character(len_trim(place))::fn
double precision::a,b

x=0d0; y=0d0

fn=trim(place)
open(28,file=trim(fn),status='old',iostat=ier,err=977)
do i=1,Nb
do j=1,Ne
read(28,*)a,b
if(i.eq.1)x(j)=a
y(i,j)=b
enddo
enddo
close(28)
return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
write(6,*)"======program stop at read_dx1_dy2_xy"
stop
end subroutine read_dx1_dy2_xy

subroutine read_dx1_cy2_xyy(Nb,Ne,x,y,place)
integer,intent(in)::Nb,Ne
double precision,intent(out)::x(1:Ne)
complex(kind(0d0)),intent(out)::y(1:Nb,1:Ne)
character(*),intent(in)::place

integer::i,j,ier
character(len_trim(place))::fn
double precision::a,b,c

x=0d0; y=0d0

fn=trim(place)
open(28,file=trim(fn),status='old',iostat=ier,err=977)
do i=1,Nb
do j=1,Ne
read(28,*)a,b,c
if(i.eq.1)x(j)=a
y(i,j)=dcmplx(b,c)
enddo
enddo
close(28)
return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
write(6,*)"======program stop at read_dx1_cy2_12"
stop
end subroutine read_dx1_cy2_xyy

subroutine read_dx2_dy2_xy(Nb,Ne,x,y,place)
integer,intent(in)::Nb,Ne
double precision,intent(out)::x(1:Nb,1:Ne),y(1:Nb,1:Ne)
character(*),intent(in)::place

integer::i,j,ier
character(len_trim(place))::fn
double precision::a,b

x=0d0; y=0d0

fn=trim(place)
open(28,file=trim(fn),status='old',iostat=ier,err=977)
do i=1,Nb
do j=1,Ne
read(28,*)a,b
x(i,j)=a
y(i,j)=b
enddo
enddo
close(28)
return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
write(6,*)"======program stop at read_dx2_dy2_11"
stop
end subroutine read_dx2_dy2_xy

subroutine read_dx2_cy2_xyy(Nb,Ne,x,y,place)
integer,intent(in)::Nb,Ne
double precision,intent(out)::x(1:Nb,1:Ne)
complex(kind(0d0))::y(1:Nb,1:Ne)
character(*),intent(in)::place

integer::i,j,ier
character(len_trim(place))::fn
double precision::a,b,c

x=0d0; y=dcmplx(0d0,0d0)

fn=trim(place)
open(28,file=trim(fn),status='old',iostat=ier,err=977)
do i=1,Nb
do j=1,Ne
read(28,*)a,b,c
x(i,j)=a
y(i,j)=dcmplx(b,c)
enddo
enddo
close(28)
return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
write(6,*)"======program stop at read_dx2_cy2_12"
stop
end subroutine read_dx2_cy2_xyy
end module read1

プログラミングと数値計算

fortranコンパイラのコマンド

2015年3月7日 sikino コメントする

ここではgfortranとifortでのコンパイル時に指定できるコマンド(最適化とデバッグ、mklの使用とopenmp)について記述します。

ここで対象とするファイル名をmain.f90にします。

最適化

gfortran
gfortran -O3 -fopenmp main.f90
ifort
ifort -mkl -openmp -xHOST -ipo -O3 main.f90

としてコンパイルするといいかと。並列計算がいらない場合は-fopenmp, -openmpのオプションを消してください。

デバッグ

デバッグとはバグを取り除いて望み通りのものにする作業です。

gfortran
gfortran -Wall -fbounds-check -O -Wuninitialized -ffpe-trap=invalid,zero,overflow -fbacktrace -g main.f90
ifort
ifort -mkl -check all -warn all -gen_interfaces -fpe0 -ftrapuv -traceback -g main.f90

とするのがいいでしょう。

※gfortranでは、オプション -ffpe-trap=invalid,zero,overflowが悪さすることがあるようです。fortranでのエラーメモにも書いたように、エラーを起こすはずがないコードでもなぜかエラーと言われることがあります。その時は

gfortran -Wall -fbounds-check -O -Wuninitialized -fbacktrace -g main.f90

でデバッグをやると良いでしょう。

[adsense1]

gfortranでMKLを使う

例えばgfortranコンパイラでIntel® Math Kernel Library (Intel® MKL) を並列計算も含めて使いたいとします。

今、こういう状況にあるとします。

・MKLがあり、/opt/intel/の中にあるディレクトリ”mkl”
が存在する状況下で
・gfortranコンパイラ
がある。

とします。MKLのバージョンによってコンパイルオプションが異なるので注意してください。
Intel® Math Kernel Library Link Line Advisorを見ればどういうコマンドを打てばいいのかを教えてくれます。

ここではMKLのバージョンが
Ver. 11.1
と
Ver. 10.2.5.035
の場合の具体例を載せます。

MKLver11.1の場合

gfortran -fdefault-integer-8 -fopenmp -m64 -I/opt/intel/mkl/include -Wl,--start-group /opt/intel/mkl/lib/intel64/libmkl_gf_ilp64.a /opt/intel/mkl/lib/intel64/libmkl_core.a /opt/intel/mkl/lib/intel64/libmkl_gnu_thread.a -Wl,--end-group -ldl -lpthread -lm main.f90

MKLver10.2.5.035の場合

gfortran -fopenmp -m64 -I/opt/intel/mkl/include /opt/intel/mkl/lib/em64t/libmkl_solver_lp64.a -Wl,--start-group /opt/intel/mkl/lib/em64t/libmkl_gf_lp64.a /opt/intel/mkl/lib/em64t/libmkl_core.a /opt/intel/mkl/lib/em64t/libmkl_gnu_thread.a -Wl,--end-group -lpthread -lm main.f90

※もしかしたら/opt/intel/mkl/ではなくて、/opt/intel/mkl/10.2.5.035/じゃないとだめかもしれません。

バージョンが違ったらどうするかはIntel® Math Kernel Library Link Line Advisorを参考にしてください。上のものを得るには下のようにすれば大丈夫です。

Select Intel® product: (MKLのバージョンは？)	Intel(R) MKL10.2
Select OS: (OSは何？)	Linux*
Select usage model of Intel® Xeon Phi™ Coprocessor: (走らせようとしているコンピュータのcpuの種類がXeon Phi っていうcpuが数百個あるような特別なもの？)	–
Select compiler: (何のコンパイラでMKLを使いたいの？, gfortranだったらGNU fortran.)	GNU Fortran
Select architecture:(良く分かりません！Help me →I don’t use any flag と進んで、”getconf LONG_BIT” と端末で打って何て表示されたかで入力しました。多分cpuが何ビットか、かな？)	Intel(R) 64
Select dynamic or static linking: (動的リンク(dynamic)がいい？静的リンク(static)がいい？　　　　動的→汎用性あり、容量は軽い、動作は遅い　　　　静的→汎用性なし、容量は重い、動作は早い )	Static
Select interface layer:(プログラムの中に整数型で2^31越えてるものはある？ない？越えてなかったら32bit integer, 越えてたら64bit integer)	LP64 (32-bit integer)
Select sequential or multi-threaded layer: (mklのサブルーチンを1つのCPUだけ(sequential)で走らせたい？それとも並列(multi-threaded)にする？)	multi-threaded
Select OpenMP library:(openMPの種類は何使う？intel か GNUか→僕は両方試したら、intelでは動かず、GNUでは動きました。)	GNU (libgomp)
Select cluster library:　ここから下は何を言っているのか分かりませんでした。空白でもいけたから多分直接は関係ないオプションなのだろうと。

次に下の方で作られたオプションをgfortranコンパイル時にくっつける。
$(MKLROOT)/libem64t/libmkl_solver_lp64.a -Wl,–start-group $(MKLROOT)/libem64t/libmkl_gf_lp64.a $(MKLROOT)/libem64t/libmkl_core.a $(MKLROOT)/libem64t/libmkl_gnu_thread.a -Wl,–end-group -lpthread -lm
と
-fopenmp -m64 -I$(MKLROOT)/include
が作られるかと思います。この２つをくっつければokです。
$(MKLROOT)はディレクトリ”mkl”へのパスです。シェルスクリプトの場合、()は要りません。

ここで、この通りコンパイルしようとすると/libem64t/というディレクトリはないよ、っていわれて出来ませんでした。
実際にたどってみると/lib/em64t/ならあって、これに変えたらコンパイルが成功しました。ケースバイケースかな。

参考

Fortranデバッグ用オプション

[adsense2]

プログラミングと数値計算, 物理学

二重振り子

2015年2月20日 sikino コメントする

座標の取り方は下図のように取ります。棒の伸び縮みは無いものとします。

どういう解き方でもいいですが、ここでは

デカルト座標$L(x,y)$でラグランジアンを記述
デカルト座標から座標変換し、$(r,\theta)$でラグランジアンを記述
新たな座標系で運動方程式を導く

という順で解いていきます。

[adsense1]

1, デカルト座標でのラグランジアンLは(運動エネルギーK)-(位置エネルギーU)と書けるため、
$
L(x_1,\dot{x}_1,y_1,\dot{y}_1,x_2,\dot{x}_2,y_2,\dot{y}_2)\\
\displaystyle =\frac{1}{2}m_1(\dot{x}_1^2+\dot{y}_1^2)+\frac{1}{2}m_1(\dot{x}_1^2+\dot{y}_2^2)-(-m_1gy_1-m_2gy_2)
$
と書けます。

2, デカルト座標から座標変換
式を簡単にするために座標変換を行います。新しい座標$(r_1,\theta_1,r_2,\theta_2)$とデカルト座標$(x_1,y_1,x_2,y_2)$の関係式は
$
\begin{align}
x_1&=r_1\sin{\theta_1}\\
y_1&=-r_1\cos{\theta_1}\\
x_2&=r_1\sin{\theta_1}+r_2\sin{\theta_2}\\
y_2&=-r_1\cos{\theta_1}-r_2\cos{\theta_2}
\end{align}
$
という関係があります。各々を時間で微分すれば、
$
\begin{align}
\dot{x}_1&=\dot{r}_1\sin{\theta_1}+r_1\dot{\theta}_1\cos{\theta_1}\\
\dot{y}_1&=-\dot{r}_1\cos{\theta_1}+r_1\dot{\theta}_1\sin{\theta_1}\\
\dot{x}_2&=\dot{r}_1\sin{\theta_1}+r_1\dot{\theta}_1\cos{\theta_1}+\dot{r}_2\sin{\theta_2}+r_2\dot{\theta}_2\cos{\theta_2}\\
\dot{y}_2&=-\dot{r}_1\cos{\theta_1}+r_1\dot{\theta}_1\sin{\theta_1}-\dot{r}_2\cos{\theta_2}+r_2\dot{\theta}_2\sin{\theta_2}
\end{align}
$
これらをラグランジアン$L(x_1,\dot{x}_1,y_1,\dot{y}_1,x_2,\dot{x}_2,y_2,\dot{y}_2)$に代入します。すると、新たな座標系でのラグランジアン$L(r_1,\dot{r}_1,\theta_1,\dot{\theta}_1,r_2,\dot{r}_2,\theta_2,\dot{\theta}_2)$が得られます。
$
\begin{align}
L(r_1,\dot{r}_1,\theta_1,\dot{\theta}_1,& r_2,\dot{r}_2,\theta_2,\dot{\theta}_2) \\
=&\frac{1}{2}m_1(\dot{r}_1^2+r_1^2\dot{\theta}_1^2)+\frac{1}{2}m_2\left[\dot{r}_1^2+\dot{r}_2^2+r_1^2\dot{\theta}_1^2+r_2^2\dot{\theta}_2^2 \right. \\
&\left.+2(\dot{r}_1 r_2 \dot{\theta}_2-r_1\dot{r}_2 \dot{\theta}_1)\sin{(\theta_1-\theta_2)}+2(\dot{r}_1 \dot{r}_2 +r_1 r_2 \dot{\theta}_1 \dot{\theta}_2)+\cos{(\theta_1-\theta_2)}\right] \\
&+m_1gr_1\cos{\theta_1}+m_2g(r_1\cos{\theta_1}+r_2\cos{\theta_2})
\end{align}
$

僕は先ほど式を簡単にするために座標変換をする、といいました。しかし、新しい座標系におけるラグランジアンはどう見ても元のデカルト座標系のラグランジアンに比べて複雑です。この理由は物理的な意味から来ています。
振り子をつないでいる棒が伸び縮みしないとすると系の自由度は角度$\theta_1,\theta_2$の２つです。
となると運動方程式は最高で2本の独立した方程式になるはずです。
しかし、デカルト座標の場合うまく自由度を落とすことができず、運動方程式は4つになってしまいます。
そこで棒が伸び縮みを簡単に表すことができる座標系に移ることでうまく方程式の数を減らせます。

新しい座標系でのラグランジアンで棒の伸び縮みがないという条件を表すには
$
\begin{align}
r_1&=l_1\ (l_1\mbox{は定数}) \\
r_2&=l_2\ (l_2\mbox{は定数})
\end{align}
$
と書けるわけで、また、
$
\begin{align}
\dot{r}_1&=0 \\
\dot{r}_2&=0
\end{align}
$
となるわけです。

$m_1=m_2=m,\ l_1=l_2=l$という場合を特に考えると、ラグランジアンは
$
\displaystyle L(\theta_1,\dot{\theta}_1,\theta_2,\dot{\theta}_2)=ml^2\left[\dot{\theta}_1^2+\frac{1}{2}\dot{\theta}_2^2+\dot{\theta}_1\dot{\theta}_2\cos{(\theta_1-\theta_2)}\right]+mgl(2\cos{\theta_1}+\theta_2)
$
と書けます。あとはラグランジュの運動方程式を当てはめて計算します。

3, 新たな座標系で運動方程式を導く
保存力場中でのラグランジュの運動方程式は
$
\begin{eqnarray}
\left\{
\begin{aligned}
\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{\theta}_1}\right)-\frac{\partial L}{\partial\theta_1}&=0 \\
\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{\theta}_2}\right)-\frac{\partial L}{\partial\theta_2}&=0
\end{aligned}
\right.
\end{eqnarray}
$
なので、代入し、$\ddot{\theta}_1,\ddot{\theta}_2$に関する運動方程式にすれば
$
\begin{eqnarray}
\left\{
\begin{aligned}
\ddot{\theta}_1&=\frac{1}{2}\left\{-\ddot{\theta}_2\cos{(\theta_1-\theta_2)}-\dot{\theta}_2^2\sin{(\theta_1-\theta_2)}-2\frac{g}{l}\sin{\theta_1} \right\} \\
\ddot{\theta}_2&=-\ddot{\theta}_1\cos{(\theta_1-\theta_2)}+\dot{\theta}_1^2\sin{(\theta_1-\theta_2)}-\frac{g}{l}\sin{\theta_2}
\end{aligned}
\right.
\end{eqnarray}
$
となります。これは非線形の２階連立微分方程式です。カオスです。解けません。
数値計算で解きます。用いるのは4次ルンゲ・クッタ法です。

実際に解いてみると、こうなります。

２重振り子の実験もされています。

[adsense2]

カオスが現れる場合、初期値鋭敏性という性質があります。
初期値がほんのちょっと変わるだけでその後の時間発展の様子が大きく変わる性質です。
ちょっと値を変えると上の動画とはまるで違う動きになります。
この動画は上の動画よりも中心に近いほうの粒子の初期速度が4%違います

機械精度の違いでさえも十分な時間発展の後には大きな違いが現れます。複雑であると同時に面白い現象です。

プログラミングと数値計算, 物理学

モンテカルロ法

2015年2月11日 amesyabody 2件のコメント

モンテカルロ法は，乱数を使って系の平衡状態にせまる一つの方法です。

ここでは二次元正方格子のイジング模型について紹介します。
イジング模型は隣り合う格子のスピン同士の相関エネルギーのみを考慮した模型です。
ハミルトニアンは以下のようになります。

$
\displaystyle
H=-J \Sigma_{< i , j >}{s_{i}s_{j}}
$
　
あるスピン状態を持つ系の状態Aについて，そのハミルトニアンを$H_A$とすると，
ボルツマン因子と分配関数$Z$を用いて状態の実現確率を計算することが出来ます。

$
P(A)=\frac{\exp{(-\frac{H_A}{k_B T})}}{Z}
$
　
Aからランダムに一つの格子におけるスピンを一つだけフリップした状態を状態Bとして，
状態Aと状態Bの実現確率比を乱数と比較します。
$P(B)/P(A)>{\rm random}$ならば状態Bを採用し，そうでなければ状態Aを採用します。
ここまでの流れを1mcs(モンテカルロステップ)といい，これを何度も繰り返すことで平衡状態に近づきます。
以下の二つのgifは，相転移温度より高い温度と低い温度で，
モンテカルロステップにより二次元正方格子のスピンがどのように振る舞うかをみているアニメーションです。
黒がダウンスピン，白がアップスピンです。

こちらは温度T=10eVで，自発磁化を持たない事が分かります。

こちらは温度T=0.001eVで，スピンが揃って自発磁化を持つ事が分かります。

画像をクリックするとmcsスタート!!

スピンのふるまいを追うと，自発磁化を持つ温度と持たない温度でスピンのふるまいが違う事が分かります。
二次元のイジング模型はオンサーガーにより厳密に解かれている(相転移温度T=2.26eV at J=1)ので，
それを踏まえて下のプログラムで遊んでみるといいですよ。

これのプログラムソースコード(保証は一切しません!!)
*注意　このプログラムを実行する前にising_anというディレクトリを作ってください。
*追記　2015 2/13　磁化とエネルギーを計算出来るようにしました。
計算方法は重み付き選択法というのを使ってます。

$
\langle A \rangle = \frac{1}{T} \sum_{t=t_0+1}^{t_0 + T}A_{\alpha_{t}}
$
　
十分平衡状態になっている$t_0$以降のモンテカルロステップで，物理量の統計平均を取る方法です。
平衡以前の物理量を計算しないmcsをからまわしと言います。

格子点51×51で，$J=2$で計算しました。
スピンの大きさを1として，相互作用をダブルカウントしているので注意しましょう。
即ち$J=1$の結果が欲しければプログラム上ではintJ=0.125にすればよいかも。

**********parameter**********
i_max::x方向の格子点の数
j_max::y方向の格子点の数
mcs_max::最大モンテカルロステップ数
mcs2step::mcs2step以降のmcsを平衡状態と見なし，物理量を計算します
intj::$J/4$
temp::温度
mov_max::mcs_max/movmax毎のスピン状態をgifアニメーションで見れるように区切ります。
**********parameter**********

単位格子あたりのエネルギーの温度依存性グラフを添付します。

相転移温度付近でエネルギーの変化が大きい事が分かります。

単位格子あたりの磁化の絶対値の温度依存性グラフを添付します。

相転移温度付近で自発磁化が発生し始めている事が分かります。

gnuplotを用いて，格子点上のスピンがどのように平衡状態へ近づくかをダイナミックに確認する事が出来ます。

gnuplot anime.gnuplot
　
を実行するとアニメが始まります。

gnuplot animegif.gnuplot
　
でgifアニメが作成出来ます。

!————————————————
! title = ising2d.f90
! developer = Amesyabody
! released = 2015 2/13
! programming language = fortran
!
!explanation
!2 dimentional square lattice Ising model's dynamical magnetism
!calculated by Monte Carlo method.
!————————————————

program main
implicit none
integer(4)::i,j,i_max,j_max,mcs,mcs_max,mov,mov_max,mcs2step
parameter(i_max=50,j_max=50,mcs_max=1000000,&
mcs2step=500000,mov_max=100)
real(8)::intj,temp
parameter(intj=0.25d0,temp=0.1d0)
real(8)::spin(0:i_max,0:j_max)
real(8)::ratio,ham1,ham2,hamres,dum,ene
integer(4)::px,mx,py,my,rndx,rndy
character(10)::movc,movc2

real(8)::mag,mag2,magres

!-----make random-----
integer(4),allocatable::seed(:)
integer cnt,seedsize
real(8)::rnd,rndi,rndj,rnds

call random_seed(size=seedsize)
allocate(seed(seedsize))
write(*,*)seedsize

dum=rndmf(0)

call system_clock(count=cnt)
seed=cnt!+idnint(dum*1000000)
! dum=rndmf(seed(0))
! seed=idnint(seed*dum)
call random_seed(put=seed)

!------------------

mag2=0.0d0
ene=0.0d0

do mov=0,mov_max

write(movc,'(i3)')mov

open(100+mov,file="ising_an/spin"//trim(adjustl(movc))//".dat"&
,status="unknown")

end do

open(5000,file="ising_an/anime.gnuplot",status="unknown")
open(5001,file="ising_an/animegif.gnuplot",status="unknown")

!-----primitive spin set-----

do i=0,i_max
do j=0,j_max
call random_number(rnd)
! write(*,*)rnd
if(rnd>0.5d0)then
spin(i,j)=1.0d0
else
spin(i,j)=-1.0d0
end if
! write(*,*)spin(i,j)!,cnt,seed(1)

write(100,'(i3,i3,f5.1)')i,j,spin(i,j)

end do

write(100,*)" "

end do

!----------------------------

do mcs=1,mcs_max

do i=0,i_max

do j=0,j_max

px=i+1
if(i==i_max)px=px-i_max-1
py=j+1
if(j==j_max)py=py-j_max-1
mx=i-1
if(i==0)mx=mx+i_max+1
my=j-1
if(j==0)my=my+j_max+1

ham1=ham1-intj*spin(i,j)*spin(px,j)&
-intj*spin(i,j)*spin(i,py)&
-intj*spin(mx,j)*spin(i,j)&
-intj*spin(i,my)*spin(i,j)

do mov=1,mov_max-1

if(mcs==mov*mcs_max/mov_max)&
write(100+mov,'(i3,i3,f5.1)')i,j,spin(i,j)

end do

if(mcs>mcs2step)then

mag2=mag2+spin(i,j)

end if

end do

do mov=1,mov_max-1

if(mcs==mov*mcs_max/mov_max)write(100+mov,*)" "

end do

end do

if(mcs>mcs2step)then

magres=magres+abs(mag2)
mag2=0.0d0
ene=ene+ham1

end if

call random_number(rndi)
call random_number(rndj)
rndx=idnint(rndi*i_max)
rndy=idnint(rndj*j_max)

! write(*,*)rnd
!write(*,*)spin(rndx,rndy)
spin(rndx,rndy)=-spin(rndx,rndy)
!write(*,*)spin(rndx,rndy)
do i=0,i_max

do j=0,j_max

px=i+1
if(i==i_max)px=px-i_max-1
py=j+1
if(j==j_max)py=py-j_max-1

mx=i-1
if(i==0)mx=mx+i_max+1
my=j-1
if(j==0)my=my+j_max+1

ham2=ham2-intj*spin(i,j)*spin(px,j)&
-intj*spin(i,j)*spin(i,py)&
-intj*spin(mx,j)*spin(i,j)&
-intj*spin(i,my)*spin(i,j)

end do

end do

ratio=exp((-ham2+ham1)/temp)
! write(*,*)ham1,ham2,ratio

ham1=0.0d0
ham2=0.0d0

call random_number(rnds)

if(ham1 > ham2)then!metroporis
goto 200
end if

if(rnds>ratio)then
spin(rndx,rndy)=-spin(rndx,rndy)
end if

200 continue

end do

do i=0,i_max
do j=0,j_max

mag=mag+spin(i,j)

px=i+1
if(i==i_max)px=px-i_max-1
py=j+1
if(j==j_max)py=py-j_max-1
mx=i-1
if(i==0)mx=mx+i_max+1
my=j-1
if(j==0)my=my+j_max+1

hamres=hamres-intj*spin(i,j)*spin(px,j)&
-intj*spin(i,j)*spin(i,py)&
-intj*spin(mx,j)*spin(i,j)&
-intj*spin(i,my)*spin(i,j)

write(100+mov_max,'(i3,i3,f5.1)')i,j,spin(i,j)

end do

write(100+mov_max,*)" "

end do

mag=mag/(dble(i_max*j_max))

! write(*,'(3f10.5,f15.5)')temp,intj,mag,hamres

ene=ene/(dble((mcs_max-mcs2step)*(i_max+1)*(j_max+1)))
magres=magres/(dble((mcs_max-mcs2step)*(i_max+1)*(j_max+1)))

write(*,'(A,f10.5,A,f10.5)')"energy=",ene," magnet=",magres

write(5000,*)"set size sq"
write(5000,*)"set pm3d map"
write(5000,*)"set palette rgbformulae 21,22,23"
write(5000,*)"set xlabel 'x'"
write(5000,*)"set ylabel 'y'"
write(5000,*)"set zra[-1:1]"
write(5000,*)"set zlabel 'spin'"
write(5000,*)"splot 'spin0.dat'"
write(5000,*)"pause 0.1"
do mov=1,mov_max
write(movc,'(i3)')mov
write(5000,*)"splot 'spin"//trim(adjustl(movc))//".dat'"
write(5000,*)"pause 0.1"

end do

write(5001,*)"set tics font 'Times,15'"
write(5001,*)"set xlabel font 'Times,15'"
write(5001,*)"set ylabel font 'Times,15'"
write(5001,*)"set size sq"
write(5001,*)"set pm3d map"
write(5001,*)"set palette rgbformulae 21,22,23"
write(5001,*)"set xlabel 'x'"
write(5001,*)"set ylabel 'y'"
write(5001,*)"set zra[-1:1]"
write(5001,*)"set cbrange [-1:1]"
write(5001,*)"set zlabel 'spin'"
write(5001,*)"set term gif animate optimize"
write(5001,*)"set output 'is_t10.gif'"
write(5001,*)"do for[i=0:100] {"
write(5001,*)'file = sprintf("spin%01d.dat", i)'
write(movc2,'(i10)')mcs_max/mov_max
write(5001,*)'time = sprintf("s=%d[mcs]",i*'//trim(adjustl(movc2))//')'
write(5001,*)"set title time"
write(5001,*)"splot file"
write(5001,*)"}"

contains

real(8) function rndmf(seeds)
implicit none

integer(4)::a,m,q,p,n,ndiv,j,k,seeds
real(8)::rm,rmax
parameter(a=16807,m=2147483647,rm=1.0/m)
parameter(q=127773,p=2836,n=32,ndiv=1+(m-1)/n)
parameter(rmax=1.0-1.2e-7)
integer(4)r(n),r0,r1

if (seeds .ne. 0)then
r1=abs(seeds)
do j=n+8,1,-1
k=r1/q
r1=a*(r1-k*q)-p*k
if(r1 .lt. 0)r1=r1+m
if(j .le. n)r(j)=r1
end do
r0=r(1)
end if

k=r1/q
r1=a*(r1-k*q) -p*k
if(r1 .lt. 0)r1=r1+m
j=1+r0/ndiv
r0=r(j)
r(j)=r1
rndmf=min(rm*r0,rmax)

end function

end program

プログラミングと数値計算

プリクラ問題

2015年1月31日 sikino コメントする

covering designのt=2に相当、
La Jolla Covering Repository Tablesに全部載ってます…
リンク先クリックすればいいんだね。
Combinatorial covering designs に文献とかあるから、なんかもう・・・
もういいや

fortran90によるのコードはこのページの下の方です
だんだんと条件を付け加えて更新していくつもりです。n=9~10でm=3,4程度は厳しいです。

力技でやってみたよ!
主眼が回数になっているので、実際の組み合わせは？と気になってやったものです。

事の発端

あのね。前々からずーーっと答えが知りたい数学の問題があって。 ①n人のグループがいる ②プリクラ機には一度にm人しか入れないメンバーがほかの各メンバーと一度は同じ写真におさまるためには最低何回撮影する必要があるか？　※ちなみに(n,m)=(7,5)の時は3 #プリクラ問題

— あしやまひろこ (@hiroko_TB) 2015, 1月 1

コンピュータによる総当たりでは今のままではn=10位が限界です。
だんだんと改良して増やせていければ、と。
コードは後ほど僕に時間ができ次第、公開します。

#プリクラ問題もすっかり大人しくなりましたが、懲りずに表（n≦"99", m≦25）をアップします。現人神ジョージ先生の解答を待ちましょう。 pic.twitter.com/5q3p1seCUa

— リッチカスタードソフト (@298itf) 2015, 1月 6

より↓

念のための確認もかねて、計算は表の値-1から確かめて総当たりを行った結果です。なので具体的な組み合わせが載っているものは正しいはずです。
組み合わせの総数はおおよそ
$\displaystyle \sum_{q=3}^{q^{Answer}} {}_{{}_nC_m}C_{q}$
という膨大な数の組み合わせとなります。ここで$q^{Answer}$はプリクラ問題の最小回数です。
$q=3$はあり得ない、全てを検証する必要はない、など確実にわかる要素は多いですが、例えばn=9,m=3,解q=12を考えると
$\displaystyle {}_{{}_9C_3}C_{12}=112,992,892,764,570 \ \ \sim 10^{14}(=100T)$回
の計算が必要です。理想的な状況を考え、cpu10GHz,コア数10で並列処理で行う場合、計算時間は
$\displaystyle \frac{10^{14}回}{10\times 10^{9}Hz\cdot 10\mbox{コア}}=1000[秒]$
なのでぎりぎり計算できそうかなーくらいになります。
しかし、このあたりが限界で、n=10,m=3,解q=17になると、
$\displaystyle {}_{{}_{10}C_3}C_{17}=189,916,591,435,396,829,640\ \ \sim 10^{20}(100E)$回
の計算となります。先ほどの設定で行うならば、$10^{9}$秒~31.7年かかります。
スパコンの京(10PHz)をフルに使えるならば1000[s]で終わります。しかしそこで終わりです。n=11は…もう無理です。

効率のいい方法を考えない限り、解は得られません。

n\m	2	3	4	5	6	7	8	9	10	11
2	1	1	1	1	1	1	1	1	1	1
3	3	1	1	1	1	1	1	1	1	1
4	6	3	1	1	1	1	1	1	1	1
5	10	4	3	1	1	1	1	1	1	1
6	15	6	3	3	1	1	1	1	1	1
7	21	7	5	3	3	1	1	1	1	1
8	28	11	6	4	3	3	1	1	1	1
9	36	12	8	5	3	3	3	1	1	1
10	45	17	9	6	4	3	3	3	1	1
11	55	19	11	7	6	4	3	3	3	1

(3,2)=3

1 2
1 3
2 3

(4,2)=6

略

(4,3)=3

1 2 3
1 2 4
1 3 4

(5,2)=10

略

(5,3)=4

1 2 3
1 2 4
1 2 5
3 4 5

(5,4)=3

1 2 3 4
1 2 3 5
1 2 4 5

(6,2)=15

略

(6,3)=6

1 2 3
1 2 4
1 5 6
2 5 6
3 4 5
3 4 6

(6,4)=3

1 2 3 4
1 2 5 6
3 4 5 6

(6,5)=3

1 2 3 4 5
1 2 3 4 6
1 2 3 5 6

(7,2)=21

略

(7,3)=7

1 2 3
1 4 5
1 6 7
2 4 6
2 5 7
3 4 7
3 5 6

(7,4)=5

1 2 3 4
1 2 3 5
1 2 3 6
1 2 3 7
4 5 6 7
or
1 2 3 5
1 4 6 7
2 3 4 5
2 4 6 7
3 5 6 7

(7,5)=3

1 2 3 4 5
1 2 3 6 7
1 4 5 6 7

(7,6)=3

1 2 3 4 5 6
1 2 3 4 5 7
1 2 3 4 6 7

(8,2)=28

略

(8,3)=11

計算中…

(8,4)=6

1 2 3 4
1 2 5 6
1 2 7 8
3 4 5 6
3 4 7 8
5 6 7 8

(8,5)=4

1 2 3 4 5
1 2 3 4 6
1 2 3 7 8
4 5 6 7 8

(8,6)=3

1 2 3 4 5 6
1 2 3 4 7 8
1 2 5 6 7 8

(8,7)=3

1 2 3 4 5 6 7
1 2 3 4 5 6 8
1 2 3 4 5 7 8

(9,2)=36

略

(9,3)=12

計算中…

(9,4)=8

1 2 3 5
1 2 3 6
1 2 4 7
1 2 8 9
1 3 4 8
1 3 7 9
4 5 6 9
5 6 7 8

info)cputime=3.0415[s], Nc=10523209, com=1 2 35 115 126

(9,5)=5

1 2 3 4 5
1 2 3 4 6
1 2 7 8 9
3 4 7 8 9
5 6 7 8 9

(9,6)=3

1 2 3 4 5 6
1 2 3 7 8 9
4 5 6 7 8 9

(9,7)=3

1 2 3 4 5 6 7
1 2 3 4 5 8 9
1 2 3 6 7 8 9

(9,8)=3

1 2 3 4 5 6 7 8
1 2 3 4 5 6 7 9
1 2 3 4 5 6 8 9

(10,2)=45

略

(10,3)=17

計算厳しい

(10,4)=9

計算中(4日後くらい?)

(10,5)=6

(10,6)=4

1 2 3 4 5 6
1 2 3 4 7 8
1 2 3 4 9 10
5 6 7 8 9 10

(10,7)=3

1 2 3 4 5 6 7
1 2 3 4 8 9 10
1 5 6 7 8 9 10

(10,8)=3

1 2 3 4 5 6 7 8
1 2 3 4 5 6 9 10
1 2 3 4 7 8 9 10

(10,9)=3

1 2 3 4 5 6 7 8 9
1 2 3 4 5 6 7 8 10
1 2 3 4 5 6 7 9 10

(11,2)=55

略

(11,3)=19

(11,4)=11

(11,5)=7

(11,6)=

(11,7)=4

1 2 3 4 5 6 7
1 2 3 4 5 6 8
1 2 3 4 9 10 11
5 6 7 8 9 10 11

(11,8)=3

1 2 3 4 5 6 7 8
1 2 3 4 5 9 10 11
1 2 6 7 8 9 10 11

(11,9)=3

1 2 3 4 5 6 7 8 9
1 2 3 4 5 6 7 10 11
1 2 3 4 5 8 9 10 11

(11,10)=3

1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 11
1 2 3 4 5 6 7 8 10 11

fortran90のコード。総当たりで調べます。
別にinputファイルが必要です。n:m=3:2以上の時は考慮していません。3回で終わるからね。
下の例はn=6, m=3に相当しています。(inputのrはmに相当します。)
qiは調べ始める最低回数、qeは調べ終わる最高回数です。
出力はこうなります。

M 35
q --> 3
q --> 4
q --> 5
q --> 6
q --> 7
number of calculation 1276464
cpu_time 0.905050993 [sec]
+------------+
1 2 3
1 4 5
1 6 7
2 4 6
2 5 7
3 4 7
3 5 6
+------------+

ここからfortran90のコード。gfortranで確かめています。

program main
implicit none
integer::n,r,qi,qe,qc
integer::i,j
real::t1,t2
integer,allocatable::groups(:,:)

n=7
r=3
qi=3
qe=10
allocate(groups(1:qe,1:r))
groups=0

call cpu_time(t1)
call pkprob(n,r,qi,qe,qc,groups)
call cpu_time(t2)
write(6,*)"cpu_time",t2-t1,"[sec]"

write(6,'(A)')"+------------+"
do i=1,qc
do j=1,r
write(6,'(i3)',advance='no')groups(i,j)
enddo
write(6,*)
enddo
write(6,'(A)')"+------------+"

stop
end program main
!==========================
subroutine pkprob(n,r,qi,qe,qc,groups)
!Covering Design t=2
implicit none
integer,intent(in)::n,r,qi,qe
integer,intent(out)::qc,groups(1:qe,1:r)

integer::q,M,ic,ncalc
integer::i,j,q2,k
integer,allocatable::ini(:),tmp(:),cgroup(:,:),pk(:,:)
integer,allocatable::choose(:)

integer::check1,nCr
integer::csign
! logical::csign,gsign
external::nCr

!Combination M=nCr
M=nCr(n,r)
write(6,*)"M",M

!csign : Sign of critical combination of array
!ic : Critical value of i which number of combination of combination
!qc : Critical value of q which number of combination of combination
!pk(1:n,1:n): Checking array if every number fill(1) combination or not(0)
!ncount(1:n): To check calculation number of times
allocate(ini(1:r),tmp(1:r),cgroup(1:M,1:r))
ini(1:r)=0
cgroup(1:M,1:r)=0
tmp(1:r)=0

do j=1,r
ini(j)=j
tmp(j)=ini(j)
enddo
!Substitute combination to cgroup.
do i=1,M
cgroup(i,1:r)=tmp(1:r)
call combination(n,r,ini,tmp)
enddo
deallocate(ini,tmp)
!csign : Sign of critical combination of array
!ic : Critical value of i which number of combination of combination
!qc : Critical value of q which number of combination of combination
!pk(1:n,1:n): Checking array if every number fill(1) combination or not(0)
!ncount(1:n): To check calculation number of times
csign=0
ic=0
qc=0
allocate(pk(1:n,1:n))
pk=0

ncalc=0

!q : We need to consider below combination.
! Example ) n=5, r=3 case.
!
!cgroup(i,1:r)
!i\r| 1 2 3
!----------
! 1 | 1 2 3
! 2 | 1 2 4
! 3 | 1 2 5
! 4 | 1 3 4
! 5 | 1 3 5
! 6 | 1 4 5
! 7 | 2 3 4
! 8 | 2 3 5
! 9 | 2 4 5
!10 | 3 4 5
!
!We have to choose q= 3 or 4 or 5 ...combination of "cgroup".
!Sum total of these combination can write D=_{nCr}C_q.
!Specifically,
!---choose q groups from cgroup
! +----q=3----+ +-----q=4-----+ +------q=5------+
! |choose(1:q)| | choose(1:q) | | choose(1:q) |
! | l | | | l | | | l | |
! |===|=======| |===|=========| |===|===========|
! | 1 | 1 2 3 | | 1 | 1 2 3 4 | | 1 | 1 2 3 4 5 |
! | 2 | 1 2 4 | | 2 | 1 2 3 5 | | 2 | 1 2 3 4 6 |
! | 3 | 1 2 5 | | 3 | 1 2 3 6 | | 3 | 1 2 3 4 7 |
! | 4 | 1 2 6 | | 4 | 1 2 4 7 | | 4 | 1 2 3 4 8 |
! | 5 | 1 2 7 | | 5 | 1 2 4 8 | | 5 | 1 2 3 4 9 |
! | 6 | 1 2 8 | | 6 | 1 2 5 9 | | 6 | 1 2 3 4 10|
! | 7 | 1 2 9 | | 7 | 1 2 5 10| | 7 | 1 2 3 5 6 |
! |...| . . . | |...| . . . . | |...| . . . . |
! |120| 4 5 6 | |210| 7 8 9 10| |252| 6 7 8 9 10|
! +-----------+ +-------------+ +---------------+
!
! If q=4, l=3
! We choose below combination, and check satisfy condition or not by "pk".
! 1 | 1 2 3
! 2 | 1 2 4
! 3 | 1 2 5
! 6 | 1 4 5

do q=qi,qe
write(6,'(A,i5)')"q -->",q
allocate(choose(1:q),ini(1:q))
choose(1:q)=0
ini(1:q)=0

do j=1,q
ini(j)=j
choose(j)=ini(j)
enddo

csign=0
do while(csign.eq.0)
!condition 1 to reduce calculation
check1=cgroup(choose(1),1)
if(check1.ne.1)exit

!condition 2 to reduce calculation
check1=cgroup(choose(2),1)
if(check1.ge.2)exit

!condition 3 to reduce calculation for less than n:r=3:2
check1=cgroup(choose(q),1)
if(check1.ge.3)then
pk=0
do q2=1,q
do j=1,r
do k=1,r
pk(cgroup(choose(q2),j),cgroup(choose(q2),k))=1
enddo
enddo
enddo
!memory calculation number of times
ncalc=ncalc+1
if(minval(pk).eq.1)then
csign=1
qc=q
endif

if(csign.eq.1)exit
endif

!go forward choose(1:q)
call combination(M,q,ini,choose)
enddo

if(csign.eq.1)exit
deallocate(choose,ini)
enddo
write(6,*)"number of calculation",ncalc

groups(1:qe,1:r)=0
do i=1,qc
do j=1,r
groups(i,j)=cgroup(choose(i),j)
enddo
enddo

return
end subroutine pkprob
!---------------------------------
subroutine combination(n,r,ini,arr)
implicit none
integer(8)::i,n,r,ini(1:r),bef(1:r),arr(1:r)
integer(8)::numx
logical::key(1:r)

bef(1:r)=arr(1:r)-ini(1:r)
arr(1:r)=0
key(1:r)=.false.

numx=n-r
do i=1,r
if(bef(i).eq.numx)key(i)=.true.
enddo

do i=1,r-1
if(key(i+1))then
if(key(i))then
if(i.ne.1)arr(i)=arr(i-1)
else
arr(i)=bef(i)+1
endif
else
arr(i)=bef(i)
endif
enddo
if(key(r))then
arr(r)=arr(r-1)
else
arr(r)=bef(r)+1
endif

arr(1:r)=arr(1:r)+ini(1:r)

return
end subroutine combination
!----------------------------------
function nCr(n,r)
implicit none
integer(8)::n,r,i,r0,nCr

r0=n-r
if(r.le.r0)r0=r

nCr=1
do i=1,r0
nCr=nCr*(n-r0+i)
nCr=nCr/i
enddo

return
end function nCr

テキストエディタ, プログラミングと数値計算

フォントサイズ変更

2015年1月31日 sikino コメントする

emacsでフォントサイズを変更したい場合は？環境はemacs24。

(調べてもよくわからないし出てこない…。)
一度だけ変更したい場合は、Ctrl-x-± で可能です。
ではずっと変更したい場合は？
emacsの設定ファイルに記述します。emacsの設定ファイルの場所は
emacswiki – InitFileより

For GnuEmacs, it is ~/.emacs or _emacs or ~/.emacs.d/init.el.

For XEmacs, it is ~/.xemacs or ~/.xemacs/init.el.

For AquamacsEmacs, it is ~/.emacs or ~/Library/Preferences/Aquamacs Emacs/Preferences.el

にあります。環境によって場所が異なります。
Linuxmint17.1,”sudo apt-get install emacs24″で行った場合、おそらく~/.emacsでしょう。

21.8 Fontsを見るとデフォルトでは
フォント(の種類) : monospace font
フォントサイズ : 10pt
として設定されているそうです。

テキストエディタで~/.emacsを開き、
(custom-set-variables
…
…
)
の括弧の外に

(add-to-list 'default-frame-alist
'(font . "DejaVu Sans Mono-12"))

という文を記述します。
現在、フォントのタイプの変更は僕は分かりません。多分DejaVu…を変えればよさそう。
12のところを10とか14とかに変えればデフォルトのフォントサイズを変更することができます。

見出しに使っている画像は8,9,10,11,12,20の場合の比較です。参考にどうぞ。

僕の使ってる~./emacsはこちら