忍者ブログ

明日の設計図

たまにロボットを考えるブログ・・・。

[PR]

×

[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。

コメント

現在、新しいコメントを受け付けない設定になっています。

ルンゲクッタ法で微分方程式(単振動)を積分

創生×⇒創成○

だそうです。
大きな問題ですね。

シミュレーションに迫られておりまして、非線形微分方程式を積分したくなっております。

解析解は出ませんどころか、点群データなわけで。

シミュレーションに迫られて います。

5年前にこのブログでもルンゲクッタ法にふれていましたが。
すっかり忘れたのでもう一度記事にします。

参考サイト様

簡単なシミュならオイラー法を使います。
dx/dt=Ax
であれば,刻みのステップをhとして
x1=x0+A*x0*h
という漸化式なので次々と計算すればいいです。
ただ、1次の傾きを積分していくので誤差が積もります。
そこで誤差の少ない方法が世間には広まっているのです。

その例としてルンゲクッタを選びました。
↑の式で進めると
a b c dという4つの係数から、理想的な傾きを求めるようです。
数値計算の教科書とかだと、一般的な常微分方程式を念頭に置いているので、単純な式だと返って分かりづらいものでした。
引数?がだいぶ減るので簡単です。
例えば、
dx/dt=A*t*x
みたいな場合であれば、もう少し教科書寄りの式に成ります。

で、どのくらい精度がいいかということですが


ただのexpカーブが解析解となりますが、このようにオイラーだと離れていくところが付いていくわけです。

でもって、簡単なところはこの辺で。

もう少し応用の利く例を置いておきます。
運動方程式であれば、2階の微分方程式だと思います。
d^2x/dt^2=Bx
みたいな。
Bが負であれば単振動です。

計算は、状態方程式のように2つに分けて、y1,y2の連立の1階微分方程式にします。
 

さっきと同様にa b c d を計算しますが、
y1とy2の流れを交互にやることに注意します。
y2のb、つまりb2はa1を使って計算する必要があります。
少し複雑です。

計算結果です。
B=-1だと思いますが。
ざっくりこんな状態です。

オイラーでもほとんど一緒です。
ですが、解析解との差分をみると
このようになっているので、ルンゲクッタの頑張りようがわかるのではないでしょうか。

ダウンロード

元エクセルを置いておきます。内容は整理していませんのですいません。

さて、2階だろうが、連立だろうができるようになりましたので、
多質点の運動方程式に手が届きます。

それはまた今度。

拍手[0回]

PR

コメント

お名前
タイトル
文字色
メールアドレス
URL
コメント
パスワード Vodafone絵文字 i-mode絵文字 Ezweb絵文字

プロフィール

HN:
Adel
年齢:
34
性別:
男性
誕生日:
1989/09/17
職業:
会社員
趣味:
モチベーション探し
自己紹介:
ロボットつくるのが夢

ブログ内検索

アクセス解析