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

ライブスクリプトを開く

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

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

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

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

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

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

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

$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))$

$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を使用して二重振子のアニメーションを再生します。