みんな知ってるルンゲクッタ法. まだ出来ない人のためのPython実装

ルンゲクッタ

最近ドリフトのシミュレータを書いている。シミュレーションにあたって車両の運動を、適当な車両運動モデルを作って時間発展を計算させている。さて、最近つくった私の貧弱モデルをオイラー法で動かすと、どうも数値振動しているのか車体が「グラグラ」して気持ち悪い。運動方程式の立て方も悪いので一概に原因が分からない。方程式の立て方で復元力によって振動しているかもしれないし、数値微分の振動かもしれない。と。原因はいろいろ考えられるわけだ。原因は一つ一つつぶしていくのがセオリーなので、ひとまずオイラー法を精度の良いルンゲクッタ法に置き換える。(オイラー法では車両が上手く動かせないことは以前に確認済み)。今回は車両の回転方程式と並進の2方程式を解く。変数は重心の横滑り角度と、角速度の二つなので、れんせい方程式のルンゲクッタを計算する練習で、ローレンツ・アトラクターをお題としてベースのコードを書く。

まぁ、こちらの記事をpythonで書き直しただけでございます。 ルンゲクッタ法で様々な連立微分方程式を解く数値計算例(C言語)

まずはそのまま書き下したコード

# -*- coding:utf-8 -*-
"""
runge kutta

http://www.geocities.jp/supermisosan/rksimultaneousequation.html


:Equation:
dxdt = 10*(y-x)             --- (1)
dydt = 28*x - y - x*z       --- (2)
dzdt = -(8./3.)*z + x*y     --- (3)

"""
import numpy as np

def f1(t,x,y,z):
    return 10.*(y - x)

def f2(t,x,y,z):
    return  (28.*x) - y - (x*z) 

def f3(t,x,y,z):
    return  (-8./3.)*z + (x*y) 


def main():
    # time step
    dt = 0.01
    tmax = dt * 10000
    # initial condition
    t = 0.0
    x = 0.
    y = 1.
    z = 1.05

    k0=[0,0,0]
    k1=[0,0,0]
    k2=[0,0,0]
    k3=[0,0,0]

    write(["t","x","y","z"], header=True)

    while t<=tmax:

        k0[0]= dt * f1(t,x,y,z);
        k0[1]= dt * f2(t,x,y,z);
        k0[2]= dt * f3(t,x,y,z);

        k1[0]= dt * f1(t+dt/2.0, x+k0[0]/2.0, y+k0[1]/2.0, z+k0[2]/2.0);
        k1[1]= dt * f2(t+dt/2.0, x+k0[0]/2.0, y+k0[1]/2.0, z+k0[2]/2.0);
        k1[2]= dt * f3(t+dt/2.0, x+k0[0]/2.0, y+k0[1]/2.0, z+k0[2]/2.0);

        k2[0]= dt * f1(t+dt/2.0, x+k1[0]/2.0, y+k1[1]/2.0, z+k1[2]/2.0);
        k2[1]= dt * f2(t+dt/2.0, x+k1[0]/2.0, y+k1[1]/2.0, z+k1[2]/2.0);
        k2[2]= dt * f3(t+dt/2.0, x+k1[0]/2.0, y+k1[1]/2.0, z+k1[2]/2.0);

        k3[0]= dt * f1(t+dt, x+k2[0], y+k2[1], z+k2[2]);
        k3[1]= dt * f2(t+dt, x+k2[0], y+k2[1], z+k2[2]);
        k3[2]= dt * f3(t+dt, x+k2[0], y+k2[1], z+k2[2]);

        dx = (k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0;
        dy = (k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0;
        dz = (k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0;

        x = x + dx
        y = y + dy
        z = z + dz

        write([t,x,y,z])
        print(t,x,y,z, k0,k1,k2,k3)
        t = t + dt

def post(filename="lorenz.csv"):
    import numpy as np
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D

    data = np.genfromtxt("lorenz.csv", delimiter="," ,filling_values=(0, 0, 0, 0))
    ts = data[:,0]
    xs = data[:,1]
    ys = data[:,2]
    zs = data[:,3]

    fig = plt.figure()
    ax = fig.gca(projection='3d')

    ax.plot(xs, ys, zs)
    ax.set_xlabel("X Axis")
    ax.set_ylabel("Y Axis")
    ax.set_zlabel("Z Axis")
    ax.set_title("Lorenz Attractor")

    plt.show()

def write(line, header=False):

    import csv

    # Header
    if header:
        with open('lorenz.csv', 'w') as f:
            writer = csv.writer(f, lineterminator='\n') # 改行コード（\n）を指定しておく
            writer.writerow(line)     # list（1次元配列）の場合
    # Body
    if header == False:
        with open('lorenz.csv', 'a') as f:
            writer = csv.writer(f, lineterminator='\n') # 改行コード（\n）を指定しておく
            writer.writerow(line)     # list（1次元配列）の場合
if __name__ == '__main__':
    main()
    post()

そして、さらに省略できたコード

mplot3d example code: lorenz_attractor.py — Matplotlib 1.5.3 documentation

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def lorenz(x, y, z, s=10, r=28, b=2.667):
    x_dot = s*(y - x)
    y_dot = r*x - y - x*z
    z_dot = x*y - b*z
    return np.array([x_dot, y_dot, z_dot])


t=0
dt = 0.01
stepCnt = 10000

# Need one more for the initial values
xs = np.empty((stepCnt + 1,))
ys = np.empty((stepCnt + 1,))
zs = np.empty((stepCnt + 1,))

# Setting initial values
xs[0], ys[0], zs[0] = (0., 1., 1.05)
k0=[0,0,0]
k1=[0,0,0]
k2=[0,0,0]
k3=[0,0,0]

# Stepping through "time".
for i in range(stepCnt):
    x,y,z=xs[i],ys[i],zs[i]

    k0 = dt * lorenz(x,y,z)
    k1 = dt * lorenz(x+k0[0]/2., y+k0[1]/2., z+k0[2]/2.)
    k2 = dt * lorenz(x+k1[0]/2., y+k1[1]/2., z+k1[2]/2.)
    k3 = dt * lorenz(x+k2[0], y+k2[1], z+k2[2])

    dx = (k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0
    dy = (k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0
    dz = (k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0

    xs[i+1] = xs[i] + dx
    ys[i+1] = ys[i] + dy
    zs[i+1] = zs[i] + dz


fig = plt.figure()
ax = fig.gca(projection='3d')

ax.plot(xs, ys, zs)
ax.set_xlabel("X Axis")
ax.set_ylabel("Y Axis")
ax.set_zlabel("Z Axis")
ax.set_title("Lorenz Attractor")

plt.show()

実装が出来たらドリフトの動画も見せたいですね。

理工学のための数値計算法 (新・数理工学ライブラリ 数学)

わかりやすい数値計算入門

やさしく学べるC言語入門―基礎から数値計算入門まで (UNIX&Information Science)

だれでもわかる数値解析入門―理論とCプログラム

妄想―彼女はなぜ狙われたのか

僕と彼女のゲーム戦争10 (電撃文庫)

ラジオでディスコを録音三昧。ラジ録１０

こんにちわふぃふぃです。

私は音響系の仕事をしていることもあり、四六時中クラブ音楽やクラシックばかり聴いています。

そんな昭和の最後に生まれた私ですが、実はDISCOソングが大好きです。

もし生まれ変われるならジュリアナやマハラジャに言ってみたかった。。

ある土曜日カーシェアで車を借りて大きな買い物をしながらドライブをしてました。CDを持っていくの忘れたのでふとラジオをつけるとなんとDISCOが流れてきた！！聴いているとどうも横浜のLuther(ルーサー)でもまわしているDJ OSSHYさん・・・。なんて日だ。昔川崎にいた時によくLutherにいっていたの思い出す。こんな所で出会うとは。

DJ OSSHY OFFICIAL SITE

Luther | 横浜駅西口のBar

家に帰って調べてみると毎週土曜日14:00-17:00に放送されているとのこと。なんてタイミングだ。音楽は探して出会うのも楽しいけど。本当にたまたま出会うことが1年に1回ぐらいある。これがたまらないー。

どうしても家で聴きたいんだけどFMラジオない。 録音して永久保存したいけどFMラジオない。。 土曜日家にいない。。。

そこでネットを調べているとどうもPCでラジオ番組を録音できるみたいだ。

ラジ録10 Windows版 [ダウンロード]

予約録音機能のついたラジオもある。

東芝 Bluetooth対応SD/USB/CDラジオTY-CWX85(W)

予約機能のついたラジオはあまり格好良くないし、

有料2000円だがラジオを買うことを考えれば安い。

体験版があるので試用レビューします。

ラジ録10

体験版をダウンロード ラジ録10 Windows版 [ダウンロード]

メイン画面

安定のWindowsグレーのUIです。最近は"おしゃれ"なUIのアプリが主流ですが、定期録音みたいな機能を信頼して任せなら、Windowsグレーだとなぜか安心します。「デザインよりも安定動作に注力してます！」みたいな。

番組一覧

番組一覧が見れます。モタモタもせずさくっとでてきます。安定のWindowsグレーです。リストでだされるよりも便利。

さすがラジ録という名前負けしてません。右クリックするとさっそく予約録音の指定ができます。これは便利。やりたいことはもう終わりですね。

非常にシンプル

アプリを起動して、なにも迷わずに3分足らずで定期録音の予約ができた。

シンプルさは最強だ。

私のなかでは「シンプルイズザベスト」はもう古い。ジョブスの言葉を借りれば「Think Simple」

シンプルに思考し、行動せよ。『Think Simple』の著者が語る、アップルだけ ...

Think Simple―アップルを生みだす熱狂的哲学

こういうときに体験、実感して腹におちますね。

音質

参考までに音質についてもレビューしておきます。

RN3.mp3

ノイズはすくないです。ラジオなので解像度が低いですけどラジオらしい音です。音声用にボーカル帯域がかなりシャンシャンしてますがイコライジングすれば問題ないでしょう。あくまでもラジオ。録音してまで聴きたい人には十分じゃないでしょうか。ラジオを感を楽しみましょう。

radiko.jp

ONKYO ネットワークCDレシーバー ハイレゾ音源対応 ブラック CR-N755(B)

ONKYO AVレシーバー 7.1ch Dolby Atmos/DTS:X/4K/HDCP2.2/ハイレゾ音源対応 ブラック TX-NR656(B)

TEAC ネットワーク/CDプレーヤー ハイレゾ音源対応 CD-P800NT-S

ゾックス パソコンでワンセグテレビとFMラジオを楽しめるUSB接続ワンセグチューナー ブラック DS-DT310BK

プリンストンテクノロジー USBラジオチューナー デジ造ラジオ版II PCA-RCU2

3ぷん。

種々の数学ソフトウェアまとめ | Maple, Maxima, GeoGebra, Knoppix/Math

こんにちわ。ふぃふぃです。

最近は物体の衝突問題に挑戦して頭を悩ませています。久しぶりに古典力学の理論を拾いながら研究をしています。論文をまとめていると数式が多すぎて疲れました。そこで、数式を扱うためのソフトを調べてみたので紹介します。

数学ソフト

Maple(メープル. マップルではありません) 有料

http://www.biwako.shiga-u.ac.jp/sensei/kumazawa/tex/maple.html

Maxima: 無料

Maxima（マキシマ）は、LISP で記述された数式処理システムである。GNU GPL に基づくフリーソフトウェアであり、現在も活発に開発が続けられている。Maple や Mathematica などの商用の数式処理システムと比べても遜色のない機能を持っている。https://ja.wikipedia.org/wiki/Maxima

GeoGebra: 無料

GeoGebra は、幾何、代数、解析を結び付けた動的な数学ソフトウェア。ソースコードはGPL の下でリリースされている。リンツ大学 (Johannes Kepler University of Linz) の Markus Hohenwarter 教授を中心とするグループにより開発が進められている。https://ja.wikipedia.org/wiki/GeoGebra

GeoGebra日本

サンプルコードがたくさんあり参考になる。

Sage

SageMath (セイジ、以前はSage、SAGEと記した) は数学の幅広い処理を扱うソフトウェアである。扱う処理は計算機代数、組み合わせ、数値計算など多岐に及ぶ。工学的応用に加え基礎科学の研究もカバーする。 SageMath は2005年2月24日にフリーソフトウェアとして GNU General Public License の元で初版が公開された。その開発目的は Magma、Maple、Mathematica (いずれも計算機代数ソフトウェア)、MATLAB の代替となるフリーかつオープンソースなソフトウェアを提供することであった^[1]。開発は、米ワシントン大学の数学准教授のウィリアム・スタイン (William Stein) が主導して始まった。 SageMath は Python プログラミング言語を使用しており，手続き型・関数型・オブジェクト指向によるプログラムの記述を行うことができる。  

MSI Nightblade MI2-001BUS Intel Skylake B150 Chipset LGA 1151 i7/i5/i3/Pentium DDR4 Gaming Desktop Barebone PC [並行輸入品]

【Amazon.co.jp限定 】 ASUS ゲーミングモニター 23型フルHDディスプレイ (応答速度1ms / HDMI×2,D-sub×1 / スピーカー内蔵 / 3年保証) VX238H-P

KOTION EACH ゲーミングヘッドホンG2200ステレオPC/Laptopゲーム用ヘッドセット 7.1chバーチャルサラウンド USB接続 重低音 振動機能 密閉型 騒音抑制機能 LEDライト 専用ソフトウェア CD付属の有線ヘッドバンド ソフトレザー式 ノートパソコン デスクトップPC Windows XP/7/8/8.1/10 など対応Gaming Headset (ブラック&ブルー)

俺の彼女と幼なじみが修羅場すぎる１２ (GA文庫)

冴えない彼女の育てかた 11 (ファンタジア文庫)

Sレア装備の似合う彼女（２） (裏少年サンデーコミックス)