λ Haskell λ で導関数を実装

今回はHaskellを用い、速度関数及び加速関数を題に導関数の実装を紹介します。

速度関数

速度関数の$v$は位置関数の$x$の導関数ですが、数学の授業に出てくる導関数と多少違います：

$$v(t) = \lim\limits_{\Delta t -> 0} \dfrac {x(t+\frac {\Delta t}2)-x(t - \frac {\Delta t}2)}{\Delta t}$$

数学なら：$f(x) = \lim\limits_{h -> 0} \dfrac {f(x+h)-f(x)}h$ですが、$v$の入力は$x+h$と$x$ではなく、$t + \frac {\Delta t}2$と$t - \frac {\Delta t}2$になります。
つまり、$f$は間隔の冒頭を重視するのに対し、$v$は$\Delta t$間隔内の速度を計測します。

$3m$から$8m$まで移動する時、速度を計算するために間隔の両端の$(t_3,3m)$と$(t_8,8m)$で平均速度を求めることができます。上記の$v$の表記はその考え方を反映しているのではないかと考えます。

velocity dt x t = (x(t + dt / 2) - x(t - dt / 2)) / dt 

ご覧の通り、極限を真似る一番簡単な方法ですからdtを指定します。上記を数学に再訳すと$t$に加え、$dt$までもらう速度関数になります：

$$v(t, dt) = \dfrac {x(t+\frac {dt} 2) - x(t - \frac {dt} 2)} {dt}$$

ただし、数学的微分と違って、この時点でdtを小さくするのではなく、大きくするときの方は精密度が高いです（手頃なためGHCiという対話環境を使います）：

ghci> velocity dt x t = (x(t + dt / 2) - x(t - dt / 2)) / dt
-- 15秒たてば速度はなーんだ？
-- 当たり前ですけど、x(t) = 30*tだと速度は常に３０
-- でも、あれ？結果は少し違います
ghci> velocity 0.001 (\t -> 30*t) 15
30.000000000086402
ghci> velocity 0.01 (\t -> 30*t) 15
30.00000000000682
-- dtを大きくすると正解が出る
ghci> velocity 0.1 (\t -> 30*t) 15
30.0
ghci> velocity 1 (\t -> 30*t) 15
30.0

期待に反してdtを大きくすればするほど正解に近づきます。なぜこうなるのかといいますと、使っている関数（\t -> 30*t）が単純すぎるからです。速度は常に等しいため、計測する間隔をどんなに伸ばしたり縮めたりしても平均速度は変わりません。
これに対して、間隔を短くすると浮動小数点数の限界に達し、精密度が下がってしまいます。

では、速度ができるだけ頻繁に変わる例を考えましょう。例えば、くねくねとした図表を持っている正弦関数なら良いでしょう。
正弦関数の導関数は余弦関数ですので、$sin'(0)=cos(0)=1$。

ghci> velocity 1 sin 0
0.958851077208406
ghci> velocity 0.1 sin 0
0.9995833854135666
ghci> velocity 0.01 sin 0
0.9999958333385416
ghci> velocity 0.001 sin 0
0.9999999583333338
ghci> velocity 0.0001 sin 0
0.9999999995833334
ghci> velocity 0.00001 sin 0
0.9999999999958332

今回はdtを小さくすると精密度がちゃんと上がってくれます。

加速と2階導関数

もちろん、微分ができると、速度だけでなく加速も計算できます。accelerationの実装は上記のvelocityと全く同じですが、名前だけ変えて多少わかりやすくしましょう：

-- 導関数
ghci> acceleration dt v t = (v(t + dt / 2) - v(t - dt / 2)) / dt
-- 速度関数
ghci> v t = 2 * sin t + 30
ghci> acceleration 0.000001 v 0
2.0000000020559128 -- 悪くはない

OK。次は位置関数を渡し、加速を返す2階導関数を作りたいです。そのために、まず積文法を使って2階導関数に渡すべき位置関数を準備します：

$$\int 2 sin(t) + 30 dt = - 2 cos (t) + 30 t + C \;\; (正しい、C = 0 とします)$$

Haskellでは：p t = -2 * cos t + 30 * t

後は、導関数を自分の入力にして2階導関数を作り出します。位置関数の2階導関数ですので一応ddxと名付けましょう：

-- 元の導関数
ghci> dx dt x t = (x(t + dt / 2) - x(t - dt / 2)) / dt 
-- 2階導関数
ghci> ddx dt x t = dx dt (dx dt x) t

最後に、上記のpを渡せば完成：

ghci> ddx 0.0001 p 0
1.9999999656405976

正解です。ただ、今回のdtは割と大きいにしないと浮動小数点数の精密度の限界を超えてしまいますので、ご注意ください。

おまけ：2階導関数の簡略化

念のため、2階導関数の式を簡略化し、dtの精密度を上げることはできないか検討してみます。

簡略化はこちら：

$$\begin{align} v(t, dt) &= \dfrac {x(t+\dfrac {dt}2)-x(t - \dfrac {dt}2)}{dt} \\\\ a(t, dt) &= \dfrac {v(t + \dfrac {dt} 2, dt) - v(t - \dfrac {dt} 2, dt)} {dt} \\\\ &= \dfrac {\dfrac {x(t + \dfrac {dt} 2 + \dfrac {dt} 2)- x (t + \dfrac {dt} 2 - \dfrac {dt} 2)} {dt} - \dfrac {x(t - \dfrac {dt} 2 + \dfrac {dt} 2)-x(t - \dfrac {dt} 2 - \dfrac {dt} 2)} {dt}} {dt} \\\\ &= \dfrac {\dfrac {x(t + dt)-x(t) - (x(t)-x(t - dt))} {dt}} {dt} \\\\ &= \dfrac {x(t + dt) + x(t - dt) - 2(x(t))} {dt^2} \\\\ \end{align}$$

Haskellに戻ると：

ghci> ddx2 dt x t = (x(t+dt)+x(t-dt)-2*x(t))/dt**2
ghci> ddx2 0.0001 p 0
1.999999987845058

dtをこれ以上小さくすると精密度が下がりますので、前と比べると結果は微妙に正解に近づいたぐらいですね。