mathematicaを用いて振り子の支点が移動（ひもの長さが変化）するときの挙動を導出したい

点Oを支点として振り子運動をしているが点Aにぶつかると、点Aを支点として振り子運動をする。
点Aを支点として振り子運動しているが点Aから離れると、点Oを支点として振り子運動をする。（下図参照）
このときの振り子の振れる角度および速度を導出、また現象のアニメーションも作成したいと考えています。
現在参考書等を用いて勉強しているところでありますが、なかなか解決できません。
知識不足で申し訳ないのですが、ご教授願います。

質問者からの補足コメント

• 

  最下点で右方向に初速v0を与えたときはどうすればよいのでしょうか…（下図参照）
  エラーが起きてしまいうまくいきません…

  
    補足日時：2016/09/02 13:21

• 

  回答No.1＆No.2の方へ
  作成して下さったプログラム内のk1、k2はそれぞれg/l1、g/l2ですよね？（l1、l2はひもの長さ）
  つまり、k1＝2.0 > k2=1.0よりl1 < l2ですよね？
  しかしNDSolve内の運動方程式が
  k1HeavisideTheta[y[t]]+k2HeavisideTheta[-y[t]]
  とひもの短い方を正としています。
  しかし、アニメーションではひもの長い方からスタートしています。。。
  なせでしょうか。。。

    補足日時：2016/09/04 02:17

No.2 回答者： sanze_hui 回答日時： 2016/09/03 11:00

#1でNDSolveに与えている初期条件

y'[0] == 0, y[0] == y0
を変えれば良いです。

最下点にいるのであればy[0] == 0ですね。
y'[0]はv0から求められますので、それを入れてください。

この回答へのお礼

NDSolveに与えている初期条件
y’[0]==a , y[0]==0
a をいろんな値にしてもエラーが起きてしまいます。。。

お礼日時：2016/09/03 12:20

No.1 回答者： sanze_hui 回答日時： 2016/08/31 21:38

> 振り子の振れる角度および速度を導出

これは位置エネルギーと運動エネルギーの保存を考えれば解けます

> 現象のアニメーションも作成
振れ角が微小ではないなら、微分方程式を知っている必要があります。この問題の運動方程式が下記の形になります。
「y(t) > 0 のとき y''(t) = k1 sin(y(t)), y(t) < 0 のときy''(t) = k2 sin(y(t)) (ただしk1 > 0, k2 > 0)」

Mathematicaでこれを解く手順としては、下記のようになります。
下記では適当にk1, k2, 初期条件y(0)を設定していますが、問題に合わせて設定してください。

k1 = 2.0; k2 = 1.0; y0 = 45 Degree;
ans = y /. NDSolve[{y''[t] == -(k1 HeavisideTheta[y[t]] + k2 HeavisideTheta[-y[t]]) Sin[y[t]], y'[0] == 0, y[0] == y0}, y, {t, 0, 20}][[1]];
Plot[ans[t], {t, 0, 20}]

条件分岐をHeavisideThetaで表現するのが良いと思います。

アニメーションを作るなら、グラフから周期を求めて、Animateを使うと作成できます。

tend = t /. FindMaximum[ans[t], {t, 11}][[2]];
Animate[Graphics[{
Circle[{0, 0}, k1], Circle[{0, k2 - k1}, k2]
, Red,
If[ans[t] >= 0,
{Line[{{0, 0}, #}], Disk[#, 0.1]} &[{k1 Sin[ans[t]], -k1 Cos[ans[t]]}],
{Line[{{0, 0}, {0, k2 - k1}, #}], Disk[#, 0.1]} &[{k2 Sin[ans[t]], -k2 Cos[ans[t]] + k2 - k1}]
]
},
PlotRange -> {{-k1 - k2, k1 + k2}, {-k1 - k2, k1 + k2}}
], {t, 0, tend}]

この回答へのお礼

実際にmathematica上でたしかめました。
思い通りのことができました。
ご丁寧にありがとうございました。
非常に助かりました。

お礼日時：2016/09/02 11:27