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

反復計算における収束判定について

2020年4月20日 sikino コメントする

収束判定

漸化式
\begin{equation}
x^{(k+1)}=x^{(k)}+\Delta^{(k)}
\end{equation}
に従って、$\lim_{k\to\infty}x^{(k)}=a$が適当な定数に収束する問題を考えます。
絶対誤差$\varepsilon_A$、相対誤差$\varepsilon_R$とすると、数値計算を終了する時は、以下の不等式が満たされる時にすると良いです。
\begin{equation}
|x^{(k+1)}-x^{(k)}|\lt \varepsilon_A + \varepsilon_R(|x^{(k)}|+|x^{(k+1)}|)
\end{equation}

判定式の意味

上記判定式の意味を考えましょう。上の判定式は絶対誤差による評価と相対誤差による評価をまとめて評価する式になっています。

相対誤差が重要な場合

本当の解が$a=100$であるとします。その時、$x^{(k)}=100.3, x^{(k+1)}=100.2$であったと仮定すると、大雑把に計算して
\begin{align}
|x^{(k+1)}-x^{(k)}| &\lt \varepsilon_A + \varepsilon_R(|x^{(k)}|+|x^{(k+1)}|) \\
0.1 &\lt \varepsilon_A + \varepsilon_R\times 100
\end{align}
となります。右辺は本来2倍ですが、2倍はそんなに重要ではないので無視しています。もし仮に、$\varepsilon\equiv\varepsilon_A=\varepsilon_R$としていれば、
\begin{align}
0.1 &\lt 100\varepsilon_R
\end{align}
と近似できます。以上から、本当の解が1より十分大きい場合、相対誤差で収束を評価する式になっていると言えます。
つまり上記の問題設定の場合、$\varepsilon_R\gt 10^{-3}$にしてあれば”収束した”として判定します。

絶対誤差が重要な場合

本当の解が$a=0.01000$であるとします。その時、$x^{(k)}=0.01002, x^{(k+1)}=0.01001$であったと仮定すると、大雑把に計算して
\begin{align}
|x^{(k+1)}-x^{(k)}| &\lt \varepsilon_A + \varepsilon_R(|x^{(k)}|+|x^{(k+1)}|) \\
0.00001 &\lt \varepsilon_A + \varepsilon_R\times 0.01
\end{align}
となります。右辺は本来2倍ですが、2倍はそんなに重要ではないので無視しています。もし仮に、$\varepsilon\equiv\varepsilon_A=\varepsilon_R$としていれば、
\begin{align}
0.00001 &\lt 0.01\varepsilon_A
\end{align}
と近似できます。以上から、本当の解が1より十分小さい場合、絶対誤差で収束を評価する式になっていると言えます。
つまり上記の問題設定の場合、$\varepsilon_A\gt 10^{-3}$にしてあれば”収束した”として判定します。

参考

[1]杉山耕一朗, OBOROの収束判定条件の設定

プログラミングと数値計算, 物理学

動かない壁に対する束縛運動と反射

2019年10月20日 sikino コメントする

動かない壁に対する束縛運動と反射を考えます。
例えば、初め跳ねてた運動が、壁に沿って動く運動に変化する、という状況です。
あんまり見たことが無いので、面白そうだと思いました。

束縛されている時と反発するときは、式(1),(2)によってあらわされます。それらは
壁に束縛されている場合

壁と反発する場合

です。壁と反発する場合、反発後の速度は式(3)に沿って動きます。

※式(1), (2)では壁が時間依存しても良いように定式化しています。この定式化は恐らく正しいです。また、本稿に載せているプログラムも壁が時間依存しても良いように作成していますが、動く壁の場合、プログラムではうまく計算が出来ません。

ここで、$e$は反発係数、$\mathbf{n}$は壁の法線ベクトルであり、

$\nabla$はナブラ演算子であり

と与えられます。また、$\hat{H}f$は関数$f$のヘッセ行列であり,

と与えられます。

定式化や数値計算手法の詳細は以下のページを参照してください。

壁との反発と束縛運動の定式化

質点と壁との反発を表す運動方程式
 束縛条件下の運動 – 自由度がうまく落とせない運動

数値計算手法

ルンゲ=クッタ法の説明と刻み幅制御
 Hyper-dual numbersによる二階偏微分の計算
 ゼロ点を探す（二分法、挟み撃ち法、Anderson-Björk法、Brent法、Newton法、Steffensen法）

解放・束縛判定

ここでいう”解放”とは、束縛されていなく、反射を繰り返している状態を表します。また、”束縛”は壁に沿って運動している状態です。

前提として、壁を通り抜けることは無いと考えます。

すなわち、時刻$t=t_0$で位置$\mathbf{r}=\mathbf{r}_0$の時、任意の時刻$t$について

が満たされるとします。
もし、$f(\mathbf{r}_0,t_0) =0$ならば判断がつかないため、プラスかマイナスはこちらから与えます。

解放→解放判定、解法→束縛判定

解放状態から壁によって単に反発する場合

関数$f(\mathbf{r},t)$の符号が変化した時、壁との反発を考えます。
壁の法線方向の速度が十分に大きい場合、壁と反発し、そうでない場合、壁に束縛されると考えます。
その前提の元、

を満たす$t=t_i, \mathbf{r}=\mathbf{r}_i$を見つけます。
その後、式(6)に従い、速度を変化させます。

数値計算的には、衝突直前の時刻を採用します。すなわち、
ゼロ点を探す際にある範囲$t_a \leq t_i \leq t_b$で挟み込んでいくのですが、$t_b$は壁を越えてしまうので採用しません。

解放状態から壁沿いに束縛される場合

もし、壁の法線方向の速度が十分に小さい場合（ある閾値を下回った場合）、壁に束縛されると考えます。
この時、壁の法線方向の速度はゼロに変更します。
すなわち、速度は時刻$t=t_i$で

を持ちますが、$v_{\parallel}=0$にしてから、束縛運動に移行するということです。
これは、$\mathbf{v}$と壁の法線方向のベクトル$\mathbf{e}_n$の内積を取ることで得られます。
また、束縛された瞬間($t=t_i$)の束縛力$C(t_i)$を計算し、その符号を記録しておきます。
束縛力$C(t)$は、

です。
この符号が変化した瞬間が壁からの束縛が無くなる時なので、そのために記録します。

束縛→解放判定

束縛力$C(t)$の符号と$C(t_i)$の符号が変わるまで、式(1)に従い、時間発展させます。
すなわち、

を満たす$t=t’, \mathbf{r}=\mathbf{r}’$を見つけます。
その後、式(2)に従い運動します。
式(2)の運動では束縛力は働かないので、符号は自然と初期条件の符号と同一になる（※）。

※この条件はあまり良くありません。この判別方法のせいで、壁が時間依存している場合、束縛力が働いていない一瞬で質点が壁を超えてしまいます。プログラム自体は壁は時間依存しても良いことになっていますが、この条件分岐は上手く動きません。以下に示す計算プログラムは、質点が束縛されている場合に壁が時間依存しなければ正しいです。

プログラム

プログラムは以下のリンク先においておきます。

https://slpr.sakura.ne.jp/qp/supplement_data/lag_ver1.0.tar.gz

適当なディレクトリで展開し、lag_ver1.0というディレクトリに移動してから以下のコマンドで実行できます。

$ sh cm.sh
$ ./a.out
&INPUT
MASS= 1.0000000000000000 ,
G= 1.0000000000000000 ,
TA= 0.0000000000000000 ,
TB= 20.000000000000000 ,
NT= 101,
ELS= 0.59999999999999998 ,
RX0= -1.0000000000000000 ,
RY0= 1.0000000000000000 ,
VX0= 3.0000000000000000 ,
VY0= 0.0000000000000000 ,
RKTOL= 1.0000000000000000E-008,
ZRTOL= 1.0000000000000000E-008,
TRTOL= 1.0000000000000000E-008,
REGION= 1,
/
0
$ 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> load "anime.plt"
gnuplot>

動かすと以下のような動画が描画されます。

デフォルトでは
ポテンシャル$V=mgy$（サブルーチンfp2d）
壁$f(x,y)=x^2+y^2-4$（サブルーチンfw2d）
に書かれています。
摩擦、空気抵抗は入っていません。唯一、反発係数(els)だけがinputファイルの中に書かれています。

初期条件の

rx0 = -2d0, ! Initial condition
ry0 = 2d0, ! position and velocity
vx0 = 1d0, ! \mathbf{r} = (rx, ry)
vy0 = 0d0, ! \mathbf{v} = (vx, vy)

だけを上の通り変更すると以下のような振る舞いをします。

確かめ

確かめを行います。
重力$g$の下で、質量$m$の物体が、半径$r$の円形の壁の内側に沿ってボールが進み、壁からの抗力が無くなり、壁から離れる状況を考えます(下の図を参照)。

円形の壁に沿っている時、垂直抗力$N$は

と書けます。エネルギー保存則より、

が成り立っています。ここで、$v_0\equiv v(t=0)$と置きました。
垂直抗力$N$がゼロになる点が壁から離れる条件ですので、$v_0$を用いて

と書けます。初速度が分かっている時、壁から離れる角度は

です。もしくは、壁から離れる角度が分かっている場合、初速度は

と与えられます。重力加速度, 半径を$g=1, r=2$とし、$y=1$と壁との交点、すなわち$\theta=2\pi/3$の場合、初速度は

です。実際に本稿のプログラムを動かしますと、$y=1$でちょうど離れていることが確認できます。

ここで、青線は壁に沿って動いて運動している状態であり、赤線は壁から離れている運動している状態です。

剛体に対する散乱

ポテンシャルを無くし(g=0)、弾性散乱(els=1)を考えると
古典的な散乱問題的なものもできます。

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

Cubic補間

2019年5月25日 sikino コメントする

多項式補間です。

Cubic補間 (1次元) $O(h^4)$
Bicubic補間 (2次元) $O(h^4)$
補間 (2次元) $O(h^2)$
参考文献

Cubic補間 (1次元)

4点を厳密に通る3次多項式によって関数を補間します。
概要は以下の図の通りです。
補間したい点よりも小さいデータ点を2点、大きいデータ点を2点使って補間します。
いわば小さい区間で区切ったラグランジュ補間です。
補間は,既知の係数$a_{i}$を用いて関数
$\displaystyle
g(x)=\sum_{i=0}^{3}a_{i}x^i
$
誤差は$O(h^4)$です。

Fortranプログラムはこちら。

program main
implicit none
integer::i,N,M
double precision::x,f,df,df2,s
double precision,allocatable::xdata(:),fdata(:)
double precision::pi=dacos(-1d0)
double precision,external::cubic

N=30
M=100
allocate(xdata(0:N),fdata(0:N))
xdata=0d0
fdata=0d0

do i=0,N
xdata(i)=dble(i)*0.1d0*pi
fdata(i)=sin(xdata(i))
write(10,*)xdata(i),fdata(i)
enddo

! Cubic-spline interpolation given position as point
do i=0,M
x=dble(i)*0.03d0*pi-1d0
write(11,'(2e20.7e2)')x,cubic(x,N,xdata,fdata)
enddo

stop
end program main

double precision function cubic(x,N,x0,f0)
implicit none
integer,intent(in)::N
double precision,intent(in)::x,f0(0:N),x0(0:N)

integer::i,i0,i1,i2,i3
double precision::tx
double precision::a,b,c,d,p,q,r,s,t,u

tx = x-x0(0)
i1 = 0
do i=1,N-2
tx = x-x0(i)
if(tx.gt.0d0)then
i1 = i
else
exit
endif
enddo
if(i1.eq.0)i1=1
i0=i1-1
i2=i1+1
i3=i1+2

a = x-x0(i0)
b = x-x0(i1)
c = x-x0(i2)
d = x-x0(i3)

p = x0(i1)-x0(i0)
q = x0(i2)-x0(i1)
r = x0(i3)-x0(i2)

s = x0(i2)-x0(i0)
t = x0(i3)-x0(i1)

u = x0(i3)-x0(i0)

cubic = a*c*( d*f0(i1)/(p*q) + b*f0(i3)/(u*r) )/t &
- b*d*( c*f0(i0)/(p*u) + a*f0(i2)/(q*r) )/s

return
end function cubic

Bicubic補間 (2次元)

16点を厳密に通る多項式によって関数を補間します。
補間したい点を囲むように存在する16点を用いて補間します。
補間は,既知の係数$a_{i,j}$を用いて関数
$\displaystyle
g(x,y)=\sum_{i=0}^{3}\sum_{j=0}^{3}a_{i,j}x^i y^j
$
によって行います[1]。誤差は$O(h^4)$です。

Fortranプログラムはこちら。

program main
implicit none
integer::i,j
integer::Nx,Ny
double precision::x,y
double precision::xa,xb,ya,yb
double precision,allocatable::x0(:),y0(:),z0(:,:)
double precision,external::f,ip16

Nx=20
Ny=20
allocate(x0(0:Nx), y0(0:Ny), z0(0:Nx,0:Ny))

xa=-3d0
xb= 3d0
ya=-3d0
yb= 3d0

! generate references points
do i=0,Nx
x0(i)=dble(i)*(xb-xa)/Nx+xa
enddo
do j=0,Ny
y0(j)=dble(j)*(yb-ya)/Ny+ya
enddo
do i=0,Nx
do j=0,Ny
z0(i,j)=f(x0(i),y0(j))
write(10,*)x0(i),y0(j),z0(i,j)
enddo
write(10,*)
enddo

! Return interpolated results
do i=0,50
do j=0,50
x=i*0.03d0-1d0
y=j*0.03d0-1d0
write(11,*)x,y,ip16(x,y,Nx,Ny,x0,y0,z0)
enddo
write(11,*)
enddo

stop
end program main

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

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

return
end function f

double precision function ip16(x,y,Nx,Ny,x0,y0,z0)
implicit none
integer,intent(in)::Nx,Ny
double precision,intent(in)::x,y,x0(0:Nx),y0(0:Ny),z0(0:Nx,0:Ny)
!
! Bicubic interpolation
! x0,y0 are sorted by ascending order,
! suppose x0 & y0 are equal interval.
! z0(x0,y0)
!
integer::i,j
integer::i0,i1,i2,i3
integer::j0,j1,j2,j3
double precision::tx,u,u2,u3
double precision::ty,v,v2,v3
double precision::a00,a01,a02,a03,a10,a11,a12,a13
double precision::a20,a21,a22,a23,a30,a31,a32,a33
double precision::p(0:3,0:3)

if(Nx.le.3.or.Ny.le.3)then
write(6,*)" *** Error at ip16"
stop
endif

tx = x-x0(0)
i1 = 0
do i=1,Nx-2
tx = x-x0(i)
if(tx.gt.0d0)then
i1 = i
else
exit
endif
enddo
if(i1.eq.0)i1=1
i0=i1-1
i2=i1+1
i3=i1+2

ty = y-y0(0)
j1 = 0
do j=1,Ny-2
ty = y-y0(j)
if(ty.gt.0d0)then
j1 = j
else
exit
endif
enddo
if(j1.eq.0)j1=1
j0=j1-1
j2=j1+1
j3=j1+2

p(0,0) = z0(i0,j0)
p(0,1) = z0(i0,j1)
p(0,2) = z0(i0,j2)
p(0,3) = z0(i0,j3)
p(1,0) = z0(i1,j0)
p(1,1) = z0(i1,j1)
p(1,2) = z0(i1,j2)
p(1,3) = z0(i1,j3)
p(2,0) = z0(i2,j0)
p(2,1) = z0(i2,j1)
p(2,2) = z0(i2,j2)
p(2,3) = z0(i2,j3)
p(3,0) = z0(i3,j0)
p(3,1) = z0(i3,j1)
p(3,2) = z0(i3,j2)
p(3,3) = z0(i3,j3)

a00 = p(1,1)
a01 = -0.5d0*p(1,0) + 0.5d0*p(1,2)
a02 = p(1,0) - 2.5d0*p(1,1) + 2*p(1,2) - 0.5d0*p(1,3)
a03 = -0.50d0*p(1,0) + 1.50d0*p(1,1) - 1.50d0*p(1,2) + 0.5d0*p(1,3)
a10 = -0.50d0*p(0,1) + 0.50d0*p(2,1)
a11 = 0.25d0*p(0,0) - 0.25d0*p(0,2) - 0.25d0*p(2,0) + 0.25d0*p(2,2)
a12 = -0.50d0*p(0,0) + 1.25d0*p(0,1) - p(0,2) + 0.25d0*p(0,3) &
+ 0.50d0*p(2,0) - 1.25d0*p(2,1) + p(2,2) - 0.25d0*p(2,3)
a13 = 0.25d0*p(0,0) - 0.75d0*p(0,1) + 0.75d0*p(0,2) - 0.25d0*p(0,3) &
- 0.25d0*p(2,0) + 0.75d0*p(2,1) - 0.75d0*p(2,2) + 0.25d0*p(2,3)
a20 = p(0,1) - 2.5d0*p(1,1) + 2.d0*p(2,1) - 0.5d0*p(3,1)
a21 = -0.50d0*p(0,0) + 0.5d0*p(0,2) + 1.25d0*p(1,0) - 1.25d0*p(1,2) &
- p(2,0) + p(2,2) + 0.25d0*p(3,0) - 0.25d0*p(3,2)
a22 = p(0,0) - 2.5d0*p(0,1) + 2.0d0*p(0,2) - 0.50d0*p(0,3) &
- 2.5d0*p(1,0) + 6.25d0*p(1,1) - 5d0*p(1,2) + 1.25d0*p(1,3) &
+ 2.0d0*p(2,0) - 5.0d0*p(2,1) + 4d0*p(2,2) - p(2,3) &
- 0.50d0*p(3,0) + 1.25d0*p(3,1) - p(3,2) + 0.25d0*p(3,3)
a23 = -0.50d0*p(0,0) + 1.50d0*p(0,1) - 1.50d0*p(0,2) + 0.5d0*p(0,3) &
+ 1.25d0*p(1,0) - 3.75d0*p(1,1) + 3.75d0*p(1,2) - 1.25d0*p(1,3) &
- p(2,0) + 3*p(2,1) - 3*p(2,2) + p(2,3) &
+ 0.25d0*p(3,0) - 0.75d0*p(3,1) + 0.75d0*p(3,2) - 0.25d0*p(3,3)
a30 = -0.50d0*p(0,1) + 1.50d0*p(1,1) - 1.50d0*p(2,1) + 0.5d0*p(3,1)
a31 = 0.25d0*p(0,0) - 0.25d0*p(0,2) - 0.75d0*p(1,0) + 0.75d0*p(1,2) &
+ 0.75d0*p(2,0) - 0.75d0*p(2,2) - 0.25d0*p(3,0) + 0.25d0*p(3,2)
a32 = -0.50d0*p(0,0) + 1.25d0*p(0,1) - p(0,2) + 0.25d0*p(0,3) &
+ 1.50d0*p(1,0) - 3.75d0*p(1,1) + 3*p(1,2) - 0.75d0*p(1,3) &
- 1.50d0*p(2,0) + 3.75d0*p(2,1) - 3*p(2,2) + 0.75d0*p(2,3) &
+ 0.50d0*p(3,0) - 1.25d0*p(3,1) + p(3,2) - 0.25d0*p(3,3)
a33 = 0.25d0*p(0,0) - 0.75d0*p(0,1) + 0.75d0*p(0,2) - 0.25d0*p(0,3) &
- 0.75d0*p(1,0) + 2.25d0*p(1,1) - 2.25d0*p(1,2) + 0.75d0*p(1,3) &
+ 0.75d0*p(2,0) - 2.25d0*p(2,1) + 2.25d0*p(2,2) - 0.75d0*p(2,3) &
- 0.25d0*p(3,0) + 0.75d0*p(3,1) - 0.75d0*p(3,2) + 0.25d0*p(3,3)

! Parallel translation u,v [0:1]
u = (x-x0(i1))/(x0(i2)-x0(i1))
v = (y-y0(j1))/(y0(j2)-y0(j1))
u2 = u*u
u3 = u2*u
v2 = v*v
v3 = v2*v
ip16= (a00 + a01 * v + a02 * v2 + a03 * v3) &
+ (a10 + a11 * v + a12 * v2 + a13 * v3) * u &
+ (a20 + a21 * v + a22 * v2 + a23 * v3) * u2 &
+ (a30 + a31 * v + a32 * v2 + a33 * v3) * u3

return
end function ip16

補間 (2次元)

Bicubic補間の低次元バージョンです。
補間したい点を囲むように存在する4点を用いて補間します。
一応Four point Formulaと名づけられています[2]。
補間は,既知の係数$a_{i,j}$を用いて関数
$\displaystyle
g(x,y)=\sum_{i=0}^{1}\sum_{j=0}^{1}a_{i,j}x^i y^j
$
によって行います[2]。誤差は$O(h^2)$です。

function ip4(x,y,Nx,Ny,x0,y0,z0)
implicit none
integer,intent(in)::Nx,Ny
double precision,intent(in)::x,y,x0(1:Nx),y0(1:Ny),z0(1:Nx,1:Ny)
double precision::ip4
!
! 2D Interpolation, 4-point formula,
! xy equal distance grid
!
integer::i,ix0,ix1,iy0,iy1
double precision::tx,ty,dx,dy,p,q

tx=1d100
ix0=1
do i=2,Nx
if(abs(x-x0(i)).le.tx)then
tx=abs(x-x0(i))
ix0=i
endif
enddo
tx=1d100
ix1=1
do i=2,Nx
if(i.ne.ix0)then
if(abs(x-x0(i)).le.tx)then
tx=abs(x-x0(i))
ix1=i
endif
endif
enddo

ty=1d100
iy0=1
do i=2,Ny
if(abs(y-y0(i)).le.ty)then
ty=abs(y-y0(i))
iy0=i
endif
enddo
ty=1d100
iy1=1
do i=2,Ny
if(i.ne.iy0)then
if(abs(y-y0(i)).le.ty)then
ty=abs(y-y0(i))
iy1=i
endif
endif
enddo

dx=x0(ix1)-x0(ix0)
dy=y0(iy1)-y0(iy0)
p=(x-x0(ix0))/dx
q=(y-y0(iy0))/dy

ip4=(1d0-p)*(1d0-q)*z0(ix0,iy0) &
+p*(1d0-q)*z0(ix1,iy0) &
+(1d0-p)*q*z0(ix0,iy1) &
+p*q*z0(ix1,iy1)

return
end function ip4

参考文献

[1]Cubic interpolation
[2]Abramowitz and Stegun.
Handbook of Mathematical Functions.

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

LU分解による連立一次方程式の解法

2019年5月24日 sikino コメントする

数値計算ライブラリLapackを使って線形連立一次方程式
$
\begin{equation}
A\mathbf{x}=\mathbf{b}
\end{equation}
$
をLU分解を用いて数値的に解くプログラムを載せます。

詳しくは説明しませんが、逆行列を用いて解く方法と比べ、
行列が疎である時(行列要素にゼロが多い時)に比較的計算量が減らせます。
また、逆行列を求める際に生じる桁落ちの問題を回避することが出来ます。

解法

連立一次方程式
$
\begin{equation}
A\mathbf{x}=\mathbf{b}
\end{equation}
$
を行列$A$のLU分解を利用して解きます。
ここで、未知なのはベクトル$\mathbf{x}$、既知なのは行列$A$とベクトル$\mathbf{b}$です。

$A$がLU分解されているとします。
すなわち、行列$A$に適当な操作を施して、

$
\begin{align}
A=LU
\end{align}
$
と書けているとします。ここで、行列$L, U$は下三角行列、上三角行列であり、
$
\begin{align}
L=
\left(
\begin{array}{*4{>{\displaystyle}c}}
1&0&\cdots &0 \\[1em]
l_{2,1}&1&\cdots &0 \\[0.5em]
\vdots&\vdots&\ddots&\vdots \\[0.3em]
l_{s,1}&l_{s,2}&\cdots &1
\end{array}
\right) ,~~~
U=
\left(
\begin{array}{*4{>{\displaystyle}c}}
u_{1,1}&u_{1,2}&\cdots &u_{1,s} \\[1em]
0&u_{2,2}&\cdots &u_{2,s} \\[0.5em]
\vdots&\vdots&\ddots&\vdots \\[0.3em]
0 & 0 &\cdots &u_{s,s}
\end{array}
\right)
\end{align}
$

という行列です。
もし、行列$A$がLU分解されていれば、

$
\begin{align}
& A\mathbf{x}=\mathbf{b} \\
& LU\mathbf{x} = \mathbf{b} \\
& L\mathbf{w} = \mathbf{b}~~~…(*)
\end{align}
$

と書けます。ここで、$\mathbf{w} \equiv U\mathbf{x}$と置きました。
まず、式(*)を解いて、$\mathbf{w}$を求めます。これは、下三角行列に対する問題なので、前進代入を用いて簡単に解くことが出来ます。

続いて、
$
U\mathbf{x}=\mathbf{w}
$
を後退代入を利用して解$\mathbf{x}$を求めます。

以上がLU分解を利用して連立一次方程式を解く方法です。
数値計算のルーチンは大きく2つのステップに分かれており、

行列$A$のLU分解を求める(ルーチン dgetrf)
LU分解された行列$A$を用いて連立一次方程式を解く(ルーチン dgetrs)

というステップです。
これは、例えば問題
$
\begin{equation}
A\mathbf{x}=\mathbf{b}
\end{equation}
$
と
$
\begin{equation}
A\mathbf{x}=\mathbf{b’}
\end{equation}
$
の右辺だけが変わる複数の問題を解きたい場合なんかに便利です。この問題の場合、行列$A$は変わらないため、LU分解した結果を両方の問題に流用することが出来るのです。

プログラム

プログラムは以下の通りです。
いつもと同じように、ワーク配列を減らす為だけに導入したルーチンを挟んでいます。

program main
implicit none
integer::i,N
double precision,allocatable::a(:,:),b(:),x(:)
integer,allocatable::ipiv(:)

N=3
allocate(a(1:N,1:N))
allocate(ipiv(1:N),b(1:N),x(1:N))
a(1,1:N)=(/1d0,3d0,5d0/)
a(2,1:N)=(/0d0,3d0,1d0/)
a(3,1:N)=(/6d0,2d0,5d0/)
b(1:N)=(/33d0,10d0,66d0/)

do i=1,N
write(6,'(3f10.5)')a(i,1:N)
enddo

call LUfact(N,a,ipiv)
call axbsolve(N,a,ipiv,b,x)

write(6,*)"---- solution ----"
do i=1,N
write(6,'(2f10.5)')x(i)
enddo

stop
end program main

subroutine LUfact(N,LU,ipiv)
implicit none
integer,intent(in)::N
double precision,intent(inout)::LU(1:N,1:N)
integer,intent(out)::ipiv(1:N)

! LU factorization for lapack
! Overwrite matrix LU by factorized LU matrix

integer::m,lda,info

m=N
lda=N
call dgetrf(m,N,LU,lda,ipiv,info)
if(info.ne.0)then
write(6,*)"**error at LUlapack"
stop
endif

return
end subroutine LUfact

subroutine axbsolve(N,LU,ipiv,b,x)
implicit none
integer,intent(in)::N
integer,intent(in)::ipiv(1:N)
double precision,intent(in)::LU(1:N,1:N),b(1:N)
double precision,intent(out)::x(1:N)
!
! Solve simultaneous equations by Lapack
!
double precision,allocatable::bl(:,:)
integer::nrhs,lda,ldb,info
character(1)::trans

nrhs=1
trans='N'
allocate(bl(1:N,1:nrhs))
bl(1:N,1)=b(1:N)
lda=N
ldb=N
call dgetrs(trans,N,nrhs,LU,lda,ipiv,bl,ldb,info)
if(info.ne.0)then
write(6,*)"**error at LUlapack"
stop
endif

x(1:N)=bl(1:N,1)

return
end subroutine axbsolve

上のプログラムをlapackと一緒にコンパイルして動かすと、端末上で

> ./a.out
1.00000 3.00000 5.00000
0.00000 3.00000 1.00000
6.00000 2.00000 5.00000
-------------
7.00000
2.00000
4.00000
>

という結果が得られるかと思います。
これは、連立一次方程式
$
\begin{align}
\left(
\begin{array}{ccc}
1& 3& 5 \\
0& 3& 1 \\
6& 2& 5
\end{array}
\right)
\left(
\begin{array}{ccc}
x \\
y \\
z
\end{array}
\right)
=
\left(
\begin{array}{ccc}
33 \\
10 \\
66
\end{array}
\right)
\end{align}
$
を解いた結果です。
解析解($x=7, y=2, z=4$)は、例えばwolframを利用して求められます。
{{1,3,5},{0,3,1},{6,2,5}}.{x,y,z}=={33,10,66} –wolfram alpha

プログラミングと数値計算, 物理学

質点と壁との反発を表す運動方程式

2019年4月30日 sikino コメントする

質点が壁に衝突し、反発することを数式で表現します。

壁の定義
運動方程式の導出
プログラム
実行結果
参考文献

壁の定義

壁とは、壁の法線方向に平行な質点の速度成分を反転させるデバイスである。
と定義します。

運動方程式の導出

まず議論を簡単にする為に一次元の運動について考え、その後多次元の運動について定式化をしていきます。

1次元の壁との衝突を表す運動方程式

ポテンシャル$V(x)$保存力下の質点の運動を考えます。
壁が
$f(x,t)=0$
で表現されているとします。
壁との衝突では位置は連続、速度は不連続な振る舞いをすると考えると、壁に衝突した時、力はデルタ関数の振る舞いをしていることが予想できます。
よって、未知の定数$c$を用いて、

と書くことが出来ます。
ここで、時刻$t＝t_i$は、
$
f(x(t),t)=0~~\to~~t=t_i
$
を満たす時刻(壁との衝突時刻)です。$f(x,t)$の時間微分は、

と書くことが出来ます。
$c$を定めるために式(2)の両辺を、時刻$t=t_j$周りを微小時間$\Delta$で積分します。すると

を得ます。ここで、保存力の時間変化は連続であることを仮定します。すなわち、

がいかなる時刻$t$で成立するとします。計算を進めると、

という結果が導かれます。ここで、$\delta_{i,j}$はクロネッカーのデルタを表します。式(6)を式(2)の右辺に代入し、$c$を消去すると、運動方程式

を得ます。式(7)はこのままでは解くことが出来ません。なぜなら、未来の時刻の速度$v(t_i+0)$が含まれているからです。これをどうにかするには、新たな条件式、すなわち反発前後の条件式が必要になります。

壁と衝突する場合、質点の衝突前後の速度$v(t_i-0), v(t_i+0)$と、壁の衝突前後の速度$v_\text{w}(t_i-0), v_\text{w}(t_i+0)$の間には、反発係数$e$を用いて

の関係があると予想します。今、壁の質量が無限大であり、衝突前後で速度変化がない場合を考えると、壁の速度は、質点との衝突によって変化しないと考えます。すなわち

が成立していると考えます。すると、衝突後の質点の速度は、$e$を用いて

と書けます。よって、式(7)に代入して、

を得ます。

最後に、位置$x(t_i)$における壁の速度$v_\text{w}(t_i)$を考えましょう。衝突の前後のごく短い時間では、質点の動きは壁に追従すると考えます。すなわち、

が成り立っているとします。衝突時には質点の位置$x$と衝突位置は同じであることを注記します。式(12)を書き換えると、衝突が起こる時刻$t=t_i$の周りで

が成立しています。このことから、衝突前後のごく短い時間では

が成立します。よって、

から、壁の速度

を得ます。よって、壁との衝突を記述する運動方程式は

となります。

多次元の壁との衝突を表す運動方程式

壁

と質点が衝突することを考えます。運動方程式は

であり、1次元の場合と同様に時刻$t=t_j$周りの微小時間で積分して

を得ます。式(20)に代入すれば運動方程式

を得ます。次の節で述べる結果(多次元の反発)を先取りすると、式(22)は

と変形することが出来ます。ここで、$e$は反発係数、$\mathbf{v}_\text{w}(t)$は壁の速度ベクトル、$\mathbf{n}$は壁の法線ベクトルであり、$\mathbf{n}$は

と表されます。
実際に数値計算を行う際には、衝突時刻$t_i$と位置$\mathbf{r}=\mathbf{r}_i$を求めた後、

に従って速度ベクトルを変更すると良いと思います（根拠は特にありません）。

多次元の反発

衝突時刻$t=t_i$, 点$\mathbf{r}(t_i)$で質点が$f(\mathbf{r},t)=0$で表される壁に衝突することを考えます。
衝突前、後の質点の速度ベクトルをそれぞれ$\mathbf{v}(t_i-0), \mathbf{v}(t_i+0),$と置きます。
位置$\mathbf{r}_i$における壁の単位法線ベクトルを$\mathbf{n}, (\mathbf{n}^2=1)$と書くと、

と書くことが出来ます。ここで、$c$は未知の定数でこれから定めていきます。
衝突では壁の法線方向の速度成分のみが変化すると仮定しているので、

が成立します。ここで、$e$は反発係数で$0\leq e\leq 1$である($e=0$:完全非弾性衝突、$e=1$:完全弾性衝突。
また、$v_\perp, v_{\text{w}\perp}$は$\mathbf{v}, \mathbf{v}_\text{w}$の、壁の法線方向の成分であり、

と書きあらわすことが出来ます。
式(27)の$c$を求めましょう。式(27)の両辺に$\mathbf{n}^\mathsf{T}$を掛けて

のように得ます。よって衝突後の質点の速度は

と表すことが出来ます。壁が$f(\mathbf{r},t)=0$を満たす線と表現されていれば、時刻$t=t_i$で$\mathbf{n}$は

と求める事が出来ます。
続いて壁面の速度$\mathbf{v}_\text{w}(t)$を考えます。1次元の場合と同様に衝突前後のごく短い時間では質点の動きは壁に追従すると予想します。すなわち、

が成立すると考えます。ここで、式(35)の$\mathbf{r}(t)$は質点の位置ベクトルです。すなわち、衝突が起こる時刻$t=t_i$の周りで

が成立しており、衝突のごく短い時間では

が得られます。よって

を得ますので、壁の速度

を得ることが出来ます。

プログラム

2次元平面の運動に対するFortran90のプログラムはこちらです。
時間発展は、刻み幅制御陽的ルンゲクッタ法
衝突時刻を求める際の根の探索には、Anderson-Björk’s法を用いています。

壁の関数(サブルーチン fw)とその偏微分(サブルーチン pwf)は手で入力しています。

壁との衝突は、関数の”符号が変わった時”で判定しているので、法線の方向はどちらでも構いません。
すなわち、
$
f(x,y,t)=\pm g(x,y,t)
$
は同じです。

プログラム

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

program main
implicit none
integer::i,j,N,info,Nt
double precision::x,h,rktol,xbound,xstart
double precision::xa,xb,fwa,fwb
double precision,allocatable::y(:),ya(:),xarr(:)
double precision,external::fw
external::grk

double precision::pfx,pfy,pft,pfabs,nx,ny,els,cv
els=1d0
rktol=1d-8

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

xstart=0d0; xbound=20d0
y(1)=4d0 ! x (0)
y(2)=0d0 ! x'(0)
y(3)=4d0 ! y (0)
y(4)=5d0 ! y'(0)

Nt=201
allocate(xarr(1:Nt))
do j=1,Nt
xarr(j) = (j-1)*(xbound-xstart)/(Nt-1) + xstart
enddo

x = xarr(1)
xa = x
ya(1:N) = y(1:N)
fwa=fw(ya(1),ya(3),xa)

write(11,'(6e25.10e3)')xarr(1),y(1),y(2)&
,y(3),y(4),0.5d0*(y(2)**2+y(4)**2)+9.8*y(3)
i=1
do j=2,Nt
h=1d-6
info=0
xbound=xarr(j)
do while(info.le.0)
call drkf45(grk,x,h,N,y,xbound,info,rktol)
xb = x
fwb = fw(y(1),y(3),x)

if(i.ge.2.and.fwa*fwb.le.0d0)then
! Time when the mass point reflected by the wall
x = xa; y(1:N) = ya(1:N)
call ABstep(fw,N,grk,x,y,fwa,xb,fwb,rktol)
i=0; h=1d-6
if(abs(xbound-x).gt.1d-14)info=0

! Change the velocity
call pfw(y(1),y(3),x,pfx,pfy,pft)
pfabs = sqrt(pfx**2+pfy**2)
nx = pfx/pfabs
ny = pfy/pfabs
cv = -(1d0+els)*(nx*y(2)+ny*y(4)+pft/pfabs)
y(2) = y(2) + cv*nx
y(4) = y(4) + cv*ny
endif

write(10,'(6e25.10e3)')x,y(1),y(2)&
,y(3),y(4),0.5d0*(y(2)**2+y(4)**2)+9.8*y(3)
xa = x
ya(1:N) = y(1:N)
fwa = fwb
i = i+1
enddo
write(11,'(6e25.10e3)')xarr(j),y(1),y(2)&
,y(3),y(4),0.5d0*(y(2)**2+y(4)**2) +9.8*y(3)
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)
double precision::g

g=9.8d0
f(1)=y(2)
f(2)=0d0
f(3)=y(4)
f(4)=-g

return
end subroutine grk

module GBLwall
implicit none
double precision,parameter::k=1d0
double precision,parameter::w=1d0
double precision,parameter::lx=5d0
double precision,parameter::ly=3d0
double precision,parameter::tphi=0d0 !dacos(-1d0)/2d0
double precision,parameter::xphi=-dacos(-1d0)/2d0
end module GBLwall

function fw(x,y,t)
use GBLwall
implicit none
double precision,intent(in)::x,y,t
double precision::fw

! example 1)
fw = y - sin(k*x+xphi)*sin(w*t+tphi)

! example 2)
!fw = y - 2d0*t

! example 3)
!fw = (x/lx)**2 + (y/ly)**2 - 1d0

! example 4)
!fw = (y-x)**2-1d0

return
end function fw

subroutine pfw(x,y,t,pfx,pfy,pft)
use GBLwall
implicit none
double precision,intent(in)::x,y,t
double precision,intent(out)::pfx,pfy,pft

! example 1)
pfx = -k*cos(k*x+xphi)*sin(w*t+tphi)
pfy = 1d0
pft = -w*sin(k*x+xphi)*cos(w*t+tphi)

! example 2)
!pfx=0d0
!pfy=1d0
!pft=-2d0

! example 3)
!pfx=2d0*x/(lx**2)
!pfy=2d0*y/(ly**2)
!pft=0d0

! example 4)
!pfx=-2*(y-x)
!pfy= 2*(y-x)
!pft=0d0

return
end subroutine pfw

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

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.2d0
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/dble(N))/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
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

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

subroutine ABstep(fw,N,grk,x,y,fwa,xb,fwb,rktol)
implicit none
integer,intent(in)::N
double precision,intent(inout)::x,y(1:N)
double precision,intent(inout)::fwa,xb,fwb
double precision,intent(in)::rktol
double precision,external::fw
external::grk

integer::k,info
integer,parameter::ktmax=50
double precision::xa,fwc,m,xc,h,tol1,xm
double precision,allocatable::ya(:)
double precision,parameter::eps=epsilon(1d0)
double precision,parameter::tol=1d-6

! Must satisfy fwa*fwb < 0
allocate(ya(1:N))
ya(1:N) = y(1:N)
xa = x

! Rootfind by (Anderson & Bj"ork's method)
do k=1,ktmax
x = xa
y(1:N) = ya(1:N)

xc = xa-fwa*(xa-xb)/(fwa-fwb)
h = xc-x
info = 0
do while(info.le.0)
call drkf45(grk,x,h,N,y,xc,info,rktol)
enddo
fwc = fw(y(1),y(3),xc)

if(fwa*fwc.gt.0d0)then
! Substitute a by c
if(fwc/fwa.lt.1d0)then
m=1d0-fwc/fwa
else
m=0.5d0
endif
fwb=m*fwb
xa = xc
ya(1:N) = y(1:N)
fwa = fwc
else
! Substitute b by c
if(fwc/fwb.lt.1d0)then
m=1d0-fwc/fwb
else
m=0.5d0
endif
fwa=m*fwa
xb = xc
fwb = fwc
endif

tol1 = 2.0d0*eps*dabs(xc)+0.5d0*tol
xm = 0.5d0*(xa-xb)
if ((dabs(xm).le.tol1).or.(fwc.eq.0d0))exit
enddo

return
end subroutine ABstep

gnuplot用のスクリプト

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

実行結果

プログラムを実行した結果です。
適当な壁を定義して実行しています。

重力なしの運動、動かない壁

楕円の焦点から放たれた質点の軌跡

$
\begin{align}
x(0)&=4,~~x'(0)=(\text{const}) \\
y(0)&=0,~~y'(0)=(\text{const}’)
\end{align}
$

壁
$
\begin{gather}
f(x,y,t)=\left(\frac{x}{l_x}\right)^2+\left(\frac{y}{l_y}\right)^2-1=0,\\
l_x=5,~l_y=3
\end{gather}
$

重力下の運動、動かない壁

$
\begin{align}
x(0)&=0,~~x'(0)=2 \\
y(0)&=0,~~y'(0)=7
\end{align}
$

壁
$
\begin{align}
f(x,y,t)=y – \sin(x)=0
\end{align}
$

重力下の運動、動く壁

$
\begin{align}
x(0)&=4,~~x'(0)=0 \\
y(0)&=4,~~y'(0)=5
\end{align}
$

壁
$
\begin{gather}
f(x,y,t)=y – \sin(kx+\phi_x)\sin(\omega t+\phi_t)=0,\\
k=1,~w=1,~\phi_x=-\pi/2,~\phi_t=0
\end{gather}
$

参考文献

3. 壁との衝突 -物理学の見つけ方

4. 動く壁との衝突 -物理学の見つけ方

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

Hyper-dual numbersによる二階偏微分の計算

2019年4月29日 sikino コメントする

Hyper-dual numbersと呼ばれる、実数を拡張した考えを導入すると導関数が計算できます。

あらかじめ、引数がHyper-dual numbersである時の数々の関数の定義を実装しておけば、その関数の組み合わせで作られる関数の一階導関数、二階導関数のほぼ厳密な答えを得ることが出来ます。
この考えに従う関数の微分方法は、Forward型の自動微分と呼ばれます。

Dual number
1. Dual numberと導関数
Hyper-dual numbers
参考文献

Dual number

Dual numberという考えがあります[2,3,5]。

これは実数を拡張する、という考えで複素数に似た考え方です。

良く知られている実数の拡張方法の一つは、複素数です。
通常の実数に$i^2=-1$を満たす$i$という数を付加するのが複素空間です。

拡張の方法は何も$i$だけではありません。
例えば、$\epsilon^2=0, (\epsilon\ne 0)$を満たす$\epsilon$という数を付加することもできます。

この$\epsilon$は複素数ではありません。
複素平面上でこの性質を満たす数は無いことからも、この$\epsilon$を追加するということは新しい方向への数の拡張です。
実際、複素平面上で$\epsilon^2=0, (\epsilon\ne 0)$を満たす数があるのかを探しても
$
\begin{align}
(a+ib)^2=a^2-b^2+i2ab=0
\end{align}
$
であるので、これを満たすのは$a=b=0$しか存在せず、複素平面上の数ではないことが分かります。

そんな新しい方向$\epsilon$の平面で定義された数をDual number(二重数、または双対数)と呼びます。

二重数$a$は実数部と非実数部から構成されており、
$
a=a_0+a_1\epsilon
$
と書くことが出来ます。ここで、$a_0, a_1$は実数であり、 $\epsilon$は、虚数単位$i$に倣って二重数単位とでも名付けておきましょう。

Dual numberと導関数

二重数は関数の導関数と大きな関係がある事を示しましょう。
テーラー展開を行います。
関数$f$の$x$周りの展開は刻み幅$\Delta$を用いて
$\displaystyle
f(x+\Delta)=f(x)+f'(x) \Delta+\frac{1}{2!}f^{\prime\prime}(x)\Delta^2+\cdots
$
と記述することが出来ます。もし仮に$\Delta$が純非実数だとしましょう。すなわち、$\Delta=h\epsilon$とします。ここで、$h$は実数です。テーラー展開の式に代入すれば、
$
\begin{eqnarray}
f(x+h\epsilon)&=& f(x)+f'(x) h\epsilon+\frac{1}{2!}f^{\prime\prime}(x){h}^2\epsilon^2+\cdots \\
&=& f(x)+f'(x) h\epsilon
\end{eqnarray}
$
となるわけです。注記しますが、上の式は$h$の1次で打ち切っているのではなく、厳密にイコールが成り立っています。
この式が言っているのは、関数$f(a), (a\text{は二重数})$を計算し、その非実数部を$h$で割ると関数の導関数になっているということです。

二重数の非実数部を取り出す関数$\text{Dual}$を定義します。
すると、導関数は
$\displaystyle
f'(x)=\frac{1}{h}\text{Dual}(f(x+h\epsilon))
$
として得られます。
言葉で書けば、二重数空間に拡張した$f(x+h\epsilon)$を計算すると、その非実数部に導関数が現れる、ということです。

二重数のプログラムはJeffrey Fikeさんによる[4]にありますので、そちらをご参考ください。

Hyper-dual numbers

さて、ここまでで二重数の概念を簡単に説明し、新しい実数の拡張を行いました。
Dual Numberのままでは高階導関数は得られません。なぜなら、導関数の二次以降は二重数の性質$\epsilon^2=0$によってゼロになるからです。

そこで、若干工夫します。
Dual Numberを二種類用意することで二階微分を得ることが出来ます[2,3]。
すなわち、
$\begin{gather}
a=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2\\
\epsilon_1^2=\epsilon_2^2=(\epsilon_1\epsilon_2)^2=0,~\epsilon_1\ne 0,~~\epsilon_2\ne 0,~~\epsilon_1\epsilon_2 \ne 0
\end{gather}
$
と実数を拡張します。この様に拡張した数$a$をHyper-dual numbersと呼びます[2,3]。

Hyper-dual numbersの演算

Hyper-dual numbersである
$
\begin{align}
a&=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2\\
b&=b_0+b_1\epsilon_1+b_2\epsilon_2+b_3\epsilon_1\epsilon_2
\end{align}
$
を用意します。和、積、商はそれぞれ
$
\begin{align}
a+b&=(a_0+b_0)+(a_1+b_1)\epsilon_1+(a_2+b_2)\epsilon_2+(a_3+b_3)\epsilon_1\epsilon_2\\
ab&=a_0b_0+(a_0b_1+a_1b_0)\epsilon_1+(a_0b_2+a_2b_0)\epsilon_2+(a_0b_3+a_1b_2+a_2b_1+a_3b_0)\epsilon_1\epsilon_2\\
\frac{1}{a}&=\frac{1}{a_0}-\frac{a_1}{a_0^2}\epsilon_1-\frac{a_2}{a_0^2}\epsilon_2+\left(-\frac{a_3}{a_0^2}+\frac{2a_1a_2}{a_0}\right)\epsilon_1\epsilon_2
\end{align}
$
と定義されます[2]。任意の関数は
$
f(x)=f(x_0)+x_1f'(x_0)\epsilon_1+x_2f'(x_0)\epsilon_2
+\left(x_3f'(x_0)+x_1x_2f^{\prime\prime}(x_0)\right)\epsilon_1\epsilon_2
$
として計算することが出来ます。なので、この結果から$\epsilon_1\epsilon_2$の係数として関数の二階微分が得られます。

Hyper-dual numbersのプログラム

実際にプログラムを組みましょう。
Hyper-dual numbersの型を持つ変数はFortran90では定義できません。
なので、構造体を利用して自分で型と、演算を定義します。

Fortran90のプログラムは以下の通りになるかと思います。
基本的な四則演算、基本的な初等関数のHyper-dual numbersの演算をモジュールとして書いています。

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

module Hyperdualmod
implicit none

type Hyperdual
! x = x0 + x1 e1 + x2 e2 + x3 e1 e2
double precision::x0
double precision::x1
double precision::x2
double precision::x3
end type Hyperdual

! Equal =
interface assignment (=)
module procedure Equal_HH
module procedure Equal_HD
end interface assignment (=)

! Unary operator +, -
interface operator (+)
module procedure Plus_HH
end interface operator (+)
interface operator (-)
module procedure Minus_HH
end interface operator (-)

! Addition operator +
interface operator (+)
module procedure Add_HH
module procedure Add_HD
module procedure Add_DH
end interface operator (+)

! Subtraction operator -
interface operator (-)
module procedure Sub_HH
module procedure Sub_HD
module procedure Sub_DH
end interface operator (-)

! Multiply operator -
interface operator (*)
module procedure Mul_HH
module procedure Mul_HD
module procedure Mul_DH
module procedure Mul_HI
module procedure Mul_IH
end interface operator (*)

! Division operator
interface operator (/)
module procedure Div_HH
module procedure Div_HD
module procedure Div_DH
end interface operator (/)

! Power operator
interface operator (**)
module procedure Pow_HI
module procedure Pow_HH
module procedure Pow_HD
module procedure Pow_DH
end interface operator (**)

! Equal logical
interface operator (.eq.)
module procedure eq_HH
module procedure eq_HD
module procedure eq_DH
module procedure eq_HI
module procedure eq_IH
end interface operator (.eq.)

! Not equal logical
interface operator (.ne.)
module procedure ne_HH
module procedure ne_HD
module procedure ne_DH
module procedure ne_HI
module procedure ne_IH
end interface operator (.ne.)

! Less than logical
interface operator (.lt.)
module procedure lt_HH
module procedure lt_HD
module procedure lt_DH
module procedure lt_HI
module procedure lt_IH
end interface operator (.lt.)

! Less or equal logical
interface operator (.le.)
module procedure le_HH
module procedure le_HD
module procedure le_DH
module procedure le_HI
module procedure le_IH
end interface operator (.le.)

! Greater than logical
interface operator (.gt.)
module procedure gt_HH
module procedure gt_HD
module procedure gt_DH
module procedure gt_HI
module procedure gt_IH
end interface operator (.gt.)

! Greater than or equal
interface operator (.ge.)
module procedure ge_HH
module procedure ge_HD
module procedure ge_DH
module procedure ge_HI
module procedure ge_IH
end interface operator (.ge.)

! abs
interface abs
module procedure abs_H
end interface abs

! int
interface int
module procedure int_H
end interface int

! nint
interface nint
module procedure nint_H
end interface nint

! real
interface real
module procedure real_H
end interface real

! sign
interface sign
module procedure sign_HH
module procedure sign_HD
module procedure sign_DH
end interface sign

! sin
interface sin
module procedure sin_H
end interface sin

! cos
interface cos
module procedure cos_H
end interface cos

! tan
interface tan
module procedure tan_H
end interface tan

! sqrt
interface sqrt
module procedure sqrt_H
end interface sqrt

! log
interface log
module procedure log_H
end interface log

! log10
interface log10
module procedure log10_H
end interface log10

! exp
interface exp
module procedure exp_H
end interface exp

! sinh
interface sinh
module procedure sinh_H
end interface sinh

! cosh
interface cosh
module procedure cosh_H
end interface cosh

! tanh
interface tanh
module procedure tanh_H
end interface tanh

! acos
interface acos
module procedure acos_H
end interface acos

! asin
interface asin
module procedure asin_H
end interface asin

! atan
interface atan
module procedure atan_H
end interface atan

! atan2
interface atan2
module procedure atan2_H
end interface atan2

contains
subroutine Equal_HH(res, inp)
implicit none
type(Hyperdual),intent(out) :: res
type(Hyperdual),intent(in) :: inp
res%x0 = inp%x0
res%x1 = inp%x1
res%x2 = inp%x2
res%x3 = inp%x3
end subroutine Equal_HH
subroutine Equal_HD(res, inp)
implicit none
type(Hyperdual),intent(out) :: res
double precision,intent(in) :: inp
res%x0 = inp
res%x1 = 0d0
res%x2 = 0d0
res%x3 = 0d0
end subroutine Equal_HD

!---------------------------------------

function Plus_HH(t1) result (t2)
type(Hyperdual),intent(in) :: t1
type(Hyperdual) :: t2
t2%x0 = t1%x0
t2%x1 = t1%x1
t2%x2 = t1%x2
t2%x3 = t1%x3
end function Plus_HH
function Minus_HH(t1) result (t2)
type(Hyperdual),intent(in) :: t1
type(Hyperdual) :: t2
t2%x0 = -t1%x0
t2%x1 = -t1%x1
t2%x2 = -t1%x2
t2%x3 = -t1%x3
end function Minus_HH

!---------------------------------------

function Add_HH(t1,t2) result (t3)
type(Hyperdual), intent(in) :: t1,t2
type(Hyperdual) :: t3
t3%x0 = t1%x0 + t2%x0
t3%x1 = t1%x1 + t2%x1
t3%x2 = t1%x2 + t2%x2
t3%x3 = t1%x3 + t2%x3
end function Add_HH
function Add_HD(t1,t2) result (t3)
type(Hyperdual), intent(in) :: t1
double precision,intent(in) :: t2
type(Hyperdual) :: t3
t3%x0 = t1%x0 + t2
t3%x1 = t1%x1
t3%x2 = t1%x2
t3%x3 = t1%x3
end function Add_HD
function Add_DH(t1,t2) result (t3)
double precision,intent(in) :: t1
type(Hyperdual), intent(in) :: t2
type(Hyperdual) :: t3
t3%x0 = t1 + t2%x0
t3%x1 = t2%x1
t3%x2 = t2%x2
t3%x3 = t2%x3
end function Add_DH

!---------------------------------------

function Sub_HH(t1,t2) result (t3)
type(Hyperdual), intent(in) :: t1,t2
type(Hyperdual) :: t3
t3%x0 = t1%x0 - t2%x0
t3%x1 = t1%x1 - t2%x1
t3%x2 = t1%x2 - t2%x2
t3%x3 = t1%x3 - t2%x3
end function Sub_HH
function Sub_HD(t1,t2) result (t3)
type(Hyperdual), intent(in) :: t1
double precision,intent(in) :: t2
type(Hyperdual) :: t3
t3%x0 = t1%x0 - t2
t3%x1 = t1%x1
t3%x2 = t1%x2
t3%x3 = t1%x3
end function Sub_HD
function Sub_DH(t1,t2) result (t3)
double precision,intent(in) :: t1
type(Hyperdual), intent(in) :: t2
type(Hyperdual) :: t3
t3%x0 = t1 - t2%x0
t3%x1 = t2%x1
t3%x2 = t2%x2
t3%x3 = t2%x3
end function Sub_DH

!---------------------------------------

function Mul_HH(t1,t2) result (t3)
type(Hyperdual), intent(in) :: t1,t2
type(Hyperdual) :: t3
t3%x0 = t1%x0 * t2%x0
t3%x1 = t1%x0 * t2%x1 + t1%x1 * t2%x0
t3%x2 = t1%x0 * t2%x2 + t1%x2 * t2%x0
t3%x3 = t1%x0 * t2%x3 + t1%x1 * t2%x2 &
+ t1%x2 * t2%x1 + t1%x3 * t2%x0
end function Mul_HH
function Mul_HD(t1,t2) result (t3)
type(Hyperdual), intent(in) :: t1
double precision,intent(in) :: t2
type(Hyperdual) :: t3
t3%x0 = t1%x0 * t2
t3%x1 = t1%x1 * t2
t3%x2 = t1%x2 * t2
t3%x3 = t1%x3 * t2
end function Mul_HD
function Mul_DH(t1,t2) result (t3)
double precision,intent(in) :: t1
type(Hyperdual), intent(in) :: t2
type(Hyperdual) :: t3
t3%x0 = t1 * t2%x0
t3%x1 = t1 * t2%x1
t3%x2 = t1 * t2%x2
t3%x3 = t1 * t2%x3
end function Mul_DH
function Mul_HI(t1,t2) result (t3)
type(Hyperdual), intent(in) :: t1
integer ,intent(in) :: t2
type(Hyperdual) :: t3
t3%x0 = t1%x0 * t2
t3%x1 = t1%x1 * t2
t3%x2 = t1%x2 * t2
t3%x3 = t1%x3 * t2
end function Mul_HI
function Mul_IH(t1,t2) result (t3)
integer ,intent(in) :: t1
type(Hyperdual), intent(in) :: t2
type(Hyperdual) :: t3
t3%x0 = t1 * t2%x0
t3%x1 = t1 * t2%x1
t3%x2 = t1 * t2%x2
t3%x3 = t1 * t2%x3
end function Mul_IH

!-----------------------------------------

function Div_HH(t1,t2) result (t3)
type(Hyperdual), intent(in) :: t1,t2
type(Hyperdual) :: s2, t3
double precision:: u2
s2%x0 = 1d0/t2%x0
u2 = s2%x0 * s2%x0
s2%x1 = - t2%x1 * u2
s2%x2 = - t2%x2 * u2
s2%x3 = (- t2%x3 + 2d0 * t2%x1 * t2%x2 * s2%x0) * u2
t3%x0 = t1%x0 * s2%x0
t3%x1 = t1%x0 * s2%x1 + t1%x1 * s2%x0
t3%x2 = t1%x0 * s2%x2 + t1%x2 * s2%x0
t3%x3 = t1%x0 * s2%x3 + t1%x1 * s2%x2 &
+ t1%x2 * s2%x1 + t1%x3 * s2%x0
end function Div_HH
function Div_HD(t1,t2) result (t3)
type(Hyperdual), intent(in) :: t1
double precision,intent(in) :: t2
type(Hyperdual) :: t3
double precision :: s2
s2 = 1d0/t2
t3%x0 = t1%x0 * s2
t3%x1 = t1%x1 * s2
t3%x2 = t1%x2 * s2
t3%x3 = t1%x3 * s2
end function Div_HD
function Div_DH(t1,t2) result (t3)
double precision,intent(in) :: t1
type(Hyperdual), intent(in) :: t2
type(Hyperdual) :: s2, t3
double precision :: u2
s2%x0 = 1d0/t2%x0
u2 = s2%x0 * s2%x0
s2%x1 = - t2%x1 * u2
s2%x2 = - t2%x2 * u2
s2%x3 = (- t2%x3 + 2d0 * t2%x1 * t2%x2 * s2%x0) * u2
t3%x0 = t1 * s2%x0
t3%x1 = t1 * s2%x1
t3%x2 = t1 * s2%x2
t3%x3 = t1 * s2%x3
end function Div_DH

!--------------------------

function Pow_HI(t1,t2) result (t3)
type(Hyperdual),intent(in) :: t1
integer,intent(in) :: t2
integer :: i, m
type(Hyperdual) :: t3
t3 = 1d0
m = abs(t2)
do i=1,m
t3 = t3*t1
enddo
if(t2 .lt. 0) t3 = 1d0/t3
end function Pow_HI
function Pow_HH(t1,t2) result (t3)
type(Hyperdual),intent(in) :: t1, t2
type(Hyperdual) :: t3, v4
v4 = log_H(t1)
t3 = exp_H(t2*v4)
end function Pow_HH
function Pow_DH(t1,t2) result (t3)
double precision,intent(in) :: t1
type(Hyperdual), intent(in) :: t2
double precision :: v4
type(Hyperdual) :: t3
v4 = log(t1)
t3 = exp_H(t2*v4)
end function Pow_DH
function Pow_HD(t1,t2) result (t3)
type(Hyperdual), intent(in) :: t1
double precision,intent(in) :: t2
type(Hyperdual) :: t3
double precision,parameter :: tol=1d-30
double precision :: tx, p1, p2
tx = t1%x0
if(abs(tx) .lt. tol) then
if(tx .ge. 0d0) then
tx = tol
else
tx = -tol
endif
endif
p1 = t2*(tx**(t2-1d0))
p2 = t2*(t2-1d0)*(tx**(t2-2d0))
t3%x0 = (t1%x0)**t2
t3%x1 = t1%x1 *p1
t3%x2 = t1%x2 *p1
t3%x3 = t1%x3 *p1 + t1%x1 * t1%x2 *p2
end function Pow_HD

!---------------------------------

! .eq.
logical function eq_HH(lhs, rhs)
type(Hyperdual),intent(in)::lhs, rhs
eq_HH = lhs%x0 == rhs%x0
end function eq_HH
logical function eq_HD(lhs, rhs)
type(Hyperdual),intent(in)::lhs
double precision,intent(in)::rhs
eq_HD = lhs%x0 == rhs
end function eq_HD
logical function eq_DH(lhs, rhs)
double precision,intent(in)::lhs
type(Hyperdual),intent(in)::rhs
eq_DH = lhs == rhs%x0
end function eq_DH
logical function eq_HI(lhs, rhs)
type(Hyperdual),intent(in)::lhs
integer,intent(in)::rhs
eq_HI = lhs%x0 == rhs
end function eq_HI
logical function eq_IH(lhs, rhs)
integer,intent(in)::lhs
type(Hyperdual),intent(in)::rhs
eq_IH = lhs == rhs%x0
end function eq_IH

!------------------------

! .ne.
logical function ne_HH(lhs, rhs)
type(Hyperdual),intent(in)::lhs, rhs
ne_HH = lhs%x0 /= rhs%x0
end function ne_HH
logical function ne_HD(lhs, rhs)
type(Hyperdual),intent(in)::lhs
double precision,intent(in)::rhs
ne_HD = lhs%x0 /= rhs
end function ne_HD
logical function ne_DH(lhs, rhs)
double precision,intent(in)::lhs
type(Hyperdual),intent(in)::rhs
ne_DH = lhs /= rhs%x0
end function ne_DH
logical function ne_HI(lhs, rhs)
type(Hyperdual),intent(in)::lhs
integer,intent(in)::rhs
ne_HI = lhs%x0 /= rhs
end function ne_HI
logical function ne_IH(lhs, rhs)
integer,intent(in)::lhs
type(Hyperdual),intent(in)::rhs
ne_IH = lhs /= rhs%x0
end function ne_IH

!-------------------

! .lt.
logical function lt_HH(lhs, rhs)
type(Hyperdual),intent(in)::lhs, rhs
lt_HH = lhs%x0 < rhs%x0
end function lt_HH
logical function lt_HD(lhs, rhs)
type(Hyperdual),intent(in)::lhs
double precision,intent(in)::rhs
lt_HD = lhs%x0 < rhs
end function lt_HD
logical function lt_DH(lhs, rhs)
double precision,intent(in)::lhs
type(Hyperdual),intent(in)::rhs
lt_DH = lhs < rhs%x0
end function lt_DH
logical function lt_HI(lhs, rhs)
type(Hyperdual),intent(in)::lhs
integer,intent(in)::rhs
lt_HI = lhs%x0 < rhs
end function lt_HI
logical function lt_IH(lhs, rhs)
integer,intent(in)::lhs
type(Hyperdual),intent(in)::rhs
lt_IH = lhs < rhs%x0
end function lt_IH

!-------------------

!.le.
logical function le_HH(lhs, rhs)
type(Hyperdual),intent(in)::lhs, rhs
le_HH = lhs%x0 <= rhs%x0
end function le_HH
logical function le_HD(lhs, rhs)
type(Hyperdual),intent(in)::lhs
double precision,intent(in)::rhs
le_HD = lhs%x0 <= rhs
end function le_HD
logical function le_DH(lhs, rhs)
double precision,intent(in)::lhs
type(Hyperdual),intent(in)::rhs
le_DH = lhs <= rhs%x0
end function le_DH
logical function le_HI(lhs, rhs)
type(Hyperdual),intent(in)::lhs
integer,intent(in)::rhs
le_HI = lhs%x0 <= rhs
end function le_HI
logical function le_IH(lhs, rhs)
integer,intent(in)::lhs
type(Hyperdual),intent(in)::rhs
le_IH = lhs <= rhs%x0
end function le_IH

!-------------------

! .gt.
logical function gt_HH(lhs, rhs)
type(Hyperdual),intent(in)::lhs, rhs
gt_HH = lhs%x0 > rhs%x0
end function gt_HH
logical function gt_HD(lhs, rhs)
type(Hyperdual),intent(in)::lhs
double precision,intent(in)::rhs
gt_HD = lhs%x0 > rhs
end function gt_HD
logical function gt_DH(lhs, rhs)
double precision,intent(in)::lhs
type(Hyperdual),intent(in)::rhs
gt_DH = lhs > rhs%x0
end function gt_DH
logical function gt_HI(lhs, rhs)
type(Hyperdual),intent(in)::lhs
integer,intent(in)::rhs
gt_HI = lhs%x0 > rhs
end function gt_HI
logical function gt_IH(lhs, rhs)
integer,intent(in)::lhs
type(Hyperdual),intent(in)::rhs
gt_IH = lhs > rhs%x0
end function gt_IH

!----------------------------

! .ge.
logical function ge_HH(lhs, rhs)
type(Hyperdual),intent(in)::lhs, rhs
ge_HH = lhs%x0 >= rhs%x0
end function ge_HH
logical function ge_HD(lhs, rhs)
type(Hyperdual),intent(in)::lhs
double precision,intent(in)::rhs
ge_HD = lhs%x0 >= rhs
end function ge_HD
logical function ge_DH(lhs, rhs)
double precision,intent(in)::lhs
type(Hyperdual),intent(in)::rhs
ge_DH = lhs >= rhs%x0
end function ge_DH
logical function ge_HI(lhs, rhs)
type(Hyperdual),intent(in)::lhs
integer,intent(in)::rhs
ge_HI = lhs%x0 >= rhs
end function ge_HI
logical function ge_IH(lhs, rhs)
integer,intent(in)::lhs
type(Hyperdual),intent(in)::rhs
ge_IH = lhs >= rhs%x0
end function ge_IH

!------------------------------

! Absolute
function abs_H(t1) result (t2)
type(Hyperdual),intent(in)::t1
type(Hyperdual)::t2
if(t1%x0 .ge. 0d0) then
t2%x0 = t1%x0
t2%x1 = t1%x1
t2%x2 = t1%x2
t2%x3 = t1%x3
else
t2%x0 = -t1%x0
t2%x1 = -t1%x1
t2%x2 = -t1%x2
t2%x3 = -t1%x3
endif
end function abs_H

! Int
function int_H(t1) result (t2)
type(Hyperdual),intent(in)::t1
integer :: t2
t2 = int(t1%x0)
end function int_H

! Nearest int
function nint_H(t1) result (t2)
type(Hyperdual),intent(in)::t1
integer :: t2
t2 = nint(t1%x0)
end function nint_H

! Real
function real_H(t1) result (t2)
type(Hyperdual),intent(in)::t1
double precision::t2
t2 = t1%x0
end function real_H

! Sign
function sign_HH(t1,t2) result (t3)
type(Hyperdual),intent(in)::t1, t2
type(Hyperdual)::t3
double precision::ssign
if(t2%x0 .lt. 0d0) then
ssign = -1d0
else
ssign = 1d0
endif
t3 = ssign*t1
end function sign_HH
function sign_HD(t1,t2) result (t3)
type(Hyperdual),intent(in)::t1
double precision,intent(in) :: t2
type(Hyperdual)::t3
double precision::ssign
if(t2 .lt. 0d0) then
ssign = -1d0
else
ssign = 1d0
endif
t3 = ssign*t1
end function sign_HD
function sign_DH(t1,t2) result (t3)
double precision,intent(in) :: t1
type(Hyperdual),intent(in)::t2
type(Hyperdual)::t3
double precision :: ssign
if(t2%x0 .lt. 0d0) then
ssign = -1d0
else
ssign = 1d0
endif
t3 = ssign*t1
end function sign_DH

! sin
function sin_H(t1) result (t2)
type(Hyperdual),intent(in)::t1
type(Hyperdual)::t2
double precision::df0, df1, df2
df0 = sin(t1%x0)
df1 = cos(t1%x0)
df2 = -df0
t2%x0 = df0
t2%x1 = t1%x1 * df1
t2%x2 = t1%x2 * df1
t2%x3 = t1%x3 * df1 + t1%x1 * t1%x2 * df2
end function sin_H

! cos
function cos_H(t1) result (t2)
type(Hyperdual),intent(in)::t1
type(Hyperdual)::t2
double precision::df0, df1, df2
df0 = cos(t1%x0)
df1 = -sin(t1%x0)
df2 = -df0
t2%x0 = df0
t2%x1 = t1%x1 * df1
t2%x2 = t1%x2 * df1
t2%x3 = t1%x3 * df1 + t1%x1 * t1%x2 * df2
end function cos_H

! tan
function tan_H(t1) result (t2)
type(Hyperdual),intent(in)::t1
type(Hyperdual)::t2
double precision::df0, df1, df2
df0 = tan(t1%x0)
df1 = 1d0 + df0*df0
df2 = 2d0 * df0 * df1
t2%x0 = df0
t2%x1 = t1%x1 * df1
t2%x2 = t1%x2 * df1
t2%x3 = t1%x3 * df1 + t1%x1 * t1%x2 * df2
end function tan_H

! sqrt
function sqrt_H(t1) result (t2)
type(Hyperdual),intent(in)::t1
type(Hyperdual)::t2
double precision::df0, df1, df2
df0 = sqrt(t1%x0)
df1 = 1d0/(2*df0)
df2 = -2d0*df1*df1*df1
t2%x0 = df0
t2%x1 = t1%x1 * df1
t2%x2 = t1%x2 * df1
t2%x3 = t1%x3 * df1 + t1%x1 * t1%x2 * df2
end function sqrt_H

! log
function log_H(t1) result (t2)
type(Hyperdual),intent(in)::t1
type(Hyperdual)::t2
double precision::df0, df1, df2
df0 = log(t1%x0)
df1 = 1d0/t1%x0
df2 = -df1*df1
t2%x0 = df0
t2%x1 = t1%x1 * df1
t2%x2 = t1%x2 * df1
t2%x3 = t1%x3 * df1 + t1%x1 * t1%x2 * df2
end function log_H

! log10
function log10_H(t1) result (t2)
type(Hyperdual),intent(in)::t1
type(Hyperdual)::t2
double precision::df0, df1, df2
df0 = log10(t1%x0)
df1 = 1d0/(log(10d0)*t1%x0)
df2 = -df1/(t1%x0)
t2%x0 = df0
t2%x1 = t1%x1 * df1
t2%x2 = t1%x2 * df1
t2%x3 = t1%x3 * df1 + t1%x1 * t1%x2 * df2
end function log10_H

! exp
function exp_H(t1) result (t2)
type(Hyperdual),intent(in)::t1
type(Hyperdual)::t2
double precision::df0, df1, df2
df0 = exp(t1%x0)
df1 = df0
df2 = df1
t2%x0 = df0
t2%x1 = t1%x1 * df1
t2%x2 = t1%x2 * df1
t2%x3 = t1%x3 * df1 + t1%x1 * t1%x2 * df2
end function exp_H

! sinh
function sinh_H(t1) result (t2)
type(Hyperdual),intent(in)::t1
type(Hyperdual)::u1, u2, t2
u1 = exp(t1)
u2 = exp(-t1)
t2 = 0.5d0*(u1-u2)
end function sinh_H

! cosh
function cosh_H(t1) result (t2)
type(Hyperdual),intent(in)::t1
type(Hyperdual)::u1, u2, t2
u1 = exp(t1)
u2 = exp(-t1)
t2 = 0.5d0*(u1+u2)
end function cosh_H

! tanh
function tanh_H(t1) result (t2)
type(Hyperdual),intent(in)::t1
type(Hyperdual)::u1, u2, t2
u1 = exp(t1)
u2 = exp(-t1)
t2 = (u1-u2)/(u1+u2)
end function tanh_H

! acos
function acos_H(t1) result (t2)
type(Hyperdual),intent(in)::t1
type(Hyperdual)::t2
double precision::df0, df1, df2
df0 = acos(t1%x0)
df1 = -1d0/sqrt(1d0 - t1%x0*t1%x0)
df2 = df1*df1*df1 * t1%x0
t2%x0 = df0
t2%x1 = t1%x1 * df1
t2%x2 = t1%x2 * df1
t2%x3 = t1%x3 * df1 + t1%x1 * t1%x2 * df2
end function acos_H

! asin
function asin_H(t1) result (t2)
type(Hyperdual),intent(in)::t1
type(Hyperdual)::t2
double precision::df0, df1, df2
df0 = asin(t1%x0)
df1 = 1d0/sqrt(1d0 - t1%x0*t1%x0)
df2 = df1*df1*df1 * t1%x0
t2%x0 = df0
t2%x1 = t1%x1 * df1
t2%x2 = t1%x2 * df1
t2%x3 = t1%x3 * df1 + t1%x1 * t1%x2 * df2
end function asin_H

! atan
function atan_H(t1) result (t2)
type(Hyperdual),intent(in)::t1
type(Hyperdual)::t2
double precision::df0, df1, df2
df0 = atan(t1%x0)
df1 = 1d0/(1.0d0 + t1%x0*t1%x0)
df2 = -2d0*t1%x0*df1*df1
t2%x0 = df0
t2%x1 = t1%x1 * df1
t2%x2 = t1%x2 * df1
t2%x3 = t1%x3 * df1 + t1%x1 * t1%x2 * df2
end function atan_H

! atan2
function atan2_H(t1, t2) result (t3)
type(Hyperdual),intent(in)::t1, t2
type(Hyperdual)::t3
double precision::a, b, c, d, e, f, g, h
double precision::df0, r2, fx, fy
!y0 y1 y2 y3
a = t1%x0; b = t1%x1; c = t1%x2; d = t1%x3
!x0 x1 x2 x3
e = t2%x0; f = t2%x1; g = t2%x2; h = t2%x3
r2 = a*a + e*e
fx = -e/r2
fy = a/r2
df0 = atan2(a,e)
t3%x0 = df0
t3%x1 = f*fx + b*fy
t3%x2 = g*fx + c*fy
t3%x3 = h*fx + d*fy + (f*c+g*b)*(fx-fy)*(fx+fy) - 2d0*(f*g+b*c)*fx*fy
end function atan2_H
end module Hyperdualmod

上のモジュールと下のメインプログラムを一緒にコンパイルすることにより、関数の二階微分が得られます。

例として
2変数関数
$\displaystyle
f(x,y)=\frac{\text{ln}(xy^2)e^x}{\sqrt{\sin^3{x}+\sin^3{y}}}
$
の偏微分
$\begin{align}
\frac{\partial f}{\partial x},~~\frac{\partial f}{\partial y},~~\frac{\partial^2 f}{\partial x\partial y}
\end{align}
$
を得ることを考えます。
$
f(x+1\epsilon_1,y+1\epsilon_2)
$
を計算すると
$
f(x+1\epsilon_1,y+1\epsilon_2)=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2
$
のように実数係数$a_0, a_1, a_2, a_3$が得られます。
すると、
$
\begin{align}
\frac{\partial f}{\partial x}=a_1 \\
\frac{\partial f}{\partial y}=a_2 \\
\frac{\partial^2 f}{\partial x\partial y} =a_3
\end{align}
$
として偏微分が得られます。
ちなみに、二階微分が欲しい場合は
$
f(x+1\epsilon_1+1\epsilon_2, y)
$
を考えると
$
\begin{align}
\frac{\partial f}{\partial x}=a_1 \\
\frac{\partial f}{\partial x}=a_2 \\
\frac{\partial^2 f}{\partial x^2} =a_3
\end{align}
$
として得られます。$a_1,~a_2$はどちらを採用しても構いません。

プログラムでは変数の型Hyperdualを持つ入力変数を$\text{xH,yH}$,出力を$\text{rH}$と置いています。

program main
use Hyperdualmod
implicit none
type(Hyperdual)::xH,yH,rH

xH%x0 = 0.3d0 ! real part
xH%x1 = 1d0 ! unreal part \epsilon_1
xH%x2 = 0d0 ! unreal part \epsilon_2
xH%x3 = 0d0 ! unreal part \epsilon_1\epsilon_2

yH%x0 = 0.4d0 ! real part
yH%x1 = 0d0 ! unreal part \epsilon_1
yH%x2 = 1d0 ! unreal part \epsilon_2
yH%x3 = 0d0 ! unreal part \epsilon_1\epsilon_2

write(6,'(4f23.16)')xH%x0,xH%x1,xH%x2,xH%x3
write(6,'(4f23.16)')yH%x0,yH%x1,yH%x2,yH%x3

rH = log(xH*yH**2)*exp(xH)/sqrt(sin(xH)**3+cos(yH)**3)
!rH = asin(2d0*xH)*acos(yH)/atan(xH*yH)
!rH = xH**yH

write(6,'(4f23.16)')rH%x0,rH%x1,rH%x2,rH%x3

stop
end program main

ヘッセ行列

二階偏微分の計算が出来たので、ヘッセ行列が簡単に計算できます。
ルーチンを作れば、以下の通りになるかと思います。
下のプログラムは3変数関数
$\displaystyle
f(x,y,z)=\exp(xy)\tan(z)
$
の$x=-2,~y=3,~z=1$におけるヘッセ行列を計算します。

program main
use Hyperdualmod
implicit none
integer::N
double precision::x,y,z,f
double precision,allocatable::nabla(:),Hesse(:,:)
double precision,allocatable::w(:)
external::func

N=3

allocate(nabla(1:N),Hesse(1:N,1:N))
nabla=0d0
Hesse=0d0
allocate(w(1:N))
w=0d0

x = -2d0
y = 3d0
z = 1d0
w(1) = x
w(2) = y
w(3) = z
call Hessian(N,w,func,nabla,Hesse)

f=exp(x*y)*tan(z)
write(6,'(1e24.16)')nabla(1)
write(6,'(1e24.16)')y*f
write(6,'(1e24.16)')nabla(2)
write(6,'(1e24.16)')x*f
write(6,'(1e24.16)')nabla(3)
write(6,'(1e24.16)')exp(x*y)/(cos(z)**2)
write(6,'(3e24.16)')Hesse(1,1:3)
write(6,'(3e24.16)')f*y**2, f*(1+x*y), f*y/(cos(z)*sin(z))
write(6,'(3e24.16)')Hesse(2,1:3)
write(6,'(3e24.16)')f*(1+x*y), x**2*f, f*x/(cos(z)*sin(z))
write(6,'(3e24.16)')Hesse(3,1:3)
write(6,'(3e24.16)')f*y/(cos(z)*sin(z)), f*x/(cos(z)*sin(z)), f*2d0/(cos(z)**2)

stop
end program main

subroutine func(N,x,f)
use Hyperdualmod
implicit none
integer::N
type(Hyperdual),intent(in)::x(1:N)
type(Hyperdual),intent(out)::f

f = exp(x(1)*x(2))*tan(x(3))

return
end subroutine func

subroutine Hessian(N,x,func,nabla,Hesse)
use Hyperdualmod
implicit none
integer,intent(in)::N
double precision,intent(in)::x(1:N)
double precision,intent(out)::nabla(1:N)
double precision,intent(out)::Hesse(1:N,1:N)
external::func

integer::i,j,k
type(Hyperdual),allocatable::xH(:)
type(Hyperdual)::rH

allocate(xH(1:N))
do i=1,N
xH(i)%x0 = 0d0
xH(i)%x1 = 0d0
xH(i)%x2 = 0d0
xH(i)%x3 = 0d0
enddo

do i=1,N
xH(i)%x0 = x(i)
xH(i)%x3 = 0d0
enddo

do i=1,N
do k=1,N
xH(k)%x1 = 0d0
xH(k)%x2 = 0d0
enddo
xH(i)%x1 = 1d0
xH(i)%x2 = 1d0

call func(N,xH,rH)
nabla(i) = rH%x1
Hesse(i,i) = rH%x3
enddo

do i=1,N
do j=i+1,N
do k=1,N
xH(k)%x1 = 0d0
xH(k)%x2 = 0d0
enddo
xH(i)%x1 = 1d0
xH(i)%x2 = 0d0
xH(j)%x1 = 0d0
xH(j)%x2 = 1d0
call func(N,xH,rH)
Hesse(i,j) = rH%x3
Hesse(j,i) = Hesse(i,j)
enddo
enddo

return
end subroutine Hessian

Hyperdual.f90にモジュールを入れ、メインプログラムをmain.f90に入れたとすると、以下の結果を得ます。

> gfortran Hyperdual.f90 main.f90
> ./a.out
0.1158128336233602E-01
0.1158128336233602E-01
-0.7720855574890680E-02
-0.7720855574890680E-02
0.8491012233306163E-02
0.8491012233306162E-02
0.3474385008700806E-01 -0.1930213893722670E-01 0.2547303669991849E-01
0.3474385008700806E-01 -0.1930213893722670E-01 0.2547303669991849E-01
-0.1930213893722670E-01 0.1544171114978136E-01 -0.1698202446661233E-01
-0.1930213893722670E-01 0.1544171114978136E-01 -0.1698202446661233E-01
0.2547303669991849E-01 -0.1698202446661233E-01 0.2644793608458059E-01
0.2547303669991849E-01 -0.1698202446661233E-01 0.2644793608458058E-01

実行結果の奇数行目はHyper-dual Numberによる計算結果、偶数行目は解析解を表します。
また、6行目までは一階微分、7行目以降はヘッセ行列を表します。

参考文献

[1]関数
$\displaystyle
f(x,y)=\frac{\text{ln}(xy^2)e^x}{\sqrt{\sin^3{x}+\sin^3{y}}}
$
の偏微分$\frac{\partial^2 f}{\partial x\partial y}$
の計算
https://www.wolframalpha.com/input/?i=D%5Bln(x*y%5E2)*e%5E(x)%2Fsqrt(sin%5E3(x)%2Bcos%5E3(y)),+y,x%5D

[2]Jeffrey A. Fike and Juan J. Alonso,~”The Development of Hyper-Dual Numbers for Exact Second-Derivative Calculations”, 49th AIAA Aerospace Sciences Meeting(2011)
http://adl.stanford.edu/hyperdual/fike_aiaa-2011-886_slides.pdf,
J. A. Fike and J. J. Alonso. The Development of Hyper-Dual Numbers for Exact Second Derivative Calculations. AIAA paper 2011-886, 49th AIAA Aerospace Sciences Meeting, January 4-7, 2011.http://adl.stanford.edu/hyperdual/Fike_AIAA-2011-886.pdf

[3]JeffreyA.Fike,~”Derivative Calculations Using Hyper-Dual Numbers”, Sandia National Laboratories (2016)
https://www.osti.gov/servlets/purl/1368722

[4]Jeffrey Fike, Aerospace Design Lab, http://adl.stanford.edu/hyperdual/

[5]松本佳彦, 新しい数をつくる, (2018) http://ymatz.net/assets/docs/20180629-jtpa-slide-mod

Lagrange, プログラミングと数値計算, 物理学

束縛条件下の運動 – 自由度がうまく落とせない運動

2019年3月24日 sikino コメントする

このページは
束縛条件下の運動 – ホロノミックな束縛と非保存力
の続きです。

拘束がある場合の2次元の運動

ラグランジュの方程式を解く際に適切な座標変換が見付からず自由度が落とせない場合を考えます。その時の運動方程式は

の形で止まります。式(75)において、未知の関数は$x(t), y(t), \lambda(t)$の3つであり、条件式は3つなので解くことは出来るはずです。
方針は、式(75)から$\lambda(t)$を消去すれば良いのです。

束縛条件$f(x,y,t)=0$が成立する点$x=x_0, y=y_0$周りで、微小時間$\Delta t$の変化を考えます。
点$(x_0, y_0)$まわりでテーラー展開すれば、

を得ます。よって

が満たされる変化でなければなりません。すなわち、$x,y$の時間変化は式(77)をいつも満たしていなければならないのです。両辺を$\Delta t$で割って極限をとれば、

が成立することが分かります。式(75)を書きかえれば、

を得ます。式(79c)もまた時間変化しても成立し続けなければなりません。よって、左辺の時間微分を取れば、

を得ます。整理して

を得ます。ここで関係式

を用いました。式(79a), (79b)から未定乗数$\lambda(t)$を消去すれば

を得ます。よって、式(81)と式(84)から連立方程式

を立てることが出来ます。これをあらわに解けば、運動方程式

を得ます。ここで、式(86), (87)の右辺の$-\frac{\partial V}{\partial x}, -\frac{\partial V}{\partial y}$は$x, y$方向の保存力であり、右辺の残りの項は束縛条件によって決まる束縛力であることが分かります。

拘束がある場合の3次元の運動

スカラー関数$f$で表される拘束がある条件のもとで, 運動方程式

またはベクトル表記で

について考えます。ここで表記を略すために

の略記を用いました($y, z$方向についても同様)。
また、$\nabla$はナブラ演算子で、

です。2次元の場合と同様に、式(88d)の時間微分から

が導けます。ここで、$\hat{H}f$は関数$f$のヘッセ行列であり、

を表します。式(92)は、式(81)を3次元に拡張した形になっています。

式(88a),(88b),(88c)から未定乗数を消去すると従属な3つの関係式

を得ることが出来ます。この3式で独立なのは2つです。
独立な式として式(94a)と式(94b)を選ぶことにし、式(92)と共に書けば、連立方程式

を得ます。この式をあらわに解けば、運動方程式

を得ることが出来ます。ここで、$n_x, n_y, n_z$は壁の法線ベクトル$\mathbf{n}=\frac{\nabla f}{|\nabla f|}$の$x,y,z$成分です。偏微分ではないことに注意してください。

ベクトル表記であれば

と書くことが出来ます。

単振り子の例

ここでは単振り子を例に挙げ、式(86)に代入した時本当に振り子を表す方程式になっているのか、また束縛条件の選び方に依らず、同じ運動を記述しているのかを例に挙げます。
保存力は重力で

を考えます。

考える束縛条件は
$
f(x,y,t)=x^2+y^2-l^2=0
$
と
$
f(x,y,t)=\sqrt{x^2+y^2}-l=0
$
です。

束縛条件1

束縛条件として

を考えます。この束縛条件に対する偏微分を計算して、

を得ます。これを式(86), (87)に代入すれば、

であり、整理すれば

を得ます。

束縛条件2

考察

同じ束縛条件なので、式(105)と(109)の示す運動は同じになるはずです。
それを示すには、2式の差異である

を示せればよい、ということです。そのために拘束力に対して座標変換

を考えます。右辺と左辺をそれぞれ計算すれば、

となり、同じ運動を記述していることが分かります。

Fortran90によるプログラム

Fortran90で2次元平面を動く質点の束縛運動のプログラムを作ります。
そのプログラムは
https://slpr.sakura.ne.jp/qp/supplement_data/constraint.tar.gz
に置いてあります。

一応補足しておきますが、初期状態の位置は$f(x, y, t)=0$を満たすx,yでなければなりませんし、初期速度も$f(x, y, t)=0$の傾きと一致していなければなりません。
上の条件を満たさずとも何も警告は出さず、計算自体は行われますが、その計算は束縛運動ではなくなります。

上のプログラムの中で、変更すると思われるのは以下の2つです。
1. fp2d
2. fw2d

1. fp2dは、ポテンシャル$V(x,y)$を記述しており、デフォルトのプログラムでは、重力による位置エネルギー
$
V(x,y)=mgy
$
に設定しています。

2. fw2dは、束縛条件$f(x,y, t)=0$を記述しており、デフォルトのプログラムでは、
$
f(x,y,t)=y – [2(cos(1.5*x)-1) + 0.1x^2]=0
$
に設定しています。

もしも非保存力を入れたい場合は、grkの最後の方を変更してください。

数値計算では、
時間発展は刻み幅制御の陽的ルンゲクッタ、
gradient, Hesse行列、束縛条件に関する時間の二階微分は、Hyper-dual numbersを利用してほぼ厳密な値を得ています。

program main
! Author : sikino
! Date : 2019/10/12 (yyyy/mm/dd)
! URL : http://slpr.sakura.ne.jp/qp/constraint-condition3/
use Parameters
use Hyperdualmod
implicit none
integer::i,N,info,Nt
double precision::t,h,tol,ta,tb
double precision,allocatable::x(:),tgrid(:)
external::grk
double precision::E,Ek,Ep

N=4
allocate(x(1:N))

! Time range [ta, tb]
ta = 0d0
tb = 30d0
! Divide time range equally among Nt
Nt = 201
! Tolerance of the RK method
tol=1d-8

! Initial condition
x(1)=0d0 ! x (ta)
x(2)=4d0 ! x'(ta)
x(3)=0d0 ! y (ta)
x(4)=0d0 ! y'(ta)

! Time grid
allocate(tgrid(1:Nt))
do i=1,Nt
tgrid(i) = (i-1)*(tb-ta)/dble(Nt-1)+ta
enddo

do i=2,Nt
info = 0
t = tgrid(i-1)
h = (tgrid(i) - t)*0.1d0
do while(info.le.0)
call drkf45(grk,t,h,N,x,tgrid(i),info,tol)
enddo

! Kinetic energy
Ek = 0.5d0*mass*(x(2)**2+x(4)**2)
! Potential energy
Ep = mass*g*x(3)
! Total energy
E = Ek + Ep
write(10,'(8e25.10e3)')t,x(1),x(2),x(3),x(4),E,Ek,Ep
enddo

stop
end program main

subroutine fp2d(N,xH,fH)
use Parameters, only:mass, g
use Hyperdualmod
implicit none
integer,intent(in)::N
type(Hyperdual),intent(in)::xH(1:N)
type(Hyperdual),intent(out)::fH

! Potential

type(Hyperdual)::x,y,t
x = xH(1)
y = xH(2)
t = xH(3)
! If conservative force,
! fH doesn't depend on the time.
!-----------------------------

fH = mass * g * y

return
end subroutine fp2d

subroutine fw2d(N,xH,fH)
use Hyperdualmod
implicit none
integer,intent(in)::N
type(Hyperdual),intent(in)::xH(1:N)
type(Hyperdual),intent(out)::fH

! Wall of the function

type(Hyperdual)::x,y,t
x = xH(1)
y = xH(2)
t = xH(3)
!-----------------

fH = y - (2d0*(cos(1.5d0*x)-1d0) + 0.1d0*x*x)
! fH = (y-5d0)**2 + x**2 - 25d0

return
end subroutine fw2d

デフォルトのまま計算すると、fort.10に以下の出力が現れます。アニメーションは、同じファイルに入っているanime.pltで実行できます。

エネルギーは計算の範囲ではほとんどありません。全力学的エネルギーは、計算時間内でおおよそ5桁目が変化していました。

参考文献

6. 振り子 -物理学の見つけ方
 9. 自由な座標　-物理学の見つけ方

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

最速のクイックソート（Fortran）

2019年2月21日 sikino コメントする

クイックソートが速いのは分かっていますが、コーディングの仕方によって速度は変わります。
本稿では、ネット上で公開されているどのクイックソートが早いのか調べていきます。

最も早いクイックソートはNUMPACのプログラムでした。

比較プログラム

比較するプログラムはネット上で公開されているクイックソート+αです。
対象は、倍精度実数をソートするプログラム、です。

1. 再帰を用いたクイックソート(f90)

t-nissieの日記：【電脳】Fortranで書いたクイックソート

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

2. 再帰を用いないクイックソート(f90)

[Fortran]再帰を使わないquicksortその２ -fortran66の日記

※ただし、私が割と変えていますので、オリジナルそのものではないということを注記しておきます。この変更によって、実行速度はほぼほぼ変わらないことは確認しています。

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

3. Netlibのクイックソートdsort.f

https://www.netlib.org/slatec/src/dsort.f

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

*DECK DSORT
SUBROUTINE DSORT (DX, DY, N, KFLAG)
C***BEGIN PROLOGUE DSORT
C***PURPOSE Sort an array and optionally make the same interchanges in
C an auxiliary array. The array may be sorted in increasing
C or decreasing order. A slightly modified QUICKSORT
C algorithm is used.
C***LIBRARY SLATEC
C***CATEGORY N6A2B
C***TYPE DOUBLE PRECISION (SSORT-S, DSORT-D, ISORT-I)
C***KEYWORDS SINGLETON QUICKSORT, SORT, SORTING
C***AUTHOR Jones, R. E., (SNLA)
C Wisniewski, J. A., (SNLA)
C***DESCRIPTION
C
C DSORT sorts array DX and optionally makes the same interchanges in
C array DY. The array DX may be sorted in increasing order or
C decreasing order. A slightly modified quicksort algorithm is used.
C
C Description of Parameters
C DX - array of values to be sorted (usually abscissas)
C DY - array to be (optionally) carried along
C N - number of values in array DX to be sorted
C KFLAG - control parameter
C = 2 means sort DX in increasing order and carry DY along.
C = 1 means sort DX in increasing order (ignoring DY)
C = -1 means sort DX in decreasing order (ignoring DY)
C = -2 means sort DX in decreasing order and carry DY along.
C
C***REFERENCES R. C. Singleton, Algorithm 347, An efficient algorithm
C for sorting with minimal storage, Communications of
C the ACM, 12, 3 (1969), pp. 185-187.
C***ROUTINES CALLED XERMSG
C***REVISION HISTORY (YYMMDD)
C 761101 DATE WRITTEN
C 761118 Modified to use the Singleton quicksort algorithm. (JAW)
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890831 Modified array declarations. (WRB)
C 891009 Removed unreferenced statement labels. (WRB)
C 891024 Changed category. (WRB)
C 891024 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
C 901012 Declared all variables; changed X,Y to DX,DY; changed
C code to parallel SSORT. (M. McClain)
C 920501 Reformatted the REFERENCES section. (DWL, WRB)
C 920519 Clarified error messages. (DWL)
C 920801 Declarations section rebuilt and code restructured to use
C IF-THEN-ELSE-ENDIF. (RWC, WRB)
C***END PROLOGUE DSORT
C .. Scalar Arguments ..
INTEGER KFLAG, N
C .. Array Arguments ..
DOUBLE PRECISION DX(*), DY(*)
C .. Local Scalars ..
DOUBLE PRECISION R, T, TT, TTY, TY
INTEGER I, IJ, J, K, KK, L, M, NN
C .. Local Arrays ..
INTEGER IL(21), IU(21)
C .. External Subroutines ..
! EXTERNAL XERMSG
C .. Intrinsic Functions ..
INTRINSIC ABS, INT
C***FIRST EXECUTABLE STATEMENT DSORT
NN = N
IF (NN .LT. 1) THEN
! CALL XERMSG ('SLATEC', 'DSORT',
! + 'The number of values to be sorted is not positive.', 1, 1)
RETURN
ENDIF
C
KK = ABS(KFLAG)
IF (KK.NE.1 .AND. KK.NE.2) THEN
! CALL XERMSG ('SLATEC', 'DSORT',
! + 'The sort control parameter, K, is not 2, 1, -1, or -2.', 2,
! + 1)
RETURN
ENDIF
C
C Alter array DX to get decreasing order if needed
C
IF (KFLAG .LE. -1) THEN
DO 10 I=1,NN
DX(I) = -DX(I)
10 CONTINUE
ENDIF
C
IF (KK .EQ. 2) GO TO 100
C
C Sort DX only
C
M = 1
I = 1
J = NN
R = 0.375D0
C
20 IF (I .EQ. J) GO TO 60
IF (R .LE. 0.5898437D0) THEN
R = R+3.90625D-2
ELSE
R = R-0.21875D0
ENDIF
C
30 K = I
C
C Select a central element of the array and save it in location T
C
IJ = I + INT((J-I)*R)
T = DX(IJ)
C
C If first element of array is greater than T, interchange with T
C
IF (DX(I) .GT. T) THEN
DX(IJ) = DX(I)
DX(I) = T
T = DX(IJ)
ENDIF
L = J
C
C If last element of array is less than than T, interchange with T
C
IF (DX(J) .LT. T) THEN
DX(IJ) = DX(J)
DX(J) = T
T = DX(IJ)
C
C If first element of array is greater than T, interchange with T
C
IF (DX(I) .GT. T) THEN
DX(IJ) = DX(I)
DX(I) = T
T = DX(IJ)
ENDIF
ENDIF
C
C Find an element in the second half of the array which is smaller
C than T
C
40 L = L-1
IF (DX(L) .GT. T) GO TO 40
C
C Find an element in the first half of the array which is greater
C than T
C
50 K = K+1
IF (DX(K) .LT. T) GO TO 50
C
C Interchange these elements
C
IF (K .LE. L) THEN
TT = DX(L)
DX(L) = DX(K)
DX(K) = TT
GO TO 40
ENDIF
C
C Save upper and lower subscripts of the array yet to be sorted
C
IF (L-I .GT. J-K) THEN
IL(M) = I
IU(M) = L
I = K
M = M+1
ELSE
IL(M) = K
IU(M) = J
J = L
M = M+1
ENDIF
GO TO 70
C
C Begin again on another portion of the unsorted array
C
60 M = M-1
IF (M .EQ. 0) GO TO 190
I = IL(M)
J = IU(M)
C
70 IF (J-I .GE. 1) GO TO 30
IF (I .EQ. 1) GO TO 20
I = I-1
C
80 I = I+1
IF (I .EQ. J) GO TO 60
T = DX(I+1)
IF (DX(I) .LE. T) GO TO 80
K = I
C
90 DX(K+1) = DX(K)
K = K-1
IF (T .LT. DX(K)) GO TO 90
DX(K+1) = T
GO TO 80
C
C Sort DX and carry DY along
C
100 M = 1
I = 1
J = NN
R = 0.375D0
C
110 IF (I .EQ. J) GO TO 150
IF (R .LE. 0.5898437D0) THEN
R = R+3.90625D-2
ELSE
R = R-0.21875D0
ENDIF
C
120 K = I
C
C Select a central element of the array and save it in location T
C
IJ = I + INT((J-I)*R)
T = DX(IJ)
TY = DY(IJ)
C
C If first element of array is greater than T, interchange with T
C
IF (DX(I) .GT. T) THEN
DX(IJ) = DX(I)
DX(I) = T
T = DX(IJ)
DY(IJ) = DY(I)
DY(I) = TY
TY = DY(IJ)
ENDIF
L = J
C
C If last element of array is less than T, interchange with T
C
IF (DX(J) .LT. T) THEN
DX(IJ) = DX(J)
DX(J) = T
T = DX(IJ)
DY(IJ) = DY(J)
DY(J) = TY
TY = DY(IJ)
C
C If first element of array is greater than T, interchange with T
C
IF (DX(I) .GT. T) THEN
DX(IJ) = DX(I)
DX(I) = T
T = DX(IJ)
DY(IJ) = DY(I)
DY(I) = TY
TY = DY(IJ)
ENDIF
ENDIF
C
C Find an element in the second half of the array which is smaller
C than T
C
130 L = L-1
IF (DX(L) .GT. T) GO TO 130
C
C Find an element in the first half of the array which is greater
C than T
C
140 K = K+1
IF (DX(K) .LT. T) GO TO 140
C
C Interchange these elements
C
IF (K .LE. L) THEN
TT = DX(L)
DX(L) = DX(K)
DX(K) = TT
TTY = DY(L)
DY(L) = DY(K)
DY(K) = TTY
GO TO 130
ENDIF
C
C Save upper and lower subscripts of the array yet to be sorted
C
IF (L-I .GT. J-K) THEN
IL(M) = I
IU(M) = L
I = K
M = M+1
ELSE
IL(M) = K
IU(M) = J
J = L
M = M+1
ENDIF
GO TO 160
C
C Begin again on another portion of the unsorted array
C
150 M = M-1
IF (M .EQ. 0) GO TO 190
I = IL(M)
J = IU(M)
C
160 IF (J-I .GE. 1) GO TO 120
IF (I .EQ. 1) GO TO 110
I = I-1
C
170 I = I+1
IF (I .EQ. J) GO TO 150
T = DX(I+1)
TY = DY(I+1)
IF (DX(I) .LE. T) GO TO 170
K = I
C
180 DX(K+1) = DX(K)
DY(K+1) = DY(K)
K = K-1
IF (T .LT. DX(K)) GO TO 180
DX(K+1) = T
DY(K+1) = TY
GO TO 170
C
C Clean up
C
190 IF (KFLAG .LE. -1) THEN
DO 200 I=1,NN
DX(I) = -DX(I)
200 CONTINUE
ENDIF
RETURN
END

4. NUMPACのクイックソート(sortdk.f)

ＳＯＲＴＰＡＣＫ（ＳＯＲＴｘＫ，ＳＯＲＴｘｙ，ＳＲＴＶｘｚ）　（スカラー又はベクトルデータの内部ソーティング）

5. ヒープソート

fortran90によるヒープソートとバブルソート -シキノート

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

以上の5つのプログラムを比較していきます。
メインプログラムはこちら↓

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

program main
implicit none
integer::i,n
double precision,allocatable::x(:),x0(:),data(:)
double precision::d
real::t0,t1

integer::j,Nt
integer::k,Nk

Nt=100
allocate(data(0:4))

do k=4,24
n=nint(10**(k*0.25d0))
allocate(x(1:n),x0(1:n))

data(0:4)=0d0
do j=1,Nt
call pre_random
do i=1,n
call random_number(d)
x0(i)=d
enddo

x=x0
call cpu_time(t0)
call heapsort(n,x)
call cpu_time(t1)
data(0)=data(0)+(t1-t0)

x=x0
call cpu_time(t0)
call quicksort1(x,1,n)
call cpu_time(t1)
data(1)=data(1)+(t1-t0)

x=x0
call cpu_time(t0)
call quicksort2(n,x)
call cpu_time(t1)
data(2)=data(2)+(t1-t0)

x=x0
call cpu_time(t0)
call dsort(x,x,n,1)
call cpu_time(t1)
data(3)=data(3)+(t1-t0)

x=x0
call cpu_time(t0)
call sortdk(n,x,0)
call cpu_time(t1)
data(4)=data(4)+(t1-t0)

enddo
write(6,*)k
write(10,*)n,data(0:4)/dble(Nt)
deallocate(x,x0)

enddo

stop
end program

subroutine pre_random
! sikinote
! Date : 2015/03/15
! : 2015/09/07
!
!How to use?
!===================
!call random_number(a)
!===================
! random_number produce
! value between 0~1
!
implicit none
integer::seedsize,c
integer,allocatable::seed(:)

!In fortran90, seed is array.
! To get seedsize, use below.
call random_seed(size=seedsize)

!Allocate seed array.
allocate(seed(1:seedsize))

!Get system time.
call system_clock(count=c)

!Substitute "seed" using system time.
seed=c
!seed=2

!Set "seed" to produce random number obey to system time.
call random_seed(put=seed)

return
end subroutine pre_random

コンパイルは

gfortran -O3 (ソートのプログラム達) main.f90

で行いました。

結果

結果を示します。
横軸にデータ数、縦軸にソートに掛かった時間を示しました。
”ソートに掛かった時間”とは、同じデータ数でソートを100回繰り返した時の1回当たりの時間です。

それぞれ、
赤線：再帰有りクイックソート
青線：再帰無しクイックソート
緑線：Netlibのクイックソート
紫線：NUMPACのクイックソート
黒線：ヒープソート
です。
最も早いのがNUMPAC, 最も遅いのがヒープソートだと分かりました。

ヒープソートもクイックソートも大体$O(n\log n)$ですので、グラフの傾きは両者でほとんど変わりません。再帰有りの方が遅い、という結果は面白いです。
先入観で”再帰は遅い”と考えていたのですが、そんなことはありませんでした。

続いて、同じデータを対数で見てみると以下の通りです。

走査したデータサイズの範囲では変化は有りません。更に大規模になれば、違いが見えてくるかもしれません。

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

ケプラー問題に対する陽的解法と陰的解法

2019年1月27日 sikino コメントする

2体のケプラー問題を数値的に解きます。
ここでは、適切な変換をして求めるのではなく、刻み幅制御されたプログラムで無理やり計算します。

そして、
陽的解法であるルンゲ=クッタ=フェールベルグ法
と
陰的解法であるルンゲ=クッタ=ガウス=ルジャンドル法
の離心率に対する計算量の違いを調べてみます。

Kepler問題

二次元で二体の問題を考えます。運動方程式は
$
\begin{align}
\frac{d^2 x}{dt^2}&=-\frac{x}{(x^2+y^2)^{3/2}} \\
\frac{d^2 y}{dt^2}&=-\frac{y}{(x^2+y^2)^{3/2}}
\end{align}
$
$t=[0,20],~~x(0)=1-e,~~x'(0)=0,~~y(0)=0,~~y'(0)=\sqrt{\frac{1+e}{1-e}}$
です。ここで、$e$は離心率を表します。

この問題に対する解は良く知られていて、
$
\displaystyle (x+e)^2+\frac{y^2}{1-e^2}=1
$
で表され、
$
0\le e\lt 1
$
の時、楕円となります。

関数の評価回数の離心率の依存性

楕円の軌道を持つ範囲において、
計算は離心率$e$が1に近づくほど難しくなります。
なぜなら、原点付近を通過する際に、質点の導関数の変化が大きくなるからです。

使用したプログラムの説明は
陽的解法はhttps://slpr.sakura.ne.jp/qp/runge-kutta-ex/
陰的解法はhttps://slpr.sakura.ne.jp/qp/irkgl-program/
です。実際のプログラムも置いておきます。

離心率を変化させたときの軌道はこんな感じです。

さて、計算速度の評価ですが、関数が何回評価されたかで比較します。
念の為、陰的解法はLU分解もあるので単純な比較は難しいことを注記しておきます。

図の一番上は関数の評価回数の離心率依存性を表しています。縦軸は評価回数、横軸は$1-e$を表します。真ん中の図は一番上と同じですが、軸を対数にとっています。一番下は、$t=20$まで計算した時に$t=0$の初期エネルギーとの相対誤差$|(E(t=20)-E(t=0))/E(t=0)|$を表しています。

特徴的な振る舞いは、離心率に対して陰的解法の評価回数は線形の振る舞い、陽的解法は指数で振る舞っている点です。
これは、質点が原点の近くを通るような問題の際に違いが際立つ事を示しています。
また、エネルギーの保存に関しても陰的解法の方が良いことが分かるでしょう。

一応注意しておきますが、ここでいう陰的解法はルンゲ=クッタ=ガウス=ルジャンドル法の振る舞いです。一般的な陰的解法については話していないことに注意してください。

プログラム

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

module GBL
implicit none
integer::count
end module GBL

program main
use GBL
implicit none
integer::N,Ns,istep,ih0,Jup,info,ih,i,Nx,iJac,j,cirk,cerk
double precision::h,h0,err0,eta0,tol,xa,xb,tx,e
double precision::y1,y2,y3,y4,eneirk,eneerk
integer,allocatable::ipiv(:),epiv(:)
double precision,allocatable::x(:),y(:),Rtol(:),Atol(:)
double precision,allocatable::Jmat(:,:),errJ(:,:),z0(:),Jac(:,:)
double precision,external::egrk
external::grk

N=4 ! Number of 1st order ODEs
allocate(y(1:N),Rtol(1:N),Atol(1:N))
y=0d0; Rtol=0d0; Atol=0d0

Nx=2
allocate(x(1:Nx))
x=0d0

!------------ Initial set up ------------
Ns=N*3 ! Ns=N*s, s means s-stage IRK.
allocate(ipiv(1:Ns),epiv(1:N))
allocate(Jmat(1:Ns,1:Ns),errJ(1:N,1:N),z0(1:Ns),Jac(1:N,1:N))
ipiv=0; epiv=0; Jmat=0d0; errJ=0d0; Jup=0; ih=0
ih0=0; z0=0d0; eta0=0d0; h0=0d0; err0=0d0; iJac=0; Jac=0d0
!----------------------------------------

xa=0d0 ! Initial value of x
xb=20d0 ! End value of x
do i=1,Nx ! Separate equal interval x.
x(i)=(i-1)*(xb-xa)/dble(Nx-1)+xa
enddo

do j=1,7
if(j.eq.1)e=0.9d0
if(j.eq.2)e=0.99d0
if(j.eq.3)e=0.999d0
if(j.eq.4)e=0.9999d0
if(j.eq.5)e=0.99999d0
if(j.eq.6)e=0.999999d0
if(j.eq.7)e=0.9999999d0

y1=1d0-e ! Initial values at x of y_1
y2=0d0 ! Initial values at x of y_2
y3=0d0 ! Initial values at x of y_1
y4=sqrt((1d0+e)/(1d0-e)) ! Initial values at x of y_2
tol=1d-4 ! Tolerance
Rtol(1:N)=tol ! Relative tolerance
Atol(1:N)=tol ! Absolute tolerance

! Imxplicit method
y(1)=y1; y(2)=y2; y(3)=y3; y(4)=y4; h=1d-6; count=0
istep=0
do i=2,Nx
info=0
tx=x(i-1)
do while(info.le.0)
call irkgl(istep,grk,N,tx,h,x(i),y,info,Atol,Rtol &
,ih,z0,ih0,h0,eta0,err0,Jup,ipiv,Jmat,epiv,errJ,iJac,Jac)
if(j.eq.4)write(21,'(7e25.15e3)')tx,y(1),y(2),y(3),y(4),h,&
0.5d0*(y(2)**2+y(4)**2) - 1d0/(sqrt(y(1)**2+y(3)**2))
enddo
h=h0
enddo
eneirk = 0.5d0*(y(2)**2+y(4)**2) - 1d0/(sqrt(y(1)**2+y(3)**2))
cirk=count
write(6,*)cirk

! Explicit method
y(1)=y1; y(2)=y2; y(3)=y3; y(4)=y4; h=1d-6; count=0
info=0
do i=2,Nx
info=-1
tx=x(i-1)
do while(info.le.0)
call drkf45(egrk,tx,h,N,y,x(i),info,tol)
if(j.eq.4)write(22,'(7e25.15e3)')tx,y(1),y(2),y(3),y(4),h, &
0.5d0*(y(2)**2+y(4)**2) - 1d0/(sqrt(y(1)**2+y(3)**2))
enddo
enddo
eneerk = 0.5d0*(y(2)**2+y(4)**2) - 1d0/(sqrt(y(1)**2+y(3)**2))
cerk = count
write(6,*)cerk

write(30,*)e,cirk,count,eneirk,eneerk

enddo

stop
end program main

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

! Write right-hand-side of 1st order ODEs

f(1)=y(2)
f(2)=-y(1)/((y(1)**2+y(3)**2)**1.5d0)
f(3)=y(4)
f(4)=-y(3)/((y(1)**2+y(3)**2)**1.5d0)

count=count+1
return
end subroutine grk

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

subroutine irkgl(istep,grk,Neq,x,h,xend,y,info,Atol,Rtol &
,ih,z0,ih0,h0,eta0,err0,Jup,ipiv,Jmat,epiv,errJ,iJac,Jac)
implicit none
integer,parameter::s=3
integer,intent(in)::Neq
integer,intent(inout)::istep,info,Jup,ipiv(1:Neq*s),epiv(1:Neq)
double precision,intent(in)::xend,Atol(1:Neq),Rtol(1:Neq)
double precision,intent(inout)::x,h,y(1:Neq)
integer,intent(inout)::ih0,ih,iJac
double precision,intent(inout)::z0(1:Neq*s),h0,eta0,err0
double precision,intent(inout)::Jmat(1:Neq*s,1:Neq*s),errJ(1:Neq,1:Neq)
double precision,intent(inout)::Jac(1:Neq,1:Neq)
external::grk
!
! Implicit Runge-Kutta method based on
! the Gauss-Legendre 3-stage 6-order
!
! Properties of this routine:
! 1. A-stable
! 2. Symplectic
! 3. Symmetric
! 4. Step size control
! Note, Gauss-Legendre IRK method is Symplectic
! even if we change the step size.
!
! Meaning of parameters
! istep : Number of IRK step
! grk : Right hand Side of ODEs
! Neq : Number of 1st-order ODEs
! x : Integral parameter (automatically updated)
! h : Step size (automatically updated)
! xend : End point of the x range
! y : Values of ODEs
! info : Information of the IRK process
! Atol : Absolute tolerance
! Rtol : Relative tolerance
!
! Other parameters are work parameters,
! referenced for istep >= 1
! ***0 : Previous parameters
! Jmat : LU decomposited J' (= I-hAJ) matrix
! ipiv : Pivot information for Jmat
! errJ : LU decomposited (I-h\gamma J) matrix
! epiv : Pivot information for errJ
! Jup : Decide the update of Jmat and errJ,
! Jup = 0 --> No update
! Jup = 1 --> Update
! Jac : Jacobian matrix of the grk function
! iJac : Decide the update of Jac,
! iJac = 0 --> No update when Jup=1
! iJac = 1 --> Update when Jup=1
!
! How to use?
! 1. Call this routine with istep = 0 and info = 0.
! 2. Call and Loop this routine while info <= 0.
!
! ==== Example =====
! istep=0
! info=0
! do while(info.le.0)
! call irk(istep,grk,Neq,x,h,xend,y,info,Atol,Rtol &
! ,ih,z0,ih0,h0,eta0,err0,Jup,ipiv,Jmat,epiv,errJ,iJac,Jac)
! enddo
! ==================
!
! After starting computation with istep=0,
! you must not touch WORK parameters.
!
!
! istep = 0 : when you start computation,
! set WORK parameters like ;
! ih = 0
! z0(1:Ns) = 0d0
! ih0 = 0
! h0 = 0d0
! eta0 = 0d0
! err0 = 0d0
! Jup = 0
! ipiv(1:Ns) = 0
! Jmat(1:Ns,1:Ns) = 0d0
! epiv(1:Neq) = 0
! errJ(1:Neq,1:Neq) = 0d0
! iJac = 0
! Jac(1:Neq,1:Neq) = 0d0
!
! Author : sikino
! URL : http://slpr.sakura.ne.jp/qp/
! Date : 2019/01/14 (yyyy/mm/dd)
! 2019/01/21 keep Jacobian matrix
! 2019/01/22 did Something
! 2019/01/23 initial value estimation
!
double precision::tx,th
double precision,allocatable::ty(:),tz0(:)

integer,parameter::kmax=7 ! Newton iteration limit
double precision,parameter::hmin=1d-13,hmax=1d0

integer::kexit,Ns,key,FLAG,Newt
double precision::theta,err,fac,del,del1,del2,teta0,th0

if(istep.le.-1)then
write(6,*)"**** Error, unexpected istep"
stop
endif

Ns=Neq*s
if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
endif
if(h.ge.abs(xend-x))h=xend-x
FLAG=1
if(abs(x-xend).le.hmin)then
info=1
FLAG=0
endif

if(istep.eq.0)then
ipiv=0; epiv=0; Jmat=0d0; errJ=0d0; Jup=0; ih=0
ih0=0; z0=0d0; eta0=0d0; h0=0d0; err0=0d0; iJac=0; Jac=0d0
call discrete_h(h,ih,th,hmin,hmax)
h=th
ih0 = ih+1
Jup = 1
iJac = 1
else
endif

key=0
allocate(ty(1:Neq),tz0(1:Ns))
ty=0d0; tz0=0d0
do while(FLAG.eq.1)

if(ih.ne.ih0)then
Jup=1
endif

tx=x
ty(1:Neq)=y(1:Neq)
tz0(1:Ns)=z0(1:Ns)
teta0 = eta0
th0=h0
call dirk6(istep,grk,Neq,tx,h,ty,Jup &
,tz0,th0,teta0,ipiv,Jmat,epiv,errJ,Atol,Rtol &
,err,kmax,kexit,Newt,theta,iJac,Jac)

! Even if the step is fail, ipiv,Jmat,epiv,errJ are updated if Jup=1.
if(Jup.eq.1)then
Jup=0
endif

if(kexit.eq.1)then
! Change step size
!fac = 0.9d0*(2d0*kmax+1d0)/(2d0*kmax+dble(Newt-1))
fac = 0.95d0*(2d0*kmax+1d0)/(2d0*kmax+dble(Newt))
if(err.ge.1d-30)then
if(istep.eq.0)then
del = fac*((1d0/err)**(0.33d0)) !(27)
else
del1 = fac*((1d0/err)**(0.33d0)) !(27)
del2 = del1*(h/h0)*((err0/err)**(0.33d0)) !(27)

del = del1
if(del2.lt.del)del=del2
if(del.gt.1d0)then
del=del1
else
del=del2
endif
endif
else
del=100d0
endif
elseif(kexit.eq.2)then
del=0d0
else
write(6,*)" **** detect unexpected kexit"
stop
endif

! Accept or Reject
if(err.le.1d0.or.key.eq.1)then
FLAG=0 ! This step with h is accepted
x=x+h
y(1:Neq)=ty(1:Neq)
z0(1:Ns)=tz0(1:Ns)
h0=h
ih0=ih
eta0=teta0
err0=err
iJac = 1
Jup=1
! Don't update Jacobian for next step
if(Newt.le.2.or.theta.lt.1d-3)Jup=0
!if(Newt.le.1.or.theta.lt.1d-3)Jup=0
endif

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

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 size alignment
if(abs(xend-x).le.abs(h))then
h=xend-x
Jup=1
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
else
call discrete_h(h,ih,th,hmin,hmax)
h=th
endif

if(h.le.0d0.and.xend-x.ge.0d0)then
info=1
FLAG=0
elseif(h.ge.0d0.and.xend-x.le.0d0)then
info=1
FLAG=0
endif
if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')"Strange point between ",x-h," and ",x
info=-9
endif
enddo

istep=istep+1

return
end subroutine irkgl

subroutine dirk6(istep,grk,Neq,x,h,y,Jup &
,z0,h0,eta0,ipiv,Jmat,epiv,errJ,Atol,Rtol &
,err,kmax,kexit,Newt,theta,iJac,Jac)
implicit none
integer,parameter::s=3
integer,intent(in)::istep,Neq,Jup,kmax
integer,intent(out)::kexit,Newt
double precision,intent(in)::h,Atol(1:Neq),Rtol(1:Neq)
double precision,intent(inout)::x,y(1:Neq),z0(1:Neq*s),h0,eta0
double precision,intent(out)::err,theta

integer,intent(inout)::ipiv(1:Neq*s),epiv(1:Neq),iJac
double precision,intent(inout)::Jmat(1:Neq*s,1:Neq*s),errJ(1:Neq,1:Neq)
double precision,intent(inout)::Jac(1:Neq,1:Neq)
external::grk
!
! istep >= 0
!
! Input
! Jup = 0 : Don't update Jacobian
! = 1 : Update Jacobian
! Output
! kexit = 1 : Newton iteration converge.
! = 2 : Newton iteration didn't converge.
!
! Below parameters have meaning when kexit = 1.
! Newt : Number of Newton iteration till converge.
! theta : Convergion rate, \theta < 1.
! err : Estimated err, if err < 1, satisfied required tol.
!
integer::i,j,k,n,m,p,q,Ns,info
double precision,parameter::kappa=5d-2
double precision,parameter::Uround=5d-16
double precision,parameter::sq15=sqrt(15d0)
! Real eigenvalue of A matrix in butcher table for Gauss-Legendre
double precision,parameter::gamma=0.215314423116112178244733530380696d0

double precision::a(1:s,1:s),c(1:s),d(1:s),d2(1:s),dc(1:s)
double precision::c12,c23,c31,xc1,xc2,xc3,xx,omega
double precision::x0,tx
double precision,allocatable::z(:),y0(:),dy(:),f(:),tf(:),ty(:),tf0(:)
double precision,allocatable::w(:,:),w0(:,:),e(:)
double precision::Ntol,sc,sdz,sdz0,tmp,eta

Ns=Neq*s
allocate(z(1:Ns))
z=0d0

! 3-stage Gauss-Legendre
c(1:3)=(/0.5d0-0.1d0*sq15, 0.5d0, 0.5d0+0.1d0*sq15/)
a(1,1:3)=(/5d0/36d0, 2d0/9d0-sq15/15d0, 5d0/36d0-sq15/30d0/)
a(2,1:3)=(/5d0/36d0+sq15/24d0, 2d0/9d0, 5d0/36d0-sq15/24d0/)
a(3,1:3)=(/5d0/36d0+sq15/30d0, 2d0/9d0+sq15/15d0, 5d0/36d0/)
d(1:3)=(/5d0/3d0,-4d0/3d0,5d0/3d0/)
d2(1:3)=(/-15d0-10d0*sqrt(5d0/3d0),12d0,-15d0+10d0*sqrt(5d0/3d0)/)
dc(1:3)=(/(5d0+sq15)*10d0/3d0, -40d0/3d0 ,(5d0-sq15)*10d0/3d0/)

if(istep.eq.0)then
z(1:Ns)=0d0
else
allocate(dy(1:Neq))
dy=0d0
do n=1,Neq
do j=1,s
dy(n)=dy(n)+d(j)*z0((n-1)*s+j)
enddo
enddo
! Lagrange interpolation
omega=h/h0
c12=1d0/(c(1)-c(2))
c23=1d0/(c(2)-c(3))
c31=1d0/(c(3)-c(1))
do n=1,Neq
do p=1,s
xx=1d0+omega*c(p)
xc1=xx-c(1)
xc2=xx-c(2)
xc3=xx-c(3)
z((n-1)*s+p)=&
-z0((n-1)*s+1)*xc2*xc3*(c12*c31)*xx/c(1) &
-z0((n-1)*s+2)*xc3*xc1*(c12*c23)*xx/c(2) &
-z0((n-1)*s+3)*xc1*xc2*(c31*c23)*xx/c(3) &
-dy(n)
enddo
enddo

deallocate(dy)
endif

if(Jup.eq.1)then
! LU decomposition of J' matrix
if(iJac.eq.1)then
call Jacobian(Neq,x,y,grk,Jac)
iJac=0
endif
do m=1,Neq
do n=1,Neq
do q=1,s
do p=1,s
Jmat(s*(n-1)+p,s*(m-1)+q) = - h*Jac(n,m)*a(p,q)
enddo
enddo
enddo
enddo
do n=1,Neq
do m=1,Neq
errJ(n,m) = - h*gamma*Jac(n,m)
enddo
enddo
do i=1,Ns
Jmat(i,i) = 1d0 + Jmat(i,i)
enddo
do i=1,Neq
errJ(i,i) = 1d0 + errJ(i,i)
enddo
! LU factorization for main part of IRK
call dgetrf(Ns,Ns,Jmat,Ns,ipiv,info)
! LU factorization for estimate error
call dgetrf(Neq,Neq,errJ,Neq,epiv,info)
endif

allocate(f(1:Ns),tf(1:Neq),ty(1:Neq),w(1:Ns,1:1))
f=0d0; tf=0d0; ty=0d0; w=0d0

!===============================
Ntol=sqrt(Rtol(1))
!if(0.03d0.le.Ntol)Ntol=0.03d0
!if(1d-12.le.Ntol)Ntol=1d-12
if(1d-12.le.Ntol)Ntol=1d-12
!if(1d-6.le.Ntol)Ntol=1d-6
!===============================

sdz0=0d0 ! Initialize
! Simple Newton iteration
kexit=0
do k=1,kmax

do j=1,s
tx = x + c(j)*h
do n=1,Neq
ty(n) = y(n) + z((n-1)*s+j)
enddo
call grk(Neq,tx,ty,tf)
do n=1,Neq
f((n-1)*s+j) = tf(n)
enddo
enddo

w(1:Ns,1:1)=0d0
do n=1,Neq
do p=1,s
do j=1,s
w((n-1)*s+p,1) = w((n-1)*s+p,1) + a(p,j)*f((n-1)*s+j)
enddo
w((n-1)*s+p,1) = z((n-1)*s+p) - h*w((n-1)*s+p,1)
enddo
enddo

! Solve J' \delta z = - w
do i=1,Ns
w(i,1) = - w(i,1)
enddo
call dgetrs('N',Ns,1,Jmat,Ns,ipiv,w,Ns,info)

! --> Now, w is \Delta z
! z^{(k+1)} = z^{(k)} + Delta z
do i=1,Ns
z(i) = z(i) + w(i,1)
enddo

sdz=0d0
do i=1,Ns
sdz = sdz + w(i,1)**2
enddo
sdz=sqrt(sdz)

! Stop iteration criteria
if(istep.eq.0.and.k.eq.1)then
! Do nothing
kexit=0
elseif(istep.eq.0.and.k.ge.2)then
if(sdz0-sdz.lt.0d0)then
! Convergion rate > 1, must change small step size.
kexit=2
elseif(sdz0*sdz**(kmax-k+1).gt.kappa*Ntol*(sdz0-sdz)*sdz0**(kmax-k))then
! Rough convergion estimation fail, must change small step size.
kexit=2
elseif(sdz*sdz.lt.kappa*Ntol*(sdz0-sdz))then
! Good behavior. Iteration finish.
kexit=1
endif
elseif(istep.ge.1.and.k.eq.1)then
tmp = Uround
if(eta0.gt.tmp)tmp=eta0
tmp=tmp**0.8d0
if(tmp*sdz.lt.kappa*Ntol)then
kexit=1
endif
elseif(istep.ge.1.and.k.ge.2)then
if(sdz0-sdz.lt.0d0)then
! Convergion rate > 1, must change small step size.
kexit=2
elseif(sdz0*sdz**(kmax-k+1).gt.kappa*Ntol*(sdz0-sdz)*sdz0**(kmax-k))then
! Rough convergion estimation fail, must change small step size.
kexit=2
elseif(sdz*sdz.lt.kappa*Ntol*(sdz0-sdz))then
! Good behavior. Iteration finish.
kexit=1
endif
else
write(6,*)" *****Unexpected parameters"
stop
endif

if(kexit.ne.0)exit
sdz0 = sdz
enddo

!
! 0 < \eta < \infty --> Good.
! if \eta ~ 0, good behavior
! -\infty < \eta < -1, --> Bad.
! Error increase as iteration increase.
!

if(k.eq.kmax+1.or.kexit.eq.2.or.kexit.eq.0)then
! Did not converge k_max iteration.
kexit=2
eta=100d0
Newt=kmax
err=100d0
theta=1000d0 ! Convergion ratio, \theta ~ 0 is Good.
! h0, z0, x, y are don't updated
return
else
! (sdz0-sdz) > 0
if(k.eq.1)then
! No sdz0 case.
eta = Uround
if(eta0.gt.eta)eta=eta0
eta = eta**0.8d0
theta = 100d0 ! Here, \theta cannot evaluate because k=1.
else
if(sdz0.eq.sdz)then
! z does not change --> converge enough.
eta = Uround
theta = 0d0 ! \theta=0.
else
! General case.
eta = sdz/(sdz0-sdz)
theta = sdz/sdz0
endif
endif
Newt = k
endif

deallocate(w)

x0=x
allocate(y0(1:Neq))
y0(1:Neq)=y(1:Neq)

! Update x and y(1:Neq)
x=x+h
do n=1,Neq
do j=1,s
y(n) = y(n) + d(j)*z((n-1)*s+j)
enddo
enddo

! Error estimation
allocate(w(1:Neq,1:1),w0(1:Neq,1:1),tf0(1:Neq),e(1:Neq))
w=0d0; w0=0d0; tf0=0d0; e=0d0
do n=1,Neq
do j=1,s
w(n,1) = w(n,1) + dc(j)*z((n-1)*s+j)
enddo
enddo
w0(1:Neq,1:1)=w(1:Neq,1:1)
call dgetrs('N',Neq,1,errJ,Neq,epiv,w,Neq,info)
do n=1,Neq
e(n)=w(n,1)
enddo

err=0d0
do n=1,Neq
sc=abs(y0(n))
if(abs(y(n)).gt.y0(n))sc=abs(y(n))
sc=Atol(n)+sc*Rtol(n)
err=err+(e(n)/sc)**2
enddo
err=sqrt(err/dble(Neq))
if(err.ge.1d0)then
call grk(Neq,x0,y0,tf0)
do n=1,Neq
ty(n) = y0(n) + e(n)
enddo
call grk(Neq,x0,ty,tf)
do n=1,Neq
w(n,1) = w0(n,1) + gamma*h*(tf(n)-tf0(n))
enddo
call dgetrs('N',Neq,1,errJ,Neq,epiv,w,Neq,info)
do n=1,Neq
e(n)=w(n,1)
enddo

err=0d0
do n=1,Neq
sc=abs(y0(n))
if(abs(y(n)).gt.y0(n))sc=abs(y(n))
sc=Atol(n)+sc*Rtol(n)
err=err+(e(n)/sc)**2
enddo
err=sqrt(err/dble(Neq))
endif

z0(1:Ns)=z(1:Ns)
h0=h
eta0=eta

return
end subroutine dirk6

subroutine discrete_h(h,ih,th,hmin,hmax)
implicit none
double precision,intent(in)::h,hmin,hmax
integer,intent(out)::ih
double precision,intent(out)::th

double precision::dx,hmin1,hmax1
integer,parameter::imax=200

hmin1=0.10d0*hmin
hmax1=1d0*hmax

dx=(log10(hmax1)-log10(hmin1))/dble(imax)

do ih=0,imax
th=10d0**(-ih*dx+log10(hmax1))
if(th.le.abs(h))exit
enddo

if(h.lt.0d0)th=-th

return
end subroutine discrete_h

subroutine Jacobian(Neq,x,y,fxy,Jac)
implicit none
integer,intent(in)::Neq
double precision,intent(in)::x,y(1:Neq)
double precision,intent(out)::Jac(1:Neq,1:Neq)
external::fxy

integer::n,m
double precision::dy
double precision,parameter::delta=2d-8
double precision,allocatable::f0(:),f1(:),ty(:)

allocate(f0(1:Neq),f1(1:Neq),ty(1:Neq))
f0=0d0; f1=0d0; ty=0d0

call fxy(Neq,x,y,f0)
do m=1,Neq
ty(1:Neq) = y(1:Neq)

dy=sqrt(abs(y(m)))
if(dy.lt.1d0)dy=1d0
dy=delta*dy

ty(m) = ty(m)+dy
call fxy(Neq,x,ty,f1)
do n=1,Neq
Jac(n,m) = (f1(n)-f0(n)) / dy
enddo
enddo

return
end subroutine Jacobian

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

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

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

egrk=0d0
if(s.eq.1)then
egrk = y(2)
elseif(s.eq.2)then
egrk = -y(1)/((y(1)**2+y(3)**2)**1.5d0)
elseif(s.eq.3)then
egrk = y(4)
elseif(s.eq.4)then
egrk = -y(3)/((y(1)**2+y(3)**2)**1.5d0)
else
write(6,*)"***Error grk"; stop
endif

if(s.eq.1)count=count+1
return
end function egrk

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

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
double precision,external::grk

integer::i,j,FLAG,key
double precision::R,delta,tx,Sy,err
double precision,allocatable::tmp(:),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(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
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=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/dble(N))/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
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(tmp,K)
return
end subroutine drkf45

きっかけ

きっかけとしてはtwitterで流れてきまして、やってみよう、と思いました。

RKFにもってこいの題材だよ。（みょんみょん度とかは無視してください（） pic.twitter.com/VWouH7IhWZ

— みーくん | itmz153 (@math153arclight) 2018年11月6日

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

陰的ルンゲ=クッタ法の高速化

2019年1月14日 sikino コメントする

陰的ルンゲ=クッタ法の基本的な考えは
陰的ルンゲ=クッタ法
にて説明しました。
基本的には
$
y(x+h)=y(x)+\Delta y
$
の右辺$y(x)+\Delta y$を計算します。
しかし、陰的ルンゲ=クッタ法の方法を数値計算で行おうとすると望ましくない部分が現れます。
それは、

$y(x)=O(h^0),~~\Delta y=O(h^1)$を同時に扱わなければならず、桁落ちが激しい
関数の評価回数が多い
ヤコビアン、LU分解の計算コストが非常に高い

という点です。
簡易ニュートン法を用いる事を前提にしておくと、上の問題は若干解決することが出来ます。

本稿では陰的ルンゲ=クッタ法を発展させ、数値計算的に陰的ルンゲ=クッタ法のアルゴリズムを工夫し、どのように計算量を減らすか？に焦点を当てます。

注意
[2,5]の同著者の新しい論文では、本稿の計算方法([1]に従う方法)ではなく、
もう一段階変換してからニュートン法を利用しています。
正直な所、私自身が追いきれなかっとのと、複素数が入ってくるプログラムでしたので、[1]の方法で止めておきます。
2つの計算方法の収束の早さや精度を比較を[2]で行っていますので、気になる方はそちらをご覧ください。

計算方法

解きたい問題は$N_{\text{eq}}$本の連立一次微分方程式

です。これを$s$段のルンゲ=クッタ法で解くことを考えます。
座標や添え字は以下のように決めました。

すると次のステップの値は

と求められます。$z_{np}^{[i]}$を求める事が問題になります。ここで、$d_j,~~(j=1,2,\cdots,s)$は既知で

のように、Butcherテーブルから求められます。
具体的に、3段6次の陰的公式であるGauss-Legendreの場合、

と求められます。
$z_{np}^{[i]}$の具体的な形は

で計算されます。$z_{np}^{[i]}$をベクトルとして表したいので、

と変換します。
この非線形方程式はニュートン法によって求める事が出来て、以下の通り、$k$回の繰り返しで解が収束するまで計算されます。

ここで、解が収束した結果が求めたい$z_{np}^{[i]}$となります。すなわち、

です。行列$J’$は

と求められます$(n,m=1,2,\cdots,N_{\text{eq}}, p,q=1,2,\cdots,s)$。ここで、$\delta_{i,j}$はクロネッカーのデルタ、$a_{i,j}$はButcherテーブルの値、$J$はヤコビアンで、

と計算されます。また、ベクトル$\mathbf{w}_{np}^{(k)}$は、

と定義します。

初期値の推定

初期値は$i=0$の初めのステップでは

として求め、それ以降では過去の結果を多項式補間して求めるとかなりよい精度です[1]。

ニュートン法の停止

ニュートン法は以下の条件が満たされた時、終了します。

まず、$i=0$の初めのステップでは必ず$k\ge 2$まで行い、
その後、以下の条件が満たされた時、終了します。

ここで、$\eta_k^{[i]}$は

で定義され、特に$\theta_k$は収束割合(convergion rate)と呼ばれます。
また、$\kappa=0.1\sim 0.01,~~Ntol=\text{min}(0.03,\sqrt{Rtol})$です[5]。ここで、$Rtol$は相対誤差を意味します。
実際に組んでみますと、理論が違うせいなのか、$Ntol=\text{min}(0.03,\sqrt{Rtol})$では十分なほど収束はしませんでした。なので、実際に組む時には$Ntol=10^{-9}$にしてしまっています。もしかしたら、$10^{-9}$では足らず、もっと必要かもしれません。

$i\ge 1$のステップでは$k=1$の時は前ステップの$\eta$の値を使い、

で判定します。ここで、倍精度演算ならば$Uround=5\times 10^{-16}$です。

$i\ge 1,~~k\ge 2$では

に従い、計算します。

ニュートン法の繰り返しの$k\ge 2$において、どこかで

が満たされる様であれば、刻み幅が大きすぎて収束しないことを意味します。なので、刻み幅を小さくする必要があります。

誤差判定

区間$i$の$n$番目の方程式の解の誤差は、以下の連立方程式を解いて得られます[1]。

ここで、$\gamma_0$はButcherテーブルの行列Aの実固有値で、ガウス=ルジャンドル陰的ルンゲクッタの3段6次であれば、
$
\gamma_0=0.215314423116112178244733530380696
$
を得ます。そして、上の方程式を解いた後に、計算結果を棄却するか判定するために、量

を計算します。
もしも$||\text{err}^{[i]}||\lt 1$であれば、その刻み幅で計算した結果を採用します。
しかし、$||\text{err}^{[i]}||\ge 1$であれば、以下の連立方程式を解きます([1]の”Hump”問題を参照)。

上を計算し、$||\text{err}^{[i]}||$を再計算した結果、もしも$||\text{err}^{[i]}||\lt 1$であれば、その刻み幅で計算した結果を採用し、そうでなければ刻み幅を次の節に従って変更します。

刻み幅制御

刻み幅制御をするためには, 2つ新しい刻み幅の推定値である$h_{i+1}^{(1)},h_{i+1}^{(2)}$を計算します。

その結果、刻み幅が小さくなるか大きくなるかに従って、どちらの刻み幅を採用するか決定します[2]。

プログラム

Fortran90で書いたプログラムはこちら。LU分解と連立方程式を解くため、LAPACKを使います。
それにリンクしてコンパイル、実行をしてください。
モジュールを使用していますが、これは関数の呼び出し回数を計測するためだけにグローバル変数として使っているため、消してもプログラムに何の影響もありません。

追記）
色々計算してみました。その結果、$tol=10^{-8}$より小さい値は使わない方が良さそうです。どうもこれ以上の精度にしてしまうと誤差の溜まり具合が増えてしまう感じがします。

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

module GBL
implicit none
integer::count
end module GBL

program main
use GBL
implicit none
integer::N,Ns,istep,ih0,Jup,info,ih,iJac
double precision::x,xend,h,h0,err0,eta0,tol
integer,allocatable::ipiv(:),epiv(:)
double precision,allocatable::y(:),Rtol(:),Atol(:)
double precision,allocatable::Jmat(:,:),errJ(:,:),z0(:),Jac(:,:)
external::grk

N=2 ! Number of 1st order ODEs
allocate(y(1:N),Rtol(1:N),Atol(1:N))
y=0d0; Rtol=0d0; Atol=0d0

x=0d0 ! Initial values of x
xend=2d0 ! End point of x
y(1)=2d0 ! Initial values at x of y_1
y(2)=-0.66d0 ! Initial values at x of y_2
tol=1d-4 ! Tolerance

h=1d-6 ! Initial step size
Rtol(1:N)=tol ! Relative tolerance
Atol(1:N)=tol ! Absolute tolerance

!------------ Initial set up ------------
Ns=N*3 ! Ns=N*s, s means s-stage IRK.
allocate(ipiv(1:Ns),epiv(1:N))
allocate(Jmat(1:Ns,1:Ns),errJ(1:N,1:N),z0(1:Ns),Jac(1:N,1:N))
ipiv=0; epiv=0; Jmat=0d0; errJ=0d0; Jup=0; ih=0
ih0=0; z0=0d0; eta0=0d0; h0=0d0; err0=0d0; iJac=0; Jac=0d0
!----------------------------------------
count=0
istep=0
info=0
do while(info.le.0)
call irkgl(istep,grk,N,x,h,xend,y,info,Atol,Rtol &
,ih,z0,ih0,h0,eta0,err0,Jup,ipiv,Jmat,epiv,errJ,iJac,Jac)
write(10,'(4e25.15e3,1i5)')x,y(1),y(2),h
enddo

write(6,*)count
stop
end program main

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

! Write right-hand-side of 1st order ODEs

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

count=count+1
return
end subroutine grk

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

subroutine irkgl(istep,grk,Neq,x,h,xend,y,info,Atol,Rtol &
,ih,z0,ih0,h0,eta0,err0,Jup,ipiv,Jmat,epiv,errJ,iJac,Jac)
implicit none
integer,parameter::s=3
integer,intent(in)::Neq
integer,intent(inout)::istep,info,Jup,ipiv(1:Neq*s),epiv(1:Neq)
double precision,intent(in)::xend,Atol(1:Neq),Rtol(1:Neq)
double precision,intent(inout)::x,h,y(1:Neq)
integer,intent(inout)::ih0,ih,iJac
double precision,intent(inout)::z0(1:Neq*s),h0,eta0,err0
double precision,intent(inout)::Jmat(1:Neq*s,1:Neq*s),errJ(1:Neq,1:Neq)
double precision,intent(inout)::Jac(1:Neq,1:Neq)
external::grk
!
! Implicit Runge-Kutta method based on
! the Gauss-Legendre 3-stage 6-order
!
! Properties of this routine:
! 1. A-stable
! 2. Symplectic
! 3. Symmetric
! 4. Step size control
! Note, Gauss-Legendre IRK method is Symplectic
! even if we change the step size.
!
! Meaning of parameters
! istep : Number of IRK step
! grk : Right hand Side of ODEs
! Neq : Number of 1st-order ODEs
! x : Integral parameter (automatically updated)
! h : Step size (automatically updated)
! xend : End point of the x range
! y : Values of ODEs
! info : Information of the IRK process
! Atol : Absolute tolerance
! Rtol : Relative tolerance
!
! Other parameters are work parameters,
! referenced for istep >= 1
! ***0 : Previous parameters
! Jmat : LU decomposited J' (= I-hAJ) matrix
! ipiv : Pivot information for Jmat
! errJ : LU decomposited (I-h\gamma J) matrix
! epiv : Pivot information for errJ
! Jup : Decide the update of Jmat and errJ,
! Jup = 0 --> No update
! Jup = 1 --> Update
! Jac : Jacobian matrix of the grk function
! iJac : Decide the update of Jac,
! iJac = 0 --> No update when Jup=1
! iJac = 1 --> Update when Jup=1
!
! How to use?
! 1. Call this routine with istep = 0 and info = 0.
! 2. Call and Loop this routine while info <= 0.
!
! ==== Example =====
! istep=0
! info=0
! do while(info.le.0)
! call irk(istep,grk,Neq,x,h,xend,y,info,Atol,Rtol &
! ,ih,z0,ih0,h0,eta0,err0,Jup,ipiv,Jmat,epiv,errJ,iJac,Jac)
! enddo
! ==================
!
! After starting computation with istep=0,
! you must not touch WORK parameters.
!
!
! istep = 0 : when you start computation,
! set WORK parameters like ;
! ih = 0
! z0(1:Ns) = 0d0
! ih0 = 0
! h0 = 0d0
! eta0 = 0d0
! err0 = 0d0
! Jup = 0
! ipiv(1:Ns) = 0
! Jmat(1:Ns,1:Ns) = 0d0
! epiv(1:Neq) = 0
! errJ(1:Neq,1:Neq) = 0d0
! iJac = 0
! Jac(1:Neq,1:Neq) = 0d0
!
! Author : sikino
! URL : http://slpr.sakura.ne.jp/qp/
! Date : 2019/01/14 (yyyy/mm/dd)
! 2019/01/21 keep Jacobian matrix
! 2019/01/22 did Something
! 2019/01/23 initial value estimation
!
double precision::tx,th
double precision,allocatable::ty(:),tz0(:)

integer,parameter::kmax=7 ! Newton iteration limit
double precision,parameter::hmin=1d-13,hmax=1d0

integer::kexit,Ns,key,FLAG,Newt
double precision::theta,err,fac,del,del1,del2,teta0,th0

if(istep.le.-1)then
write(6,*)"**** Error, unexpected istep"
stop
endif

Ns=Neq*s
if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
endif
if(h.ge.abs(xend-x))h=xend-x
FLAG=1
if(abs(x-xend).le.hmin)then
info=1
FLAG=0
endif

if(istep.eq.0)then
ipiv=0; epiv=0; Jmat=0d0; errJ=0d0; Jup=0; ih=0
ih0=0; z0=0d0; eta0=0d0; h0=0d0; err0=0d0; iJac=0; Jac=0d0
call discrete_h(h,ih,th,hmin,hmax)
h=th
ih0 = ih+1
Jup = 1
iJac = 1
else
endif

key=0
allocate(ty(1:Neq),tz0(1:Ns))
ty=0d0; tz0=0d0
do while(FLAG.eq.1)

if(ih.ne.ih0)then
Jup=1
endif

tx=x
ty(1:Neq)=y(1:Neq)
tz0(1:Ns)=z0(1:Ns)
teta0 = eta0
th0=h0
call dirk6(istep,grk,Neq,tx,h,ty,Jup &
,tz0,th0,teta0,ipiv,Jmat,epiv,errJ,Atol,Rtol &
,err,kmax,kexit,Newt,theta,iJac,Jac)

! Even if the step is fail, ipiv,Jmat,epiv,errJ are updated if Jup=1.
if(Jup.eq.1)then
Jup=0
endif

if(kexit.eq.1)then
! Change step size
!fac = 0.9d0*(2d0*kmax+1d0)/(2d0*kmax+dble(Newt-1))
fac = 0.95d0*(2d0*kmax+1d0)/(2d0*kmax+dble(Newt))
if(err.ge.1d-30)then
if(istep.eq.0)then
del = fac*((1d0/err)**(0.33d0)) !(27)
else
del1 = fac*((1d0/err)**(0.33d0)) !(27)
del2 = del1*(h/h0)*((err0/err)**(0.33d0)) !(27)

del = del1
if(del2.lt.del)del=del2
if(del.gt.1d0)then
del=del1
else
del=del2
endif
endif
else
del=100d0
endif
elseif(kexit.eq.2)then
del=0d0
else
write(6,*)" **** detect unexpected kexit"
stop
endif

! Accept or Reject
if(err.le.1d0.or.key.eq.1)then
FLAG=0 ! This step with h is accepted
x=x+h
y(1:Neq)=ty(1:Neq)
z0(1:Ns)=tz0(1:Ns)
h0=h
ih0=ih
eta0=teta0
err0=err
iJac = 1
Jup=1
! Don't update Jacobian for next step
if(Newt.le.2.or.theta.lt.1d-3)Jup=0
!if(Newt.le.1.or.theta.lt.1d-3)Jup=0
endif

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

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 size alignment
if(abs(xend-x).le.abs(h))then
h=xend-x
Jup=1
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
else
call discrete_h(h,ih,th,hmin,hmax)
h=th
endif

if(h.le.0d0.and.xend-x.ge.0d0)then
info=1
FLAG=0
elseif(h.ge.0d0.and.xend-x.le.0d0)then
info=1
FLAG=0
endif
if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')"Strange point between ",x-h," and ",x
info=-9
endif
enddo

istep=istep+1

return
end subroutine irkgl

subroutine dirk6(istep,grk,Neq,x,h,y,Jup &
,z0,h0,eta0,ipiv,Jmat,epiv,errJ,Atol,Rtol &
,err,kmax,kexit,Newt,theta,iJac,Jac)
implicit none
integer,parameter::s=3
integer,intent(in)::istep,Neq,Jup,kmax
integer,intent(out)::kexit,Newt
double precision,intent(in)::h,Atol(1:Neq),Rtol(1:Neq)
double precision,intent(inout)::x,y(1:Neq),z0(1:Neq*s),h0,eta0
double precision,intent(out)::err,theta

integer,intent(inout)::ipiv(1:Neq*s),epiv(1:Neq),iJac
double precision,intent(inout)::Jmat(1:Neq*s,1:Neq*s),errJ(1:Neq,1:Neq)
double precision,intent(inout)::Jac(1:Neq,1:Neq)
external::grk
!
! istep >= 0
!
! Input
! Jup = 0 : Don't update Jacobian
! = 1 : Update Jacobian
! Output
! kexit = 1 : Newton iteration converge.
! = 2 : Newton iteration didn't converge.
!
! Below parameters have meaning when kexit = 1.
! Newt : Number of Newton iteration till converge.
! theta : Convergion rate, \theta < 1.
! err : Estimated err, if err < 1, satisfied required tol.
!
integer::i,j,k,n,m,p,q,Ns,info
double precision,parameter::kappa=5d-2
double precision,parameter::Uround=5d-16
double precision,parameter::sq15=sqrt(15d0)
! Real eigenvalue of A matrix in butcher table for Gauss-Legendre
double precision,parameter::gamma=0.215314423116112178244733530380696d0

double precision::a(1:s,1:s),c(1:s),d(1:s),d2(1:s),dc(1:s)
double precision::c12,c23,c31,xc1,xc2,xc3,xx,omega
double precision::x0,tx
double precision,allocatable::z(:),y0(:),dy(:),f(:),tf(:),ty(:),tf0(:)
double precision,allocatable::w(:,:),w0(:,:),e(:)
double precision::Ntol,sc,sdz,sdz0,tmp,eta

Ns=Neq*s
allocate(z(1:Ns))
z=0d0

! 3-stage Gauss-Legendre
c(1:3)=(/0.5d0-0.1d0*sq15, 0.5d0, 0.5d0+0.1d0*sq15/)
a(1,1:3)=(/5d0/36d0, 2d0/9d0-sq15/15d0, 5d0/36d0-sq15/30d0/)
a(2,1:3)=(/5d0/36d0+sq15/24d0, 2d0/9d0, 5d0/36d0-sq15/24d0/)
a(3,1:3)=(/5d0/36d0+sq15/30d0, 2d0/9d0+sq15/15d0, 5d0/36d0/)
d(1:3)=(/5d0/3d0,-4d0/3d0,5d0/3d0/)
d2(1:3)=(/-15d0-10d0*sqrt(5d0/3d0),12d0,-15d0+10d0*sqrt(5d0/3d0)/)
dc(1:3)=(/(5d0+sq15)*10d0/3d0, -40d0/3d0 ,(5d0-sq15)*10d0/3d0/)

if(istep.eq.0)then
z(1:Ns)=0d0
else
allocate(dy(1:Neq))
dy=0d0
do n=1,Neq
do j=1,s
dy(n)=dy(n)+d(j)*z0((n-1)*s+j)
enddo
enddo
! Lagrange interpolation
omega=h/h0
c12=1d0/(c(1)-c(2))
c23=1d0/(c(2)-c(3))
c31=1d0/(c(3)-c(1))
do n=1,Neq
do p=1,s
xx=1d0+omega*c(p)
xc1=xx-c(1)
xc2=xx-c(2)
xc3=xx-c(3)
z((n-1)*s+p)=&
-z0((n-1)*s+1)*xc2*xc3*(c12*c31)*xx/c(1) &
-z0((n-1)*s+2)*xc3*xc1*(c12*c23)*xx/c(2) &
-z0((n-1)*s+3)*xc1*xc2*(c31*c23)*xx/c(3) &
-dy(n)
enddo
enddo

deallocate(dy)
endif

if(Jup.eq.1)then
! LU decomposition of J' matrix
if(iJac.eq.1)then
call Jacobian(Neq,x,y,grk,Jac)
iJac=0
endif
do m=1,Neq
do n=1,Neq
do q=1,s
do p=1,s
Jmat(s*(n-1)+p,s*(m-1)+q) = - h*Jac(n,m)*a(p,q)
enddo
enddo
enddo
enddo
do n=1,Neq
do m=1,Neq
errJ(n,m) = - h*gamma*Jac(n,m)
enddo
enddo
do i=1,Ns
Jmat(i,i) = 1d0 + Jmat(i,i)
enddo
do i=1,Neq
errJ(i,i) = 1d0 + errJ(i,i)
enddo
! LU factorization for main part of IRK
call dgetrf(Ns,Ns,Jmat,Ns,ipiv,info)
! LU factorization for estimate error
call dgetrf(Neq,Neq,errJ,Neq,epiv,info)
endif

allocate(f(1:Ns),tf(1:Neq),ty(1:Neq),w(1:Ns,1:1))
f=0d0; tf=0d0; ty=0d0; w=0d0

!===============================
Ntol=sqrt(Rtol(1))
!if(0.03d0.le.Ntol)Ntol=0.03d0
!if(1d-12.le.Ntol)Ntol=1d-12
if(1d-12.le.Ntol)Ntol=1d-12
!if(1d-6.le.Ntol)Ntol=1d-6
!===============================

sdz0=0d0 ! Initialize
! Simple Newton iteration
kexit=0
do k=1,kmax

do j=1,s
tx = x + c(j)*h
do n=1,Neq
ty(n) = y(n) + z((n-1)*s+j)
enddo
call grk(Neq,tx,ty,tf)
do n=1,Neq
f((n-1)*s+j) = tf(n)
enddo
enddo

w(1:Ns,1:1)=0d0
do n=1,Neq
do p=1,s
do j=1,s
w((n-1)*s+p,1) = w((n-1)*s+p,1) + a(p,j)*f((n-1)*s+j)
enddo
w((n-1)*s+p,1) = z((n-1)*s+p) - h*w((n-1)*s+p,1)
enddo
enddo

! Solve J' \delta z = - w
do i=1,Ns
w(i,1) = - w(i,1)
enddo
call dgetrs('N',Ns,1,Jmat,Ns,ipiv,w,Ns,info)

! --> Now, w is \Delta z
! z^{(k+1)} = z^{(k)} + Delta z
do i=1,Ns
z(i) = z(i) + w(i,1)
enddo

sdz=0d0
do i=1,Ns
sdz = sdz + w(i,1)**2
enddo
sdz=sqrt(sdz)

! Stop iteration criteria
if(istep.eq.0.and.k.eq.1)then
! Do nothing
kexit=0
elseif(istep.eq.0.and.k.ge.2)then
if(sdz0-sdz.lt.0d0)then
! Convergion rate > 1, must change small step size.
kexit=2
elseif(sdz0*sdz**(kmax-k+1).gt.kappa*Ntol*(sdz0-sdz)*sdz0**(kmax-k))then
! Rough convergion estimation fail, must change small step size.
kexit=2
elseif(sdz*sdz.lt.kappa*Ntol*(sdz0-sdz))then
! Good behavior. Iteration finish.
kexit=1
endif
elseif(istep.ge.1.and.k.eq.1)then
tmp = Uround
if(eta0.gt.tmp)tmp=eta0
tmp=tmp**0.8d0
if(tmp*sdz.lt.kappa*Ntol)then
kexit=1
endif
elseif(istep.ge.1.and.k.ge.2)then
if(sdz0-sdz.lt.0d0)then
! Convergion rate > 1, must change small step size.
kexit=2
elseif(sdz0*sdz**(kmax-k+1).gt.kappa*Ntol*(sdz0-sdz)*sdz0**(kmax-k))then
! Rough convergion estimation fail, must change small step size.
kexit=2
elseif(sdz*sdz.lt.kappa*Ntol*(sdz0-sdz))then
! Good behavior. Iteration finish.
kexit=1
endif
else
write(6,*)" *****Unexpected parameters"
stop
endif

if(kexit.ne.0)exit
sdz0 = sdz
enddo

!
! 0 < \eta < \infty --> Good.
! if \eta ~ 0, good behavior
! -\infty < \eta < -1, --> Bad.
! Error increase as iteration increase.
!

if(k.eq.kmax+1.or.kexit.eq.2.or.kexit.eq.0)then
! Did not converge k_max iteration.
kexit=2
eta=100d0
Newt=kmax
err=100d0
theta=1000d0 ! Convergion ratio, \theta ~ 0 is Good.
! h0, z0, x, y are don't updated
return
else
! (sdz0-sdz) > 0
if(k.eq.1)then
! No sdz0 case.
eta = Uround
if(eta0.gt.eta)eta=eta0
eta = eta**0.8d0
theta = 100d0 ! Here, \theta cannot evaluate because k=1.
else
if(sdz0.eq.sdz)then
! z does not change --> converge enough.
eta = Uround
theta = 0d0 ! \theta=0.
else
! General case.
eta = sdz/(sdz0-sdz)
theta = sdz/sdz0
endif
endif
Newt = k
endif

deallocate(w)

x0=x
allocate(y0(1:Neq))
y0(1:Neq)=y(1:Neq)

! Update x and y(1:Neq)
x=x+h
do n=1,Neq
do j=1,s
y(n) = y(n) + d(j)*z((n-1)*s+j)
enddo
enddo

! Error estimation
allocate(w(1:Neq,1:1),w0(1:Neq,1:1),tf0(1:Neq),e(1:Neq))
w=0d0; w0=0d0; tf0=0d0; e=0d0
do n=1,Neq
do j=1,s
w(n,1) = w(n,1) + dc(j)*z((n-1)*s+j)
enddo
enddo
w0(1:Neq,1:1)=w(1:Neq,1:1)
call dgetrs('N',Neq,1,errJ,Neq,epiv,w,Neq,info)
do n=1,Neq
e(n)=w(n,1)
enddo

err=0d0
do n=1,Neq
sc=abs(y0(n))
if(abs(y(n)).gt.y0(n))sc=abs(y(n))
sc=Atol(n)+sc*Rtol(n)
err=err+(e(n)/sc)**2
enddo
err=sqrt(err/dble(Neq))
if(err.ge.1d0)then
call grk(Neq,x0,y0,tf0)
do n=1,Neq
ty(n) = y0(n) + e(n)
enddo
call grk(Neq,x0,ty,tf)
do n=1,Neq
w(n,1) = w0(n,1) + gamma*h*(tf(n)-tf0(n))
enddo
call dgetrs('N',Neq,1,errJ,Neq,epiv,w,Neq,info)
do n=1,Neq
e(n)=w(n,1)
enddo

err=0d0
do n=1,Neq
sc=abs(y0(n))
if(abs(y(n)).gt.y0(n))sc=abs(y(n))
sc=Atol(n)+sc*Rtol(n)
err=err+(e(n)/sc)**2
enddo
err=sqrt(err/dble(Neq))
endif

z0(1:Ns)=z(1:Ns)
h0=h
eta0=eta

return
end subroutine dirk6

subroutine discrete_h(h,ih,th,hmin,hmax)
implicit none
double precision,intent(in)::h,hmin,hmax
integer,intent(out)::ih
double precision,intent(out)::th

double precision::dx,hmin1,hmax1
integer,parameter::imax=200

hmin1=0.10d0*hmin
hmax1=1d0*hmax

dx=(log10(hmax1)-log10(hmin1))/dble(imax)

do ih=0,imax
th=10d0**(-ih*dx+log10(hmax1))
if(th.le.abs(h))exit
enddo

if(h.lt.0d0)th=-th

return
end subroutine discrete_h

subroutine Jacobian(Neq,x,y,fxy,Jac)
implicit none
integer,intent(in)::Neq
double precision,intent(in)::x,y(1:Neq)
double precision,intent(out)::Jac(1:Neq,1:Neq)
external::fxy

integer::n,m
double precision::dy
double precision,parameter::delta=2d-8
double precision,allocatable::f0(:),f1(:),ty(:)

allocate(f0(1:Neq),f1(1:Neq),ty(1:Neq))
f0=0d0; f1=0d0; ty=0d0

call fxy(Neq,x,y,f0)
do m=1,Neq
ty(1:Neq) = y(1:Neq)

dy=sqrt(abs(y(m)))
if(dy.lt.1d0)dy=1d0
dy=delta*dy

ty(m) = ty(m)+dy
call fxy(Neq,x,ty,f1)
do n=1,Neq
Jac(n,m) = (f1(n)-f0(n)) / dy
enddo
enddo

return
end subroutine Jacobian

ある配列で定義したグリッド上の値

計算結果を、ある配列で定義したグリッド上で欲しい場合、以下のプログラムで行うことが出来ます。
下に載せたプログラムと、上のコードの中にあるサブルーチン(grk, irkgl, dirk6, discrete_h, Jacobian)を一緒にコンパイルしてください。

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

終点だけの結果が欲しい場合

終点だけの結果が欲しい場合、不要なwork配列などは省略できます。
下に載せたプログラムと、上のコードの中にあるサブルーチン(grk, irkgl, dirk6, discrete_h, Jacobian)を一緒にコンパイルしてください。

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

収束判定の余地

上のプログラムは収束判定を[1]と同じにしているため、過剰に評価しているパラメータになっているかもしれません。

その余地としては、計算回数を減らすために重要な順に

ニュートン法の収束判定 … Ntol
刻み幅の安全係数 … fac
刻み幅の離散化 … discrete_h内のimax
ヤコビアンの更新条件 … if(Newt.le.2.or.theta.lt.1d-3)Jup=0の箇所

です。現状のプログラムでは安全のために過剰評価気味にしています。

4倍精度とGECP

lapackを使わない場合、GECPというプログラムを使うことが出来ます。これは、
九州大学の渡部善隆様が公開なさっているGECP(Gaussian Elimination with Complete Pivoting, 一般実行列に対する連立1次方程式の数値解を完全ピボット選択付き Gauss の消去法によって求める Fortran サブルーチン)
で計算します。
http://yebisu.cc.kyushu-u.ac.jp/~watanabe/RESERCH/GECP/index.html。
再配布可能なので、上記プログラムで必要な物には組み込んであります。

また、倍精度、4倍精度に対応しているので、それらのプログラムを使って
陰的ルンゲクッタ法を書いたものが以下のものです。
irkgl_dge.f90
irkgl_qge.f90

倍精度 d
４倍精度 q

プログラムの評価

５つの問題について、評価した結果を上に載せます。評価の良し悪しは、
連立方程式の右辺の関数が何回評価されたかによって決めました。
比較対象は

陰的解法である本稿のプログラム（自作）
陽的解法であるルンゲ=クッタ=フェールベルグの4,5次（自作）
陰的解法であるRADAU IIAに基づくプログラム[3]

です。横軸に要求した精度、縦軸に実際に評価された回数を載せました。
まず、自作の陽的解法、陰的解法を比べますと、硬くない方程式であるEq.1,2,3では
陽的解法の方が10倍近く早いです。
しかし、硬い方程式であるEq.4,5は10倍から1000倍ほど陰的公式の方が早いという結果が得られました。望み通り、陰的解法が動いていることが確認できます。

…さて、専門家が書いたRADAU5ですが、めちゃくちゃ早いです。硬い方程式であるEq.4,5でも、自作したやつの1/10の回数で大体終わっています。
しかも、硬くない方程式であるEq.1,2,3ですら、自作の陽的解法よりも少ない回数で計算を終えています。本当にどういう事なんでしょうね…。
上には上がいるものです。

結論

ちゃんとした陰的解法が欲しいのであれば、自作せず、専門家のプログラムを使いましょう。

RADAU5の使い方

RADAU5のFORTRANプログラムはhttps://www.unige.ch/~hairer/software.htmlにあります。
使い方に手間取ったので、どうやって使うのかメモしておきます。

http://www.unige.ch/~hairer/testset/testset.htmlに移動し、Van der Pol方程式のメインプログラムをダウンロードする。場所は
・VDPOL　…　driver for RADAU5,
と書かれている所をクリックするとhttp://www.unige.ch/~hairer/testset/stiff/vdpol/driver_radau5.fに飛ぶ。これを保存し、driver_radau5.fという名前で保存する。
https://www.unige.ch/~hairer/software.htmlに移動
Stiff Differential Equations and Differential-Algebraic Problems
という項目の
・RADAU5
・DC_DECSOL
・DECSOL
のリンク先のプログラムをダウンロード。それぞれradau5.f, dc_decsol.f, decsol.fという名前で保存する。
合計4つのプログラムをダウンロードしたら、コンパイルを
gfortran decsol.f dc_decsol.f radau5.f driver_radau5.f

でコンパイルし、実行する。

参考文献

[1]E. Hairer and G. Wanner, ‘Solving Ordinary Differential Equations II’ Springer, 1996

[2]E.~Hairer and G.~Wanner. Stiff diferential equations solved by Radau methods, J. Comput. Appl. Math., 111:93-111, 1999.

[3]E. Hairer, Fortran and Matlab Codes https://www.unige.ch/~hairer/software.html

[4]10. 常微分方程式 (2)https://www.ktech.biz/jp/archives/1003, K Technologies Sites (2014)

[5]Nicola Guglielmi and E.~Hairer, User’s Guide for code RADAU5 – Version 2.1 (packed in “radar5-v2.1.tar”) http://www.unige.ch/~hairer/software.html, 2005