ODE (Open Dynamics Engine)のメーリングリストを見ていたら，積分器に４次のルンゲ-クッタが実装されそうです。これはビッグニュースです。ODEは今まで高速で安定性があり，リーズナブルな挙動を示すことが特徴でした。逆に言えば，精度が必要な分野では使用しない方が良いということです。これはODEのユーザガイドにも書いてありますし，拙著にも明記しています。
ところが，このビッグニュースにより唯一の弱点が解決される可能性があるのです。ロボット研究者も益々利用できる機会が増えることでしょう。
ここからは専門的になるので，初心者は決して読まないでください（笑い）。
ODEは積分法としては，準陰積分法(semi-implicit integration)を採用しています。準陰積分法は，陰積分（後退積分）と陽積分（前進積分）が混在している手法です。具体的には，拘束力に関しては陰積分が適用され，外力に関しては陽積分が適用されています。現時点では，陽積分として一次のオイラー積分を用いています。
ODEでは拘束力計算にAnitescuの方法を用いているため，一次の陰積分が組み込まれた運動方程式と拘束条件を連立してLCP問題に帰着させています。そのため，拘束力の計算は陰積分で求めています。この拘束ベース法はODE特有の方法というより，極めて広く使われている方法のようです。例えば，MathEngine，Havok, PhysX, OpenHRP3，Springheadなど…
ODEの陰積分（拘束計算）の計算精度に関しては，BoeingとBraunlによるとODEはPhysX, Novodex, Bullet, JigLib, Newton, Tokamak, True Axisと比較した結果，一番良い ようです。ただし，dWorldStep()を使う場合です。精度が必要な場合はdWorldQuickStep()を使わないでください。
なお，この論文では結論として，多くのテストで一番は入れ替わっているので用途に合った計算エンジンを選ぶのが良いと述べています。
私のそのとおりだと思います。ただ，その論文のテストではODEのパラメータが適切に設定されていないように思えるので，全体的にモット良い結果になると思いますが…
最後に，物理計算エンジンの性能は時と伴に変わります。研究者が優れたアルゴリズムを毎年発表しているので，結局はマン・パワーの勝負になるでしょう。2000年前半まで一世を風靡したMathEngineからODEやPhysXが産まれました。最近では，NVIDIAがAGEIA PhysXを買収し，IntelがHavokを買収しています。ODEがこれから生き残るためには，Linuxのようにユーザや利用する企業を増やすことでしょう。demura.netが少しでも，そのお役に立てれば幸いです。
参考文献
A. Boeing and T. Braunl, Evaluation of real-time physics simulation systems, Proceedings of the 5th international conference on Computer graphics and interactive techniques in Australia and Southeast Asia, P.281-288, (2007)
コメント
skasaiさん，
ご存知だと思いますが，ODEで評価するときに，CFMとERPの値に注意してください．完全な剛体としてシミュレーションするためにはERPを1.0，CFMを0.0としなければなりません．
また，精度としてはdWorldStep()がdQuickStep()より精度が高く，dWorldStep()の引数は積分の時間刻み幅です．
でむ
demu さん，アドバイスをありがとうございます．
ScilabとODEの計算結果には差異が無い，と記述しましたが，その考察は間違いのようです．
差をとっただけですと微小な差異(誤差)と言えると思いますが「誤差率 = (Scilab-ODE) / Scilab」で評価すると，とてもとても大きい値です．60秒間の計算では誤差率でプラスマイナス100％を超えてしまう時もありました．
ということで，先のエントリーは間違いと思いますので，ここに訂正させていただきます．
ODEの計算精度の評価方法とその結果についてはとても興味あるところですので，少しずつ試みたいと思っています．
skasaiさん，
コメントありがとうございます。
さて，ODEの精度が悪いとされている点は外力の計算に一次のオイラー法を使っていることだといわれています。それ以外の拘束力に関しては，他の多くの動力学計算エンジンと同じように接触問題をLCP問題に定式化し解いています。
理論的にはODEが他の計算エンジンと比較して精度が悪い点は，オイラー法とルンゲクッタとの差がメインだと考えられます。詳しくは，定量的に動力学計算エンジンの精度の差を調査した上の参考文献を読んでください。Googleで見つかると思います。それによるとPhysXよりは精度が良いとのことでしたが，OpenHRP3との比較はありませんでした。
ODEとOpenHRP3の精度に関する定量的な評価は必要だと常々考えていますが，OpenHRP3の使い方が分からないことと，ODEにルンゲクッタが実装される予定なのでそれが済んでから評価しようと考えていました。なお，ODE0.10.1ではまだルンゲクッタは実装されていません。
skasaiさんのご質問の答えになっていませんが，ODEの評価は，理論値，他の動力学計算エンジン（例えばOpenHRP3)との比較，外力と拘束力（接触力）に分け，更に，dWorldStep(), dWorldQuickStep()などに場合分けするといろいろなことが分かると思います。ルンゲクッタの実装の状況に関しては，ODEのメーリングリストに投稿して聞いてみると良いでしょう。
でむ
お世話になっています．
ODE 0.10.1 を使用しています．
その計算誤差を定量的に把握したいと思って，長さ1mの一様な棒を端点を固定軸においた振子を計算させてみました．固定軸にはHingeを使用しています．
振れ角θの初期値を1rad，角速度θdotを0とし，5秒間動作させθとθdotの結果を得ました．
一方で，ScilabでFehlbergのRunge-Kuttaを用いて運動方程式を解きθとθdotを求め，ODEの結果と比較しました．
その結果，それらには差異が無くてODEのθもθdotもちゃんと計算されているよう感じます．
ODEの計算誤差を評価する目的として，振子の運動の解析は不十分なのでしょうか，それとも，ODEのこのバージョンでは計算精度が向上している，とみるべきでしょうか．
ご意見をいただけたらありがたいです．
Cさん，
長文のコメントありがとうございます．ウェブサイトを開設しているとコメントがないと寂しいものですから．
さて，数値計算の発散を抑えるためにはステップサイズを小さくすれば良いです．その代償として速度が遅くなります．
また，主慣性モーメントの話は，外力計算が顕著に働く例です．ODEは外力計算が一次のオイラー法なので，外力が繰り替えし顕著に働く場合は数値計算が発散します．
最後にフリーの動力学計算エンジンの話ですが，Open HRP3は４月から自由にダウンロードできるようになると聞いています．ただ，速度はODEより遅いと予想されます．ODEでもオイラー積分をα版４次ルンゲクッタに変更すると，数倍遅くなるという報告があります．
精度と速度はトレードオフの関係にあるので，アプリケーションの要求性能に合わせた方法を使うというのが一般的な（あまり役に立たない）回答だと思います．
現実的には，ロボワンのヒューマノイドロボットは剛性が不足していますし，モータのギアなどのガタツキなどあります．ロボットが剛体として扱えないので剛体動力学計算エンジンで厳密に計算しても，現実とは違うのではないでしょうか（つまり，ＯＤＥでプログラムのバグや概略解を探し，後は実機で）．
産総研のPrometなどのように剛性の高いヒューマノイドはOpen HRP3を使わなければいけないと思いますが，ラジコンサーボと1mmぐらいのアルミ板，ラジコンサーボで作った手作りのヒューマノイドは実機で動かすのが一番でしょう．
でむ
お返事ありがとうございます。
そして、この話(新しい計算理論やその中身)自体が私の理解力を越えているのでこれ以上入り込むことができません。^^
なんせ、紹介していただいている初心者向けの解説すら線形代数を初めとする数学的基礎知識の無さから手に負えないで難渋しています。
（demuさんの本の式すら入り口から読みこなせない水準です。）
demuさんが紹介されている範囲内だと思うのですが、
ODEでは、モーメント差の大きな物体の回転で発散しやすいという話ですが、
精度も含めて、時間の刻みを十分に小さくすることで性能向上したり
発散を押さえたり出来ないのでしょうか？
素人考えですがモーメント差の大きな物体（以下「棒」）が回転するときに
無数に存在する質量の点の運動ベクトルは、物体の形状としての拘束が無ければ
砂のように拡散していきますよね？
そこで、砂の一粒が棒という構造物に与える力で得る高校数学程度の
単純物理を、棒を構成する砂粒の数だけ重ね合わせて、微小時間の
棒の姿勢変化を得るって感じなんでしょうか？
あと、現在、ROBO-ONE on PCなどで、ロボットの物理シミュレーションが
ありますが、使用できるツールが参加期間中限定なので、
通常の開発には応用できない不便さがありそうですね。
以前は商品に解析ソフトは無くてもCADがあったようですが、
今回はCADも無かったようです。
今後はアマチュア歩行ロボットも、現在のアニメ歩行＋古典制御理論じゃなくて、
運動学を利用した計算歩行に向かっていくと思うのですが、
末端のユーザーでも経済的障害無しに使えるツールが無いのが問題と
思われます。
ODEが棒状の回転に弱い・・の問題の程度はやってみないと分からない事だとは
思うのですが、demuさんの言われるロボットのフィールドでの粒度での
ODEを、なんとか歩行制御に運用する事は出来そうにないでしょうか？
現実問題、現実の歩行ロボットは、ガタや剛性からの揺れなどもあって、
わけのわからない外乱を押さえる最後のツメは、
振動をフィルターにかけて周期や遅れを含めた上で補正していく
泥臭い処理で押さえ込むことに成ると思われます。
そういう意味でも１００％を求める精度の運動学でなくとも、
どうせ実物は１００％とは離れたところにあるので８０％までを
ODEで予測して、残りの20%を泥臭く押さえ込むというふうなプランの上での
ODEでの可用性が気になります。
たとえば別のアプローチ。ODEで計算しやすいロボットモデルが
あれば、そうしてしまうとか、一歩丸ごとだと駄目なら、
駄目にならない水準まで時間か構造を分解して計算強度を避けてしまうとか・・。
高精度な歩行計算ではHRPなどがあるようですが、今はDL公開していないですし、
なんというか、ヒューマノイド用として縛った汎用性や、
こうして公開しないといった長期的な不安定性はちょっと嫌ですよね。
Cさん、
核心を突いたリプライありがとうございます。
ODEはまだ枯れていなくて，ODEの改良は日々進んでいるという意味です。ルンゲクッタは古典的な方法ですが，ＯＤＥの弱点の一つであった計算精度を向上させることが期待されます。
ODEは準陰積分法を使っています。拘束力計算には陰積分、外力には一次のオイラー法を今のところ使っています。陰積分は数値計算を安定化するので良いのですが、オイラー法は発散させます。しかも、計算誤差が大きいので、それが顕著になります。
ODEのマニュアルに書いてますが、主慣性モーメントの差が大きい剛体（細長いものなど）を回転させると、数値計算が発散し、回転が加速し続け、どこかへ飛んでいく現象が見られます。
ODEで陰積分を使っている理由は，StewartとTrinkleの方法を使っているからです。彼らの方法は，拘束力を求めるために運動方程式と，積分が陰に埋め込まれた速度に関する拘束条件を連立して解いているためです。
これは今では特別な方法というわけではなく，いろいろな計算エンジンで使われている一般的な方法とのことです。
Stewartの方法については以下が良く引用されています。
D.Stewart and J.Trinkle: An Implicit Time-stepping Scheme for Rigid Body Dynamics with inelastic collisions and Coulomb Friction, International Journal of Neumerical Methods in Engineering, 1996.
また，
Kenny Erlebenのドクター論文が参考になると思います。
Kenny Erleben, Stable, Robust, and Versatile Multibody Dynamics Animation, Ph.D. dissertation, University of Copenhagen, 2004.
なお，demura.netは初心者向けにODEの情報を発信していますので，これ以上専門的なことは，ご本家様のメーリングリストでお聞きになられると良いと思います。
でむ
はじめまして。
全然分かってないのですが、４次のルンゲ-クッタが実装されたことと、
後半の準陰積分法、陰積分・陽積分などの様々な計算方式にはそれぞれに
最適な方法があってどんどん進化してますねー。というのはどう繋がっているのでしょうか？
４次のルンゲ-クッタが実装されたことで、最適な計算方式が１つ増えたという事でしょうか？
４次のルンゲ-クッタが実装されたことで、準陰積分法より先進的な性能を獲得したのでようか？
ちょっと検索したら日本語の論文が見つかりました。
http://haselab.hi.mce.uec.ac.jp/files/pubs/01078_pdf_Tazaki_Hasegawa_gcad2006.pdf
この分野は全くの素人なんですが、シミュレーションでの安定性とは
離散時間と誤差の共振みたいなもののようですね。
昔、BASICでブロック崩しを作ったときにボールが狭いところで反射したときに、
離散時間後には、壁やブロックに埋まった状態から次の時間が発生して
その次もおかしな所で・・って感じでイメージ通りに動かないバグを思い出しました。
そんな感じなのでしょうか。
後退積分というは、この埋まる要素を無くすものなのでしょうか？
理論を追い込む力がないので良く分からないのですが、
拘束条件に対して準陰積分法（条件で陰積分と陽積分を使い分ける）を使うのは、
安定化のためには塗り絵で枠からはみ出ずに、塗りシロが残ってもいいから塗る。
みたいなものなのかなと創造しながら読みました。
アホな投稿ですみません。笑ってやってください。m(_ _)m