多重振り子のルンゲクッタを使った計算ができません。

例えば、二重振り子をルンゲクッタを使って計算したいのですが。
dy1/dt=x1 (1)
dy2/dt=x2 (2)
dx1/dt=-A*dx2/dt*cos(y2-y1)-A*x2*sin(y2-y1)-U1 (3)
dx2/dt=-A*dx1/dt*cos(y2-y1)-A*x1*sin(y2-y1)-U2 (4)

この式のまま計算すると発散してしまいます。しかし、(3)式と(4)式を連立してdx1/dt、dx2/dtについて解いた式を使うとうまく計算できます。
この方法だと振り子が増えていった場合計算できません。(1)～(4)式のままルンゲクッタを使って計算したいのですが、いい方法知りませんか？
お願いします！教えてください。

No.5 回答者： noname#108554 回答日時： 2010/01/08 00:27

M y ' = f(x,y)（Mはある行列、yはあるベクトル）のような微分方程式の場合、

陰解法を用いることができます。与式の場合には、
(3)(4)右辺第一項を左辺に移動することで、M y ' = f(x,y)の形にできます。
陰解法アルゴリズムとして、私が知っているのは、RadauIIA法で、
http://www.unige.ch/~hairer/software.html
からFortranコードが入手できます。

とはいってもFortranは私はほとんど忘れてしまって、
実行したことはないのでアドバイスはできません。
C++のコードなら、最近、stiff.tar.gz をダウンロードして実行しましたので、
不明な点があればコメントできるかもしれません。

No.4 回答者： kfnorisu 回答日時： 2009/09/23 15:15

　こんにちわ。

自分も同じような、多重振り子のシュミレーションを作れないかと考えています。計算の使用用途が分からないのですが、多少時間を用いればクラメル自体の計算はできると思いますが、それでは不十分なのでしょうか？

No.3 回答者： eatern27 回答日時： 2008/12/17 03:34

>しかし、(3)式と(4)式を連立してdx1/dt、dx2/dtについて解いた式を使うとうまく計算できます。

ってのは、(3),(4)を手計算などで(dx/dtについて)『解析的に』解いて(つまり,dx/dtをx,yの関数として『具体的に』書き下して)、その式をプログラムに打ち込んでいるんですか？

もしそうなら、わざわざ『解析的に』解く必要性は何かあるのでしょうか？

No.2 回答者： boobee0125 回答日時： 2008/12/16 22:44

私に答えられそうにもない高度な質問に口出ししてごめんなさい。

ちなみに(3)、(4)のままでルンゲクッタ法を使う場合は右辺の微分 dx~/dt = g~(z~) の中間値：
g~(zi~+αj*k0~+βj*k1~+...) （z~は x~、y~ の合成ベクトル）
はどのように算出、更新されているんでしょうか？これってエクセルの循環参照のような形になりませんか？

>連立して右辺に微分が含まれない形にすると１つの運動方程式の中に100の階乗個の項が入ってしまいます．

100の階乗個とは天文学的な数字ですが、これは100の２乗程度ではないんでしょうか？であればコンピューターブン回しの術もありえますが。しかし100個の多重振り子となるとほとんど連続体の世界ですから（可能ならば）偏微分方程式の数値解法の方が良さそうな気もしますね。

この回答への補足

ご意見ありがとうございます。中間値は今手計算してみているのですが、確かにエクセルの循環参照のようになります。このようにならないような方法をさがしています。

また、運動方程式の中に１００の階乗個の項が入るといいましたが、勘違いでした。すみません！ただ、クラメルの公式（クラメルの公式の条件を満たしていると仮定しています）を使って右辺に微分が含まれないような形にしようとしているのですが、その計算過程で１００の階乗個の計算をしなければならないと思うのです。

偏微分方程式の数値解法見てみます！

補足日時：2008/12/17 21:21

No.1 回答者： boobee0125 回答日時： 2008/12/15 20:54

(1)、(2) はともかく、(3)、(4) は右辺に微分が含まれていてルンゲクッタ法を始めとする微分方程式のデジタルシミュレーション（あるいは状態方程式）の形式になっていませんのでムリな話ではないでしょうか？たとえ発散しなくても正しい結果は得られないと思います。

なぜ dx1/dt、dx2/dtについて解いた式が使えないんでしょうか？

この回答への補足

ご意見ありがとうございます．最終的には振り子が100個繋がった多重振り子の計算をしなければなりません．その場合，100個の運動方程式が導かれ，連立して右辺に微分が含まれない形にすると１つの運動方程式の中に100の階乗個の項が入ってしまいます．ですので，(1)～(4)式のままで計算したいのです．おそらく，ルンゲクッタの式を変える必要があると思うのですが，どうしていいかわかりません．

補足日時：2008/12/16 18:03