- 締切済み
多重振り子のルンゲクッタを使った計算ができません。
例えば、二重振り子をルンゲクッタを使って計算したいのですが。 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)式のままルンゲクッタを使って計算したいのですが、いい方法知りませんか？ お願いします！教えてください。
- lks162
- お礼率0% (0/5)
- 物理学
- 回答数5
- ありがとう数1
- みんなの回答 （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 をダウンロードして実行しましたので、 不明な点があればコメントできるかもしれません。
- kfnorisu
- ベストアンサー率25% (2/8)
こんにちわ。自分も同じような、多重振り子のシュミレーションを作れないかと考えています。計算の使用用途が分からないのですが、多少時間を用いればクラメル自体の計算はできると思いますが、それでは不十分なのでしょうか？
- eatern27
- ベストアンサー率55% (635/1135)
>しかし、(3)式と(4)式を連立してdx1/dt、dx2/dtについて解いた式を使うとうまく計算できます。 ってのは、(3),(4)を手計算などで(dx/dtについて)『解析的に』解いて(つまり,dx/dtをx,yの関数として『具体的に』書き下して)、その式をプログラムに打ち込んでいるんですか？ もしそうなら、わざわざ『解析的に』解く必要性は何かあるのでしょうか？
- boobee0125
- ベストアンサー率72% (35/48)
私に答えられそうにもない高度な質問に口出ししてごめんなさい。 ちなみに(3)、(4)のままでルンゲクッタ法を使う場合は右辺の微分 dx~/dt = g~(z~) の中間値： g~(zi~+αj*k0~+βj*k1~+...) （z~は x~、y~ の合成ベクトル） はどのように算出、更新されているんでしょうか？これってエクセルの循環参照のような形になりませんか？ >連立して右辺に微分が含まれない形にすると１つの運動方程式の中に100の階乗個の項が入ってしまいます． 100の階乗個とは天文学的な数字ですが、これは100の２乗程度ではないんでしょうか？であればコンピューターブン回しの術もありえますが。しかし100個の多重振り子となるとほとんど連続体の世界ですから（可能ならば）偏微分方程式の数値解法の方が良さそうな気もしますね。
- boobee0125
- ベストアンサー率72% (35/48)
(1)、(2) はともかく、(3)、(4) は右辺に微分が含まれていてルンゲクッタ法を始めとする微分方程式のデジタルシミュレーション（あるいは状態方程式）の形式になっていませんのでムリな話ではないでしょうか？たとえ発散しなくても正しい結果は得られないと思います。なぜ dx1/dt、dx2/dtについて解いた式が使えないんでしょうか？
補足
ご意見ありがとうございます．最終的には振り子が100個繋がった多重振り子の計算をしなければなりません．その場合，100個の運動方程式が導かれ，連立して右辺に微分が含まれない形にすると１つの運動方程式の中に100の階乗個の項が入ってしまいます．ですので，(1)～(4)式のままで計算したいのです．おそらく，ルンゲクッタの式を変える必要があると思うのですが，どうしていいかわかりません．
関連するQ&A
- 常微分方程式、４次のルンゲクッタ法
(d^2x/dt^2)-2(dy/dt)=f(x) (d^2y/dt^2)+2(dx/dt)=g(y) この連立常微分方程式を４次のルンゲクッタ法で解くためにはどうすればいいのでしょうか？
- 締切済み
- 数学・算数
- ルンゲクッタで非線形連立微分方程式を解きたい
数値計算は素人です。 ルンゲクッタで下記の3本のようなの方程式系を解きたいのですが、 うまくいきません。 dx/dt=3x+(x-c)dz/dt dy/dt=y*(dz/dt - 2(x-1)/(x-c)) dz/dt=z* ((a/(x-c))(y/z)+ 2x^2 - 3b)/(y/z-x^2) a,b,cは定数です。 計算自体は途中で発散等せずにできているのようですが、 結果が間違っています。 ルンゲクッタのコード自体は、他の既知の関数で実験し、 誤りはないかことを確認したつもりです。 そもそも、関数系がこれくらい煩雑になれば、 ルンゲクッタでは解けないのでしょうか？ 年の暮れのお忙しい時期ですが、アドバイスお願いします。
- ベストアンサー
- 数学・算数
- 振り子の問題について （平衡点など
d^2x/dt^2 + a*dx/dt +sin(x)=0 この問題は非線形微分方程式です dx/dt=y dy/dt=-a*y -sin(x) とおいて連立させて平衡点を求めれば （0,0）（0,π）？になるとお思います この微分方程式を平衡点のまわりで線形化させたいのですが解法が全然分かりません アドバイスをいただけるか参考サイトを教えていただけけないでしょうか
- 締切済み
- 物理学
- 複素関数論
以下の証明ができません。 解答方針だけでもいいので、どなたか教えてください。 u1=u(x1(t),y1(t)) v1=v(x1(t),y1(t)) u2=u(x2(t),y2(t)) v2=v(x2(t),y2(t)) とし、 cosθ={(dx1/dt*dx2/d2)+(dy1/dt*dy2/dt)}/√[{(dx1/dt)^2+(dy1/dt)^2}+{(dx2/dt)^2+(dy2/dt)^2}] cosθ'={(du1/dt*du2/d2)+(dv1/dt*dv2/dt)}/√[{(du1/dt)^2+(dv1/dt)^2}+{(du2/dt)^2+(dv2/dt)^2}] とすると、u(x,y)とv(x,y)がコーシーリーマンの関係式を満たす時 |cosθ|＝|cosθ'| となることを示せ。
- ベストアンサー
- 数学・算数
- 陰関数そのものを使った積分の計算法
いろいろな曲線の表示において、微分や積分の計算法を整理してみました。 x^2+y^2=4上の点(x,y)=(1,√3)でのdy/dxの値の求め方。 陽関数。y=√(4-x^2)よりdy/dx=-x/√(4-x^2)。x=1のとき、dy/dx=-1/√3。 媒介変数。x=2cos(θ),y=2sin(θ)とすると、dy/dx=dy/dθ÷dx/dθ=-cos(θ)/sin(θ)。 θ=π/3のとき、dy/dx=-1/√3。 逆関数。x=√(4-y^2)よりdy/dx=1÷dx/dy=-√(4-y^2)/y。y=√3のとき、dy/dx=-1/√3。 極座標に変数変換。(x,y)→(r,θ) (ただし、x=rcos(θ),y=rsin(θ))とすると、(1,√3)→(2,π/3)。 x^2+y^2=4→r=2。dx=cos(θ)dr-rsin(θ)dθ、dy=sin(θ)dr+rcos(θ)dθ。dr/dθ=0。 よって、dy/dx=-cos(θ)/sin(θ)。θ=π/3のとき、dy/dx=-1/√3。 陰関数。2x+2y(dy/dx)=0より、dy/dx=-x/y=1/√3。 y≧0,x^2+y^2≦4の面積の求め方。 陽関数。境界はy=√(4-x^2)より∫[-2,2]ydx=∫[-2,2]√(4-x^2)dx=[(1/2)√(4-x^2)+2arcsin(x/2)] [-2,2] = 2π 媒介変数。境界をx=2cos(θ),y=2sin(θ)とすると、∫[-2,2]ydx=∫[π,0]2sin(θ){-2sin(θ)}dθ = 2π 逆関数。境界はx=√(4-y^2)より∫[-2,2]ydx=2∫[0,1]y(dx/dy)dy=2∫[2,0]y(-y/√(4-y^2))dy=2π 極座標に変数変換。(x,y)→(r,θ)(ただし、x=rcos(θ),y=rsin(θ))とすると、 [y≧0,x^2+y^2≦4]→[0≦r≦1,0≦θ≦π]、ヤコビアンはr。よって、 ∫[y≧0,x^2+y^2≦4]dxdy=∫[0≦r≦2,0≦θ≦π]rdrdθ=2π 以上のように計算法を比べてみると、陰関数そのものを使った積分の計算法を僕は知りません。 数学の理論はボタンをかけるように、パラレルな理論があると信じているのですが、 一方を知らないので気になります。 陰関数そのものを使った積分の計算法があれば教えていただけますようお願いいたします。
- ベストアンサー
- 数学・算数
- 三角関数微分の問題です
===================================================== 【問題】 (１) x=a(t-sin t) y=a(1-cos t) (a>0) (0 <= t <= 2π) dy/dxを求めよ。 (２) y=f(x)=sin(α arcsin x) f^(n) (0)を求めよ。 ↑ f(0)をn回微分したもの ======================================================== という問題で、(１)はなんとか解けたと思うのですが、(２)が行き詰ってしまいました。私の回答を載せさせてもらいますので、ご指摘や模範解答のほう宜しくお願いします。 =========================================================== 【自分の回答】 (1) dx/dt=a(1-cos t),dy/dt=a*sin t ∴dy/dx=(a*sin t)/{a(1-cos t)}=(sin t) /(1-cos t) (2) y'=1 / √(1- α^2 * sin^-2 x)=(sin x)/ √(sin^2 x - α^2) ∴y'*√(sin^2 x - α^2)/(sin x)=1 両辺をxについて微分し両辺√(sin^2 x - α^2)を掛けて整理すると、 y"*sin x +y'*α^2 * (cos x) / (sin x) =0 ⇒(1/α^2)* y" *(sin^2 x) /(cos x)+ y'=0 **************************************************** ここでライプニッツの定理や数学的帰納法を使って計算していくのですが、 f'(0),f"(0),f^(3) (0),..........といった感じに出来ません。 **************************************************** ===========================================================
- ベストアンサー
- 数学・算数
- 三角関数の微分
IIIＣをやってて少し気になったので 質問させてください “y=sin(3x) と表されるとき(dy/dx)を求めよ” という問題で私は２つの解答例が思い浮かびました [解答例１] u=3xと置くと (dy/du)=3 (du/dx)=u*cos(u) となり、合成関数の微分法の公式から (dy/dx)=(dy/du)*(du/dx) =(3)*{u*cos(u)} =3*3x*cos(3x) =9x*cos(3x) (答) [解答例２] ３倍角の公式から sin(3x) =3sin(x)-4{sin(x)}^3 よって (dy/dx) =[3sin(x)]'-[4{sin(x)}^3]' =3cos(x)-12[{sin(x)}^2]*[cos(x)] (答) となってしまい、同じ式を微分したのに 異なる解答が出てきます。 この場合どちらが正しいのでしょうか。 あるいはどちらも正しいのでしょうか。 回答をお願いします
- ベストアンサー
- 数学・算数
- 数III 微分
微分について、2点わからないことがあるので質問します。 問題1. x^2-y^2=a^2のとき、d^2/y/dx^2をx,yを用いてあらわせ。ただし、aは定数とする。 x^2-y^2=a^2の両辺をxで微分して、2x-2ydy/dx=0より、dy/dx＝x/y(y≠0)さらに d^2/y/dx^2=(d/dx)*(dy/dx)・・・(1)として計算することには納得できました。しかし (1)で{1*y-x*(dy/dx)}/y^2 となっているのは、ｙをｘの関数 ｙ＝±√(x^2-a^2)として書けるからでしょうか。ｙが定数にならない理由を教えてください。 問題2.(dy/dx)を求めよ、ただしaは0でない定数とする。 x=a(cost+tsint),y=a(sint-tcost)のとき、 dx/dt=atcost,dy/dt=atsint よって dy/dx=sint/costなぜ,tantにならないのでしょうか。ほかの問題の答えでも、dy/dx=-sinθ/cosθでした。tanθになおさない理由をおしえてください。お願いします。
- 締切済み
- 数学・算数
補足
ご意見ありがとうございます。中間値は今手計算してみているのですが、確かにエクセルの循環参照のようになります。このようにならないような方法をさがしています。 また、運動方程式の中に１００の階乗個の項が入るといいましたが、勘違いでした。すみません！ただ、クラメルの公式（クラメルの公式の条件を満たしていると仮定しています）を使って右辺に微分が含まれないような形にしようとしているのですが、その計算過程で１００の階乗個の計算をしなければならないと思うのです。 偏微分方程式の数値解法見てみます！