例えば、二重振り子をルンゲクッタを使って計算したいのですが。
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)式のままルンゲクッタを使って計算したいのですが、いい方法知りませんか?
お願いします!教えてください。
A 回答 (5件)
- 最新から表示
- 回答順に表示
No.5
- 回答日時:
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
- 回答日時:
こんにちわ。
自分も同じような、多重振り子のシュミレーションを作れないかと考えています。計算の使用用途が分からないのですが、多少時間を用いればクラメル自体の計算はできると思いますが、それでは不十分なのでしょうか?No.3
- 回答日時:
>しかし、(3)式と(4)式を連立してdx1/dt、dx2/dtについて解いた式を使うとうまく計算できます。
ってのは、(3),(4)を手計算などで(dx/dtについて)『解析的に』解いて(つまり,dx/dtをx,yの関数として『具体的に』書き下して)、その式をプログラムに打ち込んでいるんですか?
もしそうなら、わざわざ『解析的に』解く必要性は何かあるのでしょうか?
No.2
- 回答日時:
私に答えられそうにもない高度な質問に口出ししてごめんなさい。
ちなみに(3)、(4)のままでルンゲクッタ法を使う場合は右辺の微分 dx~/dt = g~(z~) の中間値:
g~(zi~+αj*k0~+βj*k1~+...) (z~は x~、y~ の合成ベクトル)
はどのように算出、更新されているんでしょうか?これってエクセルの循環参照のような形になりませんか?
>連立して右辺に微分が含まれない形にすると1つの運動方程式の中に100の階乗個の項が入ってしまいます.
100の階乗個とは天文学的な数字ですが、これは100の2乗程度ではないんでしょうか?であればコンピューターブン回しの術もありえますが。しかし100個の多重振り子となるとほとんど連続体の世界ですから(可能ならば)偏微分方程式の数値解法の方が良さそうな気もしますね。
この回答への補足
ご意見ありがとうございます。中間値は今手計算してみているのですが、確かにエクセルの循環参照のようになります。このようにならないような方法をさがしています。
また、運動方程式の中に100の階乗個の項が入るといいましたが、勘違いでした。すみません!ただ、クラメルの公式(クラメルの公式の条件を満たしていると仮定しています)を使って右辺に微分が含まれないような形にしようとしているのですが、その計算過程で100の階乗個の計算をしなければならないと思うのです。
偏微分方程式の数値解法見てみます!
No.1
- 回答日時:
(1)、(2) はともかく、(3)、(4) は右辺に微分が含まれていてルンゲクッタ法を始めとする微分方程式のデジタルシミュレーション(あるいは状態方程式)の形式になっていませんのでムリな話ではないでしょうか?たとえ発散しなくても正しい結果は得られないと思います。
なぜ dx1/dt、dx2/dtについて解いた式が使えないんでしょうか?この回答への補足
ご意見ありがとうございます.最終的には振り子が100個繋がった多重振り子の計算をしなければなりません.その場合,100個の運動方程式が導かれ,連立して右辺に微分が含まれない形にすると1つの運動方程式の中に100の階乗個の項が入ってしまいます.ですので,(1)~(4)式のままで計算したいのです.おそらく,ルンゲクッタの式を変える必要があると思うのですが,どうしていいかわかりません.
補足日時:2008/12/16 18:03お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
似たような質問が見つかりました
- 数学 x=r・cosθの2回微分 θ=ωtとすると? 5 2022/05/10 23:53
- 数学 単振り子とルンゲ・タック法 1 2022/07/15 00:05
- 数学 dl と 、 Δl はどのように違いますか? 2 2022/11/30 15:48
- 数学 dx/dt = |y| , dy/dt = x (-∞<t<∞) をとけ 1 2022/09/17 09:56
- その他(プログラミング・Web制作) Pythonにおける物理のシミュレーションでの単位変換について 2 2023/06/02 17:11
- その他(プログラミング・Web制作) Pythonでのかんたんな物理シミュレーションについての書籍 5 2023/06/02 07:37
- 数学 連立微分方程式 dx/dt = |y| , dy/dt = x (-∞<t<∞) について、 (1) 3 2022/09/16 21:59
- 数学 dx/dt=x-2y +e^t dy/dt=-3x +2y+1 初期値[1,0] [x,y] この連 3 2023/05/15 18:23
- 数学 微分の解けない問題があるので誰か教えて頂きたいです。 dx/dt+2x=cos(4t) x(0)=- 3 2023/06/20 21:15
- 物理学 物体に一定の大きさfの力をx軸の正の向きに加える。またこの物体には抵抗係数がγの速度に比例する抵抗力 2 2023/07/06 04:01
おすすめ情報
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
加速度 a=dv/dt = (d^2 x) /dt^2
-
高校・大学の物理(力学)教え...
-
v^2-v0^2=2ax 今日この式を習っ...
-
力学について質問です。 1.棒の...
-
シュヴァルツシルト時空での時...
-
dH/dtとH(t)の関係
-
EXCEL上の数字を自動で振り分け...
-
なぜ力積を時間に関して積分す...
-
質量流量の記号「・ の読み方を...
-
仕事率の表し方について。
-
Id²θ/dt²=-mghsinθの厳密解の...
-
今基礎物理学の問題を解いてい...
-
運動方程式の微分積分の計算
-
物理で微積をつかう。
-
高校物理 授業でこういうのをや...
-
摩擦を考慮したサイクロイド曲線
-
d^2r/dt^2の意味
-
プランクの放射式
-
ルンゲクッタ法で数値解析(C...
-
運動エネルギー
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
v^2-v0^2=2ax 今日この式を習っ...
-
d^2r/dt^2の意味
-
電流の時間微分、電圧の時間微分
-
質量流量の記号「・ の読み方を...
-
EXCEL上の数字を自動で振り分け...
-
dx/dt=√(1-x^2)の一般解の求め...
-
Debug.Printで表示される内容を...
-
力学について質問です。 1.棒の...
-
物理で微積をつかう。
-
微分積分のdの意味
-
最後のdv/dtは何でしょうか。
-
加速度 a=dv/dt = (d^2 x) /dt^2
-
運動方程式の微分積分の計算
-
微分記号“d”について
-
Id²θ/dt²=-mghsinθの厳密解の...
-
雨滴の運動質量が変化する落体...
-
運動方程式を求めてください
-
機械力学の問題です!!!
-
d/dx=dt/dx * d/dt =d/dt * dt/...
-
地動加速度が単位インパルスの...
おすすめ情報