二重振子運動のアニメーションと解

この例では、MATLAB® と Symbolic Math Toolbox™ を使用して二重振子の動きをモデル化する方法を示します。

二重振子の運動方程式を解き、アニメーションを作成して二重振子運動をモデル化します。

手順 1: 二重振子の質量の変位、速度、加速度の定義

次の図は二重振子のモデルを示したものです。二重振子は 2 つの振り玉と 2 つの剛体棒で構成されます。

Double pendulum consisting of two bobs with masses m1 and m2, connected by two rigid rods with lengths L1 and L2. The angles between the first and second bobs with respect to the vertical line are theta1 and theta2.

状態変数を定義して二重振子の動きを記述します。

1 つ目の玉の角度位置 $θ_{1} (t)$
2 つ目の玉の角度位置 $θ_{2} (t)$

変数を定義して二重振子のプロパティを記述します。

1 つ目の棒の長さ $L_{1}$
2 つ目の棒の長さ $L_{2}$
1 つ目の玉の質量 $m_{1}$
2 つ目の玉の質量 $m_{2}$
重力定数 $g$

簡単にするために、2 つの剛体棒の質量は無視します。変数はすべてsymsを使用して指定します。

syms theta_1(t) theta_2(t) L_1 L_2 m_1 m_2 g

直交座標における二重振子の変位を定義します。

x_1 = L_1*sin(theta_1);
y_1 = -L_1*cos(theta_1);
x_2 = x_1 + L_2*sin(theta_2);
y_2 = y_1 - L_2*cos(theta_2);

関数diffを使用して、変位を時間で微分することで速度を求めます。

vx_1 = diff(x_1);
vy_1 = diff(y_1);
vx_2 = diff(x_2);
vy_2 = diff(y_2);

速度を時間で微分することで加速度を求めます。

ax_1 = diff(vx_1);
ay_1 = diff(vy_1);
ax_2 = diff(vx_2);
ay_2 = diff(vy_2);

手順 2: 運動方程式の定義

ニュートンの法則に基づいて運動方程式を定義します。

最初に、1 つ目の棒の張力を $T_{1}$ として指定し、2 つ目の棒の張力を $T_{2}$ として指定します。

syms T_1 T_2

次に、両方の質量に作用する力の自由体図を作成します。

Free-body diagram of the forces that act on the first bob

$m_{1}$ に作用する力を評価します。水平方向と垂直方向の力の成分を均衡化して 1 つ目の玉の運動方程式を定義します。それらの 2 つの方程式をシンボリック方程式 eqx_1 および eqy_1 として指定します。

eqx_1 = m_1*ax_1(t) == -T_1*sin(theta_1(t)) + T_2*sin(theta_2(t))

eqx_1 = 
 $- m_{1} (L_{1} \sin (θ_{1} (t)) {(\frac{\partial}{\partial t} θ_{1} (t))}^{2} - L_{1} \cos (θ_{1} (t)) \frac{\partial^{2}}{\partial t^{2}} θ_{1} (t)) = T_{2} \sin (θ_{2} (t)) - T_{1} \sin (θ_{1} (t))$

eqy_1 = m_1*ay_1(t) == T_1*cos(theta_1(t)) - T_2*cos(theta_2(t)) - m_1*g

eqy_1 = 
 $m_{1} (L_{1} \sin (θ_{1} (t)) \frac{\partial^{2}}{\partial t^{2}} θ_{1} (t) + L_{1} \cos (θ_{1} (t)) {(\frac{\partial}{\partial t} θ_{1} (t))}^{2}) = T_{1} \cos (θ_{1} (t)) - g m_{1} - T_{2} \cos (θ_{2} (t))$

Free-body diagram of the forces that act on the second bob

$m_{2}$ に作用する力を評価します。水平方向と垂直方向の力の成分を均衡化して 2 つ目の玉の運動方程式を定義します。それらの 2 つの方程式をシンボリック方程式 eqx_2 および eqy_2 として指定します。

eqx_2 = m_2*ax_2(t) == -T_2*sin(theta_2(t))

eqx_2 = 
 $- m_{2} (L_{1} \sin (θ_{1} (t)) {(\frac{\partial}{\partial t} θ_{1} (t))}^{2} + L_{2} \sin (θ_{2} (t)) {(\frac{\partial}{\partial t} θ_{2} (t))}^{2} - L_{1} \cos (θ_{1} (t)) \frac{\partial^{2}}{\partial t^{2}} θ_{1} (t) - L_{2} \cos (θ_{2} (t)) \frac{\partial^{2}}{\partial t^{2}} θ_{2} (t)) = - T_{2} \sin (θ_{2} (t))$

eqy_2 = m_2*ay_2(t) == T_2*cos(theta_2(t)) - m_2*g

eqy_2 = 
 $m_{2} (L_{1} \cos (θ_{1} (t)) {(\frac{\partial}{\partial t} θ_{1} (t))}^{2} + L_{2} \cos (θ_{2} (t)) {(\frac{\partial}{\partial t} θ_{2} (t))}^{2} + L_{1} \sin (θ_{1} (t)) \frac{\partial^{2}}{\partial t^{2}} θ_{1} (t) + L_{2} \sin (θ_{2} (t)) \frac{\partial^{2}}{\partial t^{2}} θ_{2} (t)) = T_{2} \cos (θ_{2} (t)) - g m_{2}$

手順 3: 力の評価とシステム方程式の簡約

二重振子の運動を 4 つの運動方程式で記述しています。棒に作用する力を評価し、一連の 4 つの方程式を 2 つの方程式に簡約します。

運動方程式には、 $θ_{1}$ 、 $θ_{2}$ 、 $T_{1}$ 、および $T_{2}$ の 4 つの未知数があります。 $T_{1}$ と $T_{2}$ の 2 つの未知数を eqx_1 と eqy_1 から評価します。solve関数を使用して $T_{1}$ と $T_{2}$ を求めます。

Tension = solve([eqx_1 eqy_1],[T_1 T_2]);

$T_{1}$ と $T_{2}$ の解を eqx_2 と eqy_2 に代入します。

eqRed_1 = subs(eqx_2,[T_1 T_2],[Tension.T_1 Tension.T_2]);
eqRed_2 = subs(eqy_2,[T_1 T_2],[Tension.T_1 Tension.T_2]);

簡約した 2 つの方程式で振子運動が完全に記述されます。

手順 4: システム方程式の求解

システム方程式を解いて振子運動を記述します。

最初に、質量 ( $kg$ )、棒の長さ ( $m$ )、および重力 ( $m / s^{2}$ ) を SI 単位で定義します。それらの値を 2 つの簡約した方程式に代入します。

L_1 = 1;
L_2 = 1.5;
m_1 = 2;
m_2 = 1;
g = 9.8;
eqn_1 = subs(eqRed_1)

eqn_1 = 
 $\begin{array}{l} \cos (θ_{1} (t)) σ_{1} - \frac{3 \sin (θ_{2} (t)) {(\frac{\partial}{\partial t} θ_{2} (t))}^{2}}{2} - \sin (θ_{1} (t)) {(\frac{\partial}{\partial t} θ_{1} (t))}^{2} + \frac{3 \cos (θ_{2} (t)) \frac{\partial^{2}}{\partial t^{2}} θ_{2} (t)}{2} = - \frac{2 \sin (θ_{2} (t)) ({\cos (θ_{1} (t))}^{2} σ_{1} + {\sin (θ_{1} (t))}^{2} σ_{1} + \frac{49 \sin (θ_{1} (t))}{5})}{\cos (θ_{1} (t)) \sin (θ_{2} (t)) - \cos (θ_{2} (t)) \sin (θ_{1} (t))} \\ where \\ σ_{1} = \frac{\partial^{2}}{\partial t^{2}} θ_{1} (t) \end{array}$

eqn_2 = subs(eqRed_2)

eqn_2 = 
 $\begin{array}{l} \cos (θ_{1} (t)) {(\frac{\partial}{\partial t} θ_{1} (t))}^{2} + \frac{3 \cos (θ_{2} (t)) {(\frac{\partial}{\partial t} θ_{2} (t))}^{2}}{2} + \sin (θ_{1} (t)) σ_{1} + \frac{3 \sin (θ_{2} (t)) \frac{\partial^{2}}{\partial t^{2}} θ_{2} (t)}{2} = \frac{2 \cos (θ_{2} (t)) ({\cos (θ_{1} (t))}^{2} σ_{1} + {\sin (θ_{1} (t))}^{2} σ_{1} + \frac{49 \sin (θ_{1} (t))}{5})}{\cos (θ_{1} (t)) \sin (θ_{2} (t)) - \cos (θ_{2} (t)) \sin (θ_{1} (t))} - \frac{49}{5} \\ where \\ σ_{1} = \frac{\partial^{2}}{\partial t^{2}} θ_{1} (t) \end{array}$

2 つの方程式は非線形 2 階微分方程式です。これらの方程式を解くには、関数odeToVectorFieldを使用して 1 階微分方程式に変換します。

[V,S] = odeToVectorField(eqn_1,eqn_2);

ベクトル V の要素は、S の要素の時間微分に等しい 1 階微分方程式を表します。S の要素は、状態変数 $θ_{2}$ 、 $d θ_{2} / dt$ 、 $θ_{1}$ 、および $d θ_{1} / dt$ です。これらの状態変数で二重振子の角変位と角速度を記述します。

S = 
 $(\begin{array}{c} θ_{2} \\ {Dtheta}_{2} \\ θ_{1} \\ {Dtheta}_{1} \end{array})$

次に、1 階微分方程式をハンドル M をもつ MATLAB 関数に変換します。

M = matlabFunction(V,Vars={'t','Y'});

状態変数の初期条件を [pi/4 0 pi/6 0] と定義します。関数ode45を使用して状態変数の解を求めます。解は区間 [0 10] の時間の関数になります。

initCond = [pi/4 0 pi/6 0];
sols = ode45(M,[0 10],initCond);

状態変数の解をプロットします。

plot(sols.x,sols.y)
legend("\theta_2","d\theta_2/dt","\theta_1","d\theta_1/dt")
title("Solutions of State Variables")
xlabel("Time (s)")
ylabel("Solutions (rad or rad/s)")

$Figure contains an axes object. The axes object with title Solutions of State Variables, xlabel Time (s), ylabel Solutions (rad or rad/s) contains 4 objects of type line. These objects represent \theta_2, d\theta_2/dt, \theta_1, d\theta_1/dt.$

手順 5: 振動する二重振子のアニメーションの作成

振動する二重振子のアニメーションを作成します。

最初に、両方の振子の座標を解 sols から deval を使用して評価する 4 つの関数を作成します。

x_1 = @(t) L_1*sin(deval(sols,t,3));
y_1 = @(t) -L_1*cos(deval(sols,t,3));
x_2 = @(t) L_1*sin(deval(sols,t,3))+L_2*sin(deval(sols,t,1));
y_2 = @(t) -L_1*cos(deval(sols,t,3))-L_2*cos(deval(sols,t,1));

次に、関数fanimatorを使用して、1 つ目の振り玉のストップモーションアニメーションオブジェクトを作成します。既定では、fanimator を使用すると、単位時間あたりのフレーム数を 10 として t の範囲が 0 から 10 までのアニメーションオブジェクトが作成されます。関数 plot を使用して座標をプロットします。x 軸と y 軸が同じ長さになるように設定します。

fanimator(@(t) plot(x_1(t),y_1(t),"ro",MarkerSize=m_1*10,MarkerFaceColor="r"));
axis equal;

次に、1 つ目の剛体棒、2 つ目の振り玉、および 2 つ目の剛体棒のアニメーションオブジェクトを追加します。

hold on;
fanimator(@(t) plot([0 x_1(t)],[0 y_1(t)],"r-"));
fanimator(@(t) plot(x_2(t),y_2(t),"go",MarkerSize=m_2*10,MarkerFaceColor="g"));
fanimator(@(t) plot([x_1(t) x_2(t)],[y_1(t) y_2(t)],"g-"));

関数textを使用して、経過時間をカウントするテキストを追加します。num2strを使用して時間パラメーターを文字列に変換します。

fanimator(@(t) text(-0.3,0.3,"Timer: "+num2str(t,2)));
hold off;

Figure contains an axes object. The axes object contains 5 objects of type line, text. One or more of the lines displays its values using only markers

コマンドplayAnimationを使用して二重振子のアニメーションを再生します。

Animation of the double pendulum motion. The timer runs from 0 to 10 seconds.

参考

関数

diff | subs | solve | matlabFunction | fanimator | playAnimation | ode45