pythonでローレンツアトラクタ（オイラー法）

wikipediaのローレンツ方程式に書いてあるそのままのパラメータで解いてみた。画像はgnuplotで描画したもの。
同じものをいくつか他の解法で解いて比較すると面白そう。

# coding: utf-8
# lorenz attractor
# scheme: euler

fx = lambda x,y,z,r,p,b: -p*x + p*y
fy = lambda x,y,z,r,p,b: -x*z + r*x - y
fz = lambda x,y,z,r,p,b: +x*y - b*z

def main():
    # steps
    nStep = 100000
    dt = 1e-3 # time step
    x = 1
    y = 1
    z = 1
    t = 0
    # parameters
    p = 10
    r = 28
    b = 8./ 3

    for i in xrange(nStep):
        xx = dt * fx(x, y, z, r, p, b)
        yy = dt * fy(x, y, z, r, p, b)
        zz = dt * fz(x, y, z, r, p, b)
        x += xx
        y += yy
        z += zz
        t += dt
        print x,y,z


if __name__ == '__main__':
    main()

追記：縮めてみた

p,r,b = 10,28,8./3
t,dt = 0,1e-3
q = [1.,1.,1.]
while t < 100:
    x,y,z=q
    print x,y,z
    q = [q[i] + dt*m for i,m in enumerate([-p*x + p*y, -x*z + r*x - y, +x*y - b*z])]
    t += dt