最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

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

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

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

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

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

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

  • 1 つ目の玉の角度位置 θ1(t)

  • 2 つ目の玉の角度位置 θ2(t)

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

  • 1 つ目の棒の長さ L1

  • 2 つ目の棒の長さ L2

  • 1 つ目の玉の質量 m1

  • 2 つ目の玉の質量 m2

  • 重力定数 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 つ目の棒の張力を T1 として指定し、2 つ目の棒の張力を T2 として指定します。

syms T_1 T_2

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

m1 に作用する力を評価します。水平方向と垂直方向の力の成分を均衡化して 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 = 

-m1L1sin(θ1(t))t θ1(t)2-L1cos(θ1(t))2t2 θ1(t)=T2sin(θ2(t))-T1sin(θ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 = 

m1L1sin(θ1(t))2t2 θ1(t)+L1cos(θ1(t))t θ1(t)2=T1cos(θ1(t))-gm1-T2cos(θ2(t))

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

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

-m2L1sin(θ1(t))t θ1(t)2+L2sin(θ2(t))t θ2(t)2-L1cos(θ1(t))2t2 θ1(t)-L2cos(θ2(t))2t2 θ2(t)=-T2sin(θ2(t))

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

m2L1cos(θ1(t))t θ1(t)2+L2cos(θ2(t))t θ2(t)2+L1sin(θ1(t))2t2 θ1(t)+L2sin(θ2(t))2t2 θ2(t)=T2cos(θ2(t))-gm2

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

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

運動方程式には、θ1θ2T1、および T2 の 4 つの未知数があります。T1T2 の 2 つの未知数を eqx_1eqy_1 から評価します。関数solveを使用して T1T2 を求めます。

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

T1T2 の解を eqx_2eqy_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/s2) を 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 = 

cos(θ1(t))σ1-3sin(θ2(t))t θ2(t)22-sin(θ1(t))t θ1(t)2+3cos(θ2(t))2t2 θ2(t)2=-2sin(θ2(t))cos(θ1(t))2σ1+sin(θ1(t))2σ1+49sin(θ1(t))5cos(θ1(t))sin(θ2(t))-cos(θ2(t))sin(θ1(t))where  σ1=2t2 θ1(t)

eqn_2 = subs(eqRed_2)
eqn_2 = 

cos(θ1(t))t θ1(t)2+3cos(θ2(t))t θ2(t)22+sin(θ1(t))σ1+3sin(θ2(t))2t2 θ2(t)2=2cos(θ2(t))cos(θ1(t))2σ1+sin(θ1(t))2σ1+49sin(θ1(t))5cos(θ1(t))sin(θ2(t))-cos(θ2(t))sin(θ1(t))-495where  σ1=2t2 θ1(t)

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

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

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

S
S = 

(θ2Dtheta2θ1Dtheta1)

次に、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)')

手順 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;

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