演習: 振り子のアニメーション

演習:振り子のアニメーション
担当:柳澤優香
「第８回演習問題」として,CourseN@viで提出.締め切り 6/10（金）13時まで
振り子運動
• 振り子とは,一点で支えられた糸や棒（伸び縮みしな
い）の先に小さいおもりをつけたもの
• 一度力を加えると揺れを繰り返し、摩擦や空気抵抗
がない場合永久に揺れ続ける．時計やメトロノームに
応用.
• 二重振り子は角度が大きい場合,カオス現象（初期値
に対する鋭敏な依存性）が観察される
1
単振り子の運動方程式
•
•
運動方程式より,重力を分解した力Fが復元力として単振動する．
赤矢印の向きは,半径lの円周の接線方向となり,円弧を描いて往復運動を行う
𝐹 = −𝑚𝑔sin𝜃
𝑚𝑎 = −𝑚𝑔sin𝜃
,
θ
l
𝑑 𝑠
𝑚 , = −𝑚𝑔sin𝜃
x
𝑑𝑡
s
𝑑, 𝜃
𝑔
= − sin𝜃 （ s=lθより）
,
𝑑𝑡
𝑙
l:糸の長さ
m:質量
θ：角度
t：経過時間
g：重力加速度（9.8）
θが微小ならx∼s
l
θ
糸の長さ：l
質量：m
mgsinθ
原点
θ
mg
2
振れ角θの周期を計算
・高次の微分方程式を１階の連立微分方程式に変換
𝑑𝜃
𝑦1 = 𝜃, 𝑦, =
𝑑𝑡 とすれば,
𝑦1 =𝑦,
𝑔
𝑦, = − sin(𝑦1 )
𝑙
function dy = furiko(t,y,flag,l)
g = 9.8;
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = - (g/l).*sin(y(1));
スクリプトで実行
ts=[0 20]; y0=[1;0];l=1.5;
[T,Y]=ode45('furiko',ts,y0,[],l)
;
plot(T,Y(:,1))
3
アニメーションの作成
1,糸の長さとθから, おもりのx座標とy座標を計算
2,糸の描画は,plot([0 x],[0 y])
3,描画範囲を設定
axis(plotsize);
plotsize = l*[-1 1 -1 0];
xの範囲
yの範囲
4,おもりの描画
plot(px,py,'o','MarkerSize',20)
5,少し動きを止めたい場合は
pauseを使う（ボール投げ参照）
6,1から5を繰り返す
(0,0)
θ
糸の長さ：l
質量：m
(x,y)
4
時間がある人は
応用:二重振り子の運動方程式
(0,0)
• ラグランジュの運動方程式より導出
θ
糸の長さ：l１
(x,y)
𝑚 1 + 𝑚 , 𝑙1
𝑙 1 𝑙,
67 8
7
69
67 8
6 97
+ 𝑚, 𝑙,
cos 𝜃 − 𝜑
67 :
6 97
7
,6 :
+𝑙 , 7
69
cos 𝜃 − 𝜑
6:
+𝑚 , 𝑙, ( ), 𝑠𝑖𝑛(𝜃 − 𝜑)+(𝑚1
69
質量：m１
+ 𝑚 , ) 𝑔𝑠𝑖𝑛(𝜃)=0
68
− 𝑙 1 𝑙, ( ), 𝑠𝑖𝑛(𝜃 − 𝜑)+ 𝑔𝑙 , 𝑠𝑖𝑛(𝜑)=0
69
長い。。。
φ
糸の長さ：l２
(x2,y2)
質量：m２
5
まずはこちらでやってみよう!
• 𝜃 − 𝜑が著しく０に近いと考えると,
𝑠𝑖𝑛 𝜃 − 𝜑 = 0, cos 𝜃 − 𝜑 = 1,
𝑚 1 + 𝑚 , 𝑙1
67 8
6 97
+ 𝑚, 𝑙,
67 :
6 97
𝑠𝑖𝑛 𝜃 ~𝜃
とでき,
+(𝑚 1 + 𝑚 ,) 𝑔𝜃=0
67 8
7
,6 :
𝑙 1 𝑙, 7 +𝑙 ,
+𝑔𝑙 , 𝜑=0
69
6 97
𝜃 = 𝑦,
68
69
= 𝑦, , 𝜑=𝑦@ ,
6:
69
= 𝑦A
とすれば,
𝑦1 =𝑦,
𝑚1 + 𝑚, 𝑙1𝑦,+𝑚, 𝑙, 𝑦A=- (𝑚1 + 𝑚,) g𝑦1
𝑦@ =𝑦A
𝑙1 𝑦,+𝑙, 𝑦A =- g𝑦@
6
二重振り子のアニメーション作成
１：振れ角の周期を計算
-
４元一次方程式を解き,y1とy3を求める
×
=
function dy = dblfuriko(t,y,flag,l1,m1,l2,m2)
g = 9.8;
F = zeros(4); %4×4のゼロ行列
F(1,1) = 1;
F(2,2) = (m1 + m2)*l1; F(2,4) = m2*l2;
F(3,3) = 1;
F(4,2) = l1; F(4,4) = l2;
b = zeros(4,1);
b(1) = y(2);
b(2) = -(m1 + m2)*g*y(1);
b(3) = y(4);
b(4) = -g*y(3);
dy = F¥b; %ガウスの消去法
２：アニメーションの作成
y1とy3を基に作成．単振り子と同様の方法
7
参考:二重振り子のアニメーション
dbl_furiko.m
ts=[0 10]; y0=[0.1 0 0.1 1]';
l1=1; m1=8; l2=1; m2=1;
[T,Y]=ode45('dblfuriko',ts,y0,[],l1,m1,l2,m2);
dblfurikoplot(T,Y,l1,l2) %角度の周期を計算
plotsize = 1.2*(l1 + l2)*[-1 1 -1 0]; % 描画範囲
for i=1:length(T)-1
px1 = l1*sin(Y(i,1)); % ボール1のX座標
py1 = -l1*cos(Y(i,1)); % ボール1のY座標
px2 = px1 + l2*sin(Y(i,3)); % ボール2のX座標
py2 = py1 - l2*cos(Y(i,3)); % ボール2のY座標
hold on
plot([0 px1],[0 py1]) % 糸1の描画
plot(px1,py1,'o','MarkerSize',20) % ボール1の描画
plot([px1 px2],[py1 py2]) % 糸1の描画
plot(px2,py2,'o','MarkerSize',20) % ボール1の描画
axis(plotsize); % 描画範囲の設定
pause(0.01)
hold off
plot(0)
end
8
参考文献
• 青山貴伸,蔵本一峰,森口肇:「最新 使え
る!MATLAB 」
9