科学と家事とプログラミング (python を中心に)

python 温度計測 湿度計測 DS18B20 USB9097

超構造化(09) プログラムは手順の集まりか?

手続き vs 宣言

「プログラムは手順の集まりである」という表現は、賞味期限つきの標語の良い例だ。 初めてプログラムを書くときに、この助言が役立つ場面は確かにあると思う。 しかしその一方で、さまざまなプログラミング言語が目指すものは 「人間から見た分かりやすさ/扱いやすさ」であり、そこに通底する主題は、 「いかに計算機のように考えないで済ませるか」ということに他ならない。 通常の構造化プログラミングの範囲は超えるが、 python を用いた関数型プログラミングっぽい手法を紹介する。

無限ストリーム

例として、減衰振動のシミュレーションを取り上げる。 力学系を表す微分方程式を関数としてを動的に生成して、 状態量を無限ストリームとして与える。 こういう技法を適用すると、いたって宣言的なプログラムが出来上がる。 python の場合、無限ストリームはジェネレータと呼ばれる。 文法的には、return の代わりに yield という構文を使って値を返す。 例題中の関数 states() は yield t や x を返すが、これらの t, x は、 参照される時に(必要になって初めて)生成される。 一方、使う側では for t, x in stream: という形で、stream が既に出来上がっている (完成している存在である)かのように記述できる。 つまり、未来を先取りできるというわけだ。

コードは、scheme の文法書(r4rs)の末尾にある例題を真似たものである。

"""減衰振動(LCR 回路)の時間シミュレーション"""

class SVEC:
    """ 状態ベクトル (rk4 が要求する演算のみ実装) """
    def __init__(self, x): self.x = x[:]
    def __getitem__(self, j): return self.x[j]
    def __add__(self, o): return SVEC([a + b for a, b in zip(self.x, o.x)])
    def __rmul__(self, f): return SVEC([f*a for a in self.x])
    def __div__(self, f): return SVEC([a/f for a in self.x])

def rk4(t, x, f, h):
    """ 積分ステップ h 後の状態を返す。f は力学系(微分方程式)。"""
    f1 = h*f(t, x )
    f2 = h*f(t + h/2, x + f1/2)
    f3 = h*f(t + h/2, x + f2/2)
    f4 = h*f(t + h , x + f3 )
    return x + (f1 + 2*f2 + 2*f3 + f4)/6

def damped_oscillator(L, C, R):
    """ 減衰振動する系(微分方程式を表す関数)を生成して返す。"""
    def f(t, s):
        v, i = s # s = [C_電圧, L_電流]
        return SVEC([-(i + v/R)/C, # d(s[])/dt = dv/dt
        v/L ]) # di/dt
    return f # 上で生成した関数を返す

def propagation(intg_method, f, h):
    """ intg_method を使った刻み h の伝播 => next_state() """
    def next_state(t, x):
        x = intg_method(t, x, f, h)
        t = t + h
        return t, x
    return next_state
    
def states(derivative, t, x, h):
    """ states は、状態量の無限ストリーム(微分方程式 derivative, 初期値t, x, 刻み h) """
    next_state = propagation(rk4, derivative, h)
    while 1:
        t, x = next_state(t, x) # 次の状態は、直前の状態で決まる。
        yield t, x

if __name__ == '__main__':
    """ 状態量の無限ストリーム """
    s = states(damped_oscillator(L=1e3, C=1e-3, R=10e3),
                0.0, # 初期時刻
                SVEC([1.0, 0.0]), # 初期状態
                0.1) # 刻み幅

    # ストリームを印刷
    for t, x in s:
        if t >= 100: break
        print t, x[0], x[1]

f:id:sken20k:20180212175350j:plain