分数階微積分の数値計算

Loverroの "Fractional Calculus: History, Definitions and Applications for the Engineer" に沿って、数値計算がしやすそうなGrünwald–Letnikovの分数階微積分について見ていきます。

分数階微積分とは

ここではざっくりとした導入だけを行います。基本的な理論については次の資料が参考になります。

分数階微分について考えます。微分演算子 DD を導入します。

Df(x)=ddxf(x) D f(x) = \frac{d}{dx} f(x)

例えば 1/21 / 2 階微分を 22 回繰り返すと1階微分になります。

D12(D12f(x))=Df(x) D^{\frac{1}{2}} \left( D^{\frac{1}{2}} f(x) \right) = D f(x)

1/31 / 3 階微分を 33 回繰り返すと1階微分になります。

D13(D13(D13f(x)))=ddxf(x) D^{\frac{1}{3}} \left( D^{\frac{1}{3}} \left( D^{\frac{1}{3}} f(x) \right) \right) = \frac{d}{dx} f(x)

つまり 1/q1 / q 階微分を qq 回繰り返すと1階微分になります。

D1q((D1qf(x)))1/q differentiation for q times=ddxf(x) \underbrace{ D^{\frac{1}{q}} \left( \ldots \left( D^{\frac{1}{q}} f(x) \right) \right) }_\text{1/q differentiation for q times} = \frac{d}{dx} f(x)

分数階積分も同様に捉えることができます。積分演算子 JJ を定義します。

Jf(x)=0xf(s)ds J f(x) = \int_{0}^{x} f(s) ds

1/q1 / q 階積分を qq 回繰り返すと1階積分になります。

J1q((J1qf(x)))1/q integration for q times=0xf(s)ds \underbrace{ J^{\frac{1}{q}} \left( \ldots \left( J^{\frac{1}{q}} f(x) \right) \right) }_\text{1/q integration for q times} = \int_{0}^{x} f(s) ds

分数階微積分が物理的にどういう意味を持つのかはよくわかりませんが、波動方程式への応用が提案されています。

Grünwald–Letnikovの分数階微分

nn 階微分の差分形です。

dnf(x)=limh01hnm=0n(1)m(nm)f(xmh),nN d^n f(x) = \lim_{h \to 0} \frac{1}{h^n} \sum_{m=0}^{n} (-1)^{m} \binom{n}{m} f(x - mh) \,,\qquad n \in \Bbb{N}

m>nm > n のとき (nm)=0\binom{n}{m} = 0 なので総和の上限を無限にできます。

dnf(x)=limh01hnm=0(1)m(nm)f(xmh),nN d^n f(x) = \lim_{h \to 0} \frac{1}{h^n} \sum_{m=0}^{\infty} (-1)^{m} \binom{n}{m} f(x - mh) \,,\qquad n \in \Bbb{N}

自然数 nn を実数 α\alpha に置き換えます。

dαf(x)=limh01hαm=0(1)m(αm)f(xmh),αR d^\alpha f(x) = \lim_{h \to 0} \frac{1}{h^\alpha} \sum_{m=0}^{\infty} (-1)^{m} \binom{\alpha}{m} f(x - mh) \,,\qquad \alpha \in \Bbb{R}

この形がGrünwald–Letnikovの分数階微分です。

二項係数の計算

実数 α\alpha が入力される二項係数 (αm)\binom{\alpha}{m} を計算できる形にします。まずは自然数についての二項係数を展開します。

(nm)=n!m!(nm)!,nN,mN \binom{n}{m} = \frac{n!}{m!(n - m)!} \,,\qquad n \in \Bbb{N}\,,\quad m \in \Bbb{N}

階乗はガンマ関数 Γ(n)\Gamma(n) に置き換えられます。

Γ(n)=(n1)!,nN \Gamma(n) = (n - 1)!, \qquad n \in \Bbb{N}

ガンマ関数を用いて二項係数を変形します。

(αm)=Γ(α+1)Γ(m+1)Γ(αm+1),αR,mR \binom{\alpha}{m} = \frac{\Gamma(\alpha + 1)}{\Gamma(m + 1)\Gamma(\alpha - m + 1)} \,,\qquad \alpha \in \Bbb{R}\,,\quad m \in \Bbb{R}

Grünwald–Letnikovの分数階微分に代入します。

dαf(x)=limh01hαm=0(1)mΓ(α+1)Γ(m+1)Γ(αm+1)f(xmh),αR d^{\alpha} f(x) = \lim_{h \to 0} \frac{1}{h^{\alpha}} \sum_{m=0}^{\infty} (-1)^{m} \frac{\Gamma(\alpha + 1)}{\Gamma(m + 1)\Gamma(\alpha - m + 1)} f(x - mh) \,, \qquad \alpha \in \Bbb{R}

これで実数 α\alpha が入力されても計算できる形になりました。

ガンマ関数に入力される値が 00 あるいは負の整数のときは無限大になるので、計算できない場合があります。Grünwald–Letnikovの分数階微分の式では α>0\alpha > 0 かつ m>0m > 0 かつ α>m1\alpha > m - 1 のときは問題なく計算できます。

次のガンマ関数の計算と組み合わせた二項係数の計算はあまり精度が良くありません。できれば信頼できるライブラリを使ってください。

ガンマ関数の計算

ガンマ関数の数値計算にはLanczos approximationが使えます。次のコードはJavaScriptでの実装です。

function gamma(z) { var p = [ 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 ] if (z < 0.5) { return Math.PI / (Math.sin(Math.PI * z) * gamma(1 - z)) } z -= 1 var a = 0.99999999999980993 for (var i = 0; i < p.length; ++i) { a += p[i] / (z + i + 1) } var t = z + p.length - 0.5 return Math.sqrt(2 * Math.PI) * Math.pow(t, z + 0.5) * Math.exp(-t) * a }

Grünwald–Letnikovの分数階積分

α\alpha 階の分数階積分も計算できます。 α-\alpha 階微分の形ではガンマ関数に負の値が入ることになるので計算できない場合が出てきます。

dαf(x)=limh01hαm=0(1)m(αm)f(xmh)=limh01hαm=0(1)mΓ(α+1)Γ(m+1)Γ(αm+1)f(xmh) \begin{aligned} d^{-\alpha} f(x) &= \lim_{h \to 0} \frac{1}{h^{-\alpha}} \sum_{m=0}^{\infty} (-1)^{m} \binom{-\alpha}{m} f(x - mh)\\ &= \lim_{h \to 0} \frac{1}{h^{-\alpha}} \sum_{m=0}^{\infty} (-1)^{m} \frac{\Gamma(-\alpha + 1)}{\Gamma(m + 1)\Gamma(-\alpha - m + 1)} f(x - mh) \end{aligned}

負の整数について一般化された二項係数を使って式を変形します。

負の整数について一般化された二項係数 (nm)\binom{-n}{m} を展開します。

(nm)=n(n1)(n2)(n3)(nm+1)m!=(1)mn(n+1)(n+2)(n+3)(n+m1)m!=(1)m(n+m1)!m!(n1)!=(1)mΓ(n+m)Γ(m+1)Γ(n) \begin{aligned} \binom{-n}{m} &= \frac{-n(-n - 1)(-n - 2)(-n - 3) \ldots (-n - m + 1)}{m!}\\ &= (-1)^{m} \frac{n(n + 1)(n + 2)(n + 3) \ldots (n + m - 1)}{m!}\\ &= (-1)^{m} \frac{(n + m - 1)!}{m!(n - 1)!}\\ &= (-1)^{m} \frac{\Gamma(n + m)}{\Gamma(m + 1)\Gamma(n)} \end{aligned}

これで α\alpha 階積分が定義できます。

dαf(x)=limh01hαm=0Γ(α+m)Γ(m+1)Γ(α)f(xmh),αR d^{-\alpha} f(x) = \lim_{h \to 0} \frac{1}{h^{-\alpha}} \sum_{m=0}^{\infty} \frac{\Gamma(\alpha + m)}{\Gamma(m + 1)\Gamma(\alpha)} f(x - mh) \,, \qquad \alpha \in \Bbb{R}

分数階微積分の実装

Grünwald–Letnikovの分数階微分と分数階積分を再掲します。

dαf(x)=limh01hαm=0(1)mΓ(α+1)Γ(m+1)Γ(αm+1)f(xmh)dαf(x)=limh01hαm=0Γ(α+m)Γ(m+1)Γ(α)f(xmh) \begin{aligned} d^{\alpha} f(x) &= \lim_{h \to 0} \frac{1}{h^{\alpha}} \sum_{m=0}^{\infty} (-1)^{m} \frac{\Gamma(\alpha + 1)}{\Gamma(m + 1)\Gamma(\alpha - m + 1)} f(x - mh)\\ d^{-\alpha} f(x) &= \lim_{h \to 0} \frac{1}{h^{-\alpha}} \sum_{m=0}^{\infty} \frac{\Gamma(\alpha + m)}{\Gamma(m + 1)\Gamma(\alpha)} f(x - mh) \end{aligned}

実装します。

// ガンマ関数の計算で定義した関数 gamma を使用。 function fractionalDerivative(func, x, alpha, h) { var sum = 0 var sign = 1 for (var m = 0; m < MAX_SUM; ++m) { sum += sign / gamma(m + 1) / gamma(alpha - m + 1) * func(x - m * h) sign *= -1 } return sum * gamma(alpha + 1) / h ** alpha } function fractionalIntegral(func, x, alpha, h) { alpha = -alpha var sum = 0 for (var m = 0; m < MAX_SUM; ++m) { sum += gamma(alpha + m) / gamma(m + 1) * func(x - m * h) } return sum / gamma(alpha) * h ** alpha } // 計算する関数とパラメータ。 // alpha が正の値なら微分、負の値なら積分。 var MAX_SUM = 128 var exampleFunc = (x) => { return x * Math.sin(x) } var alpha = 0.5 var h = 0.01 var maxI = 1024 - 1 var maxX = 10 var result = [] for (var i = 0; i <= maxI; ++i) { var x = maxX * i / maxI var deriv = fractionalDerivative(exampleFunc, x, alpha, h) var integ = fractionalIntegral(exampleFunc, x, alpha, h) result.push([deriv, integ]) } console.log(result)

デモ

黒のプロットが fractionalDerivative 、赤が fractionalIntegral です。

ブラウザの開発者コンソールから関数 func を上書きすることで描画する関数を変更できます。 func の第一引数には横軸の座標が入力されます。コンソールから再描画するには関数 refresh を呼び出してください。

func = (x) => x**5 refresh()

問題点

分数階微分についての理解が足りていないのでテストケースが用意できていません。

0<α<10 < \alpha < 1 の範囲を超えると計算結果がおかしいように見えます。積分方向は見当違いの値がでているようにしか見えません。何とかして整数階の微分と小数階の微分で分けて計算することが考えられます。

参考文献