grunwald_letnikov

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

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

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

\[ D^{\frac{1}{3}} \left( D^{\frac{1}{3}} \left( D^{\frac{1}{3}} f(x) \right) \right) = \frac{d}{dx} f(x) \]

つまり \(1 / q\) 階微分を \(q\) 回繰り返すと1階微分になります。

\[ \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) \]

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

\[ \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 \]

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

\[ 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 > n\) のとき \(\binom{n}{m} = 0\) なので総和の上限を無限にできます。

\[ 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} \]

\[ 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} \]

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

\[ \binom{n}{m} = \frac{n!}{m!(n - m)!} \,,\qquad n \in \Bbb{N}\,,\quad m \in \Bbb{N} \]

\[ \binom{\alpha}{m} = \frac{\Gamma(\alpha + 1)}{\Gamma(m + 1)\Gamma(\alpha - m + 1)} \,,\qquad \alpha \in \Bbb{R}\,,\quad m \in \Bbb{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\) が入力されても計算できる形になりました。

ガンマ関数に入力される値が \(0\) あるいは負の整数のときは無限大になるので、計算できない場合があります。Grünwald–Letnikovの分数階微分の式では \(\alpha > 0\) かつ \(m > 0\) かつ \(\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
}

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

\[ \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} \]

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

\[ \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} \]

\[ 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} \]

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