振子の周期的揺れの動きのシミュレーション

この例では、Symbolic Math Toolbox™ を使用して単純な振子の動きをモデル化する方法を示します。振子の運動方程式を導出し、その方程式を小角について解析的に解き、任意の角度について数値的に解きます。

手順 1: 運動方程式の導出

振子は、微分方程式に従う単純な機械的システムです。振子は、初期状態では垂直位置で停止しています。振子がある角度 $θ$ だけ動かされ放されると、重力によりその静止位置へ引っ張られます。その運動量により行き過ぎて角度 $- θ$ に達し (摩擦力がない場合)、これを繰り返します。重力による振子の動きに伴う復元力は $- m g \sin θ$ です。したがって、ニュートンの第 2 法則に従って、質量と加速度の積は $- m g \sin θ$ に等しくなければなりません。

syms m a g theta(t)
eqn = m*a == -m*g*sin(theta)

eqn(t) =  $a m = - g m \sin (θ (t))$

長さが $r$ の振子では、振子の重りの加速度は角加速度に $r$ を乗じたものです。

$a = r \frac{d^{2} θ}{d t^{2}}$ .

subs を $a$ の代わりに使用します。

syms r
eqn = subs(eqn,a,r*diff(theta,2))

eqn(t) = 
 $m r \frac{\partial^{2}}{\partial t^{2}} θ (t) = - g m \sin (θ (t))$

isolate を使用して eqn 内の角加速度を分離します。

eqn = isolate(eqn,diff(theta,2))

eqn = 
 $\frac{\partial^{2}}{\partial t^{2}} θ (t) = - \frac{g \sin (θ (t))}{r}$

定数 $g$ および $r$ を "固有振動数" としても知られる単一のパラメーターにまとめます。

$ω_{0} = \sqrt{\frac{g}{r}}$ .

syms omega_0
eqn = subs(eqn,g/r,omega_0^2)

eqn = 
 $\frac{\partial^{2}}{\partial t^{2}} θ (t) = - {ω_{0}}^{2} \sin (θ (t))$

手順 2: 運動方程式の線形化

この運動方程式は非線形であるため解析的な求解が困難です。角度が小さいと仮定して、 $\sin θ$ のテイラー展開を使用して方程式を線形化します。

syms x
approx = taylor(sin(x),x,'Order',2);
approx = subs(approx,x,theta(t))

approx =  $θ (t)$

運動方程式は線形方程式になります。

eqnLinear = subs(eqn,sin(theta(t)),approx)

eqnLinear = 
 $\frac{\partial^{2}}{\partial t^{2}} θ (t) = - {ω_{0}}^{2} θ (t)$

手順 3: 運動方程式の解析的な求解

方程式 eqnLinear を dsolve を使用して解きます。初期条件を 2 番目の引数として指定します。assume を使用して $ω_{0}$ が実数であると仮定して、解を単純化します。

syms theta_0 theta_t0
theta_t = diff(theta);
cond = [theta(0) == theta_0, theta_t(0) == theta_t0];
assume(omega_0,'real')
thetaSol(t) = dsolve(eqnLinear,cond)

thetaSol(t) = 
 $θ_{0} \cos (ω_{0} t) + \frac{θ_{t0} \sin (ω_{0} t)}{ω_{0}}$

手順 4: $ω_{0}$ の物理的意味

項 $ω_{0} t$ は、"位相" と呼ばれます。余弦関数と正弦関数は、 $2 π$ ごとに繰り返します。 $ω_{0} t$ を $2 π$ だけ変化させるのに必要な時間は時間周期と呼ばれます。

$T = \frac{2 π}{ω_{0}} = 2 π \sqrt{\frac{r}{g}}$ .

時間周期 $T$ は振子の長さの平方根に比例し、質量には依存しません。線形運動方程式では、時間周期は初期条件に依存しません。

手順 5: 振子運動のプロット

小角近似における振子の運動をプロットします。

物理パラメーターを定義します。

重力加速度 $g = 9.81 {m/s}^{2}$
振子の長さ $r = 1 m$

gValue = 9.81;
rValue = 1;
omega_0Value = sqrt(gValue/rValue);
T = 2*pi/omega_0Value;

初期条件を設定します。

theta_0Value  = 0.1*pi; % Solution only valid for small angles.
theta_t0Value = 0;      % Initially at rest.

物理パラメーターおよび初期条件を一般解に代入します。

vars   = [omega_0      theta_0      theta_t0];
values = [omega_0Value theta_0Value theta_t0Value];
thetaSolPlot = subs(thetaSol,vars,values);

調和振子運動をプロットします。

fplot(thetaSolPlot(t*T)/pi, [0 5]);
grid on;
title('Harmonic Pendulum Motion');
xlabel('t/T');
ylabel('\theta/\pi');

Figure contains an axes object. The axes object with title Harmonic Pendulum Motion, xlabel t/T, ylabel theta / pi contains an object of type functionline.

$θ (t)$ の解を求めた後、振子の運動を可視化します。

x_pos = sin(thetaSolPlot);
y_pos = -cos(thetaSolPlot);
fanimator(@fplot,x_pos,y_pos,'ko','MarkerFaceColor','k','AnimationRange',[0 5*T]);
hold on;
fanimator(@(t) plot([0 x_pos(t)],[0 y_pos(t)],'k-'),'AnimationRange',[0 5*T]);
fanimator(@(t) text(-0.3,0.3,"Timer: "+num2str(t,2)+" s"),'AnimationRange',[0 5*T]);

Figure contains an axes object. The axes object contains 3 objects of type parameterizedfunctionline, line, text.

コマンド playAnimation を入力して振子運動のアニメーションを再生します。

手順 6: 定エネルギー軌道を使用した非線形振子運動の判別

振子の非線形運動を理解するには、全エネルギーの方程式を使用して振子の軌道を可視化します。全エネルギーは保存されます。

$E = \frac{1}{2} m r^{2} {(\frac{d θ}{dt})}^{2} + m g r (1 - \cos θ)$

三角恒等式 $1 - \cos θ = 2 \sin^{2} (θ / 2)$ および関係 $ω_{0} = \sqrt{g / r}$ を使用して、スケーリングされたエネルギーを書き換えます。

$\frac{E}{m r^{2}} = \frac{1}{2} [{(\frac{d θ}{dt})}^{2} + {(2 ω_{0} \sin \frac{θ}{2})}^{2}]$

エネルギーは保存されるため、振子の運動は位相空間内の定エネルギー軌道によって記述できます。位相空間は座標 $θ$ および $d θ / d t$ をもつ抽象空間です。これらの軌道を fcontour を使用して可視化します。

syms theta theta_t omega_0
E(theta, theta_t, omega_0) = (1/2)*(theta_t^2+(2*omega_0*sin(theta/2))^2);
Eplot(theta, theta_t) = subs(E,omega_0,omega_0Value);

figure;
fc = fcontour(Eplot(pi*theta, 2*omega_0Value*theta_t), 2*[-1 1 -1 1], ...
    'LineWidth', 2, 'LevelList', 0:5:50, 'MeshDensity', 1+2^8);
grid on;
title('Constant Energy Contours in Phase Space ( \theta vs. \theta_t )');
xlabel('\theta/\pi');
ylabel('\theta_t/2\omega_0');

Figure contains an axes object. The axes object with title Constant Energy Contours in Phase Space ( blank theta blank vs. blank theta indexOf t baseline blank ), xlabel theta / pi, ylabel theta indexOf t/2 baseline omega indexOf 0 baseline contains an object of type functioncontour.

一定エネルギーの等高線は $θ$ 軸および $d θ / d t$ 軸に対して対称であり、 $θ$ 軸に沿って周期的です。2 つの動作領域を図に示します。

等高線図の下側のエネルギーは自身に近くなります。振子は 2 つの最大角と最大速度の間で揺れます。振子の運動エネルギーは重力エネルギーに打ち勝つには十分ではなく、振子を完全に回転させることはできません。

等高線図の上側のエネルギーは自身に近づきません。振子は常に 1 つの角度方向で運動します。振子の運動エネルギーは重力エネルギーに打ち勝つのに十分であり、振子を完全に回転させることができます。

手順 7:非線形運動方程式の求解

非線形運動方程式は、2 階微分方程式です。ode45 ソルバーを使用して、これらの方程式を数値的に解きます。ode45 は 1 次の系だけを受け入れるため、系を 1 次の系に簡約します。次に、ode45 への入力となる関数ハンドルを生成します。

2 階 ODE を 1 階 ODE 系として書き換えます。

syms theta(t) theta_t(t) omega_0
eqs = [diff(theta)   == theta_t;
       diff(theta_t) == -omega_0^2*sin(theta)]

eqs(t) = 
 $(\begin{array}{c} \frac{\partial}{\partial t} θ (t) = θ_{t} (t) \\ \frac{\partial}{\partial t} θ_{t} (t) = - {ω_{0}}^{2} \sin (θ (t)) \end{array})$

eqs  = subs(eqs,omega_0,omega_0Value);
vars = [theta, theta_t];

系の質量行列 M および方程式 F の右辺を含むベクトルを求めます。

[M,F] = massMatrixForm(eqs,vars)

M = 
 $(\begin{array}{c} 1 & 0 \\ 0 & 1 \end{array})$

F = 
 $(\begin{array}{c} θ_{t} (t) \\ - \frac{981 \sin (θ (t))}{100} \end{array})$

M および F は次の形で表されます。

$M (t, x (t)) \frac{dx}{dt} = F (t, x (t)) .$

後の計算を単純化するために、系を次の形に書き換えます。 $d x / d t = f (t, x (t))$ 。

f = M\F

f = 
 $(\begin{array}{c} θ_{t} (t) \\ - \frac{981 \sin (θ (t))}{100} \end{array})$

odeFunction を使用して、f を MATLAB 関数ハンドルに変換します。結果の関数ハンドルは MATLAB ODE ソルバー ode45 の入力になります。

f = odeFunction(f, vars)

f = function_handle with value:
    @(t,in2)[in2(2,:);sin(in2(1,:)).*(-9.81e+2./1.0e+2)]

手順 8:閉じたエネルギー等高線における運動方程式の求解

ode45 を使用して、閉じたエネルギー等高線の ODE を解きます。

エネルギー等高線図から、閉じた等高線は条件 $θ_{0} = 0$ 、 $θ_{t 0} / 2 ω_{0} \leq 1$ を満たします。 $θ$ および $d θ / d t$ の初期条件を変数 x0 に保存します。

x0 = [0; 1.99*omega_0Value];

解を求めるため、0 秒から 10 秒までの時間間隔を指定します。この間隔で、振子は 2 周期分動くことができます。

tInit  = 0;
tFinal = 10;

ODE を解きます。

sols = ode45(f,[tInit tFinal],x0)

sols = struct with fields:
     solver: 'ode45'
    extdata: [1x1 struct]
          x: [0 3.2241e-05 1.9344e-04 9.9946e-04 0.0050 0.0252 0.1259 0.3449 0.6020 0.8591 1.1161 1.3597 1.5996 1.8995 2.2274 2.4651 2.7028 2.9567 3.2138 3.4709 3.7150 3.9511 4.2483 4.5759 4.8239 5.0719 5.3182 5.5764 5.8346 6.0803 ... ] (1x45 double)
          y: [2x45 double]
      stats: [1x1 struct]
      idata: [1x1 struct]

sols.y(1,:) は角変位 $θ$ を表し、sols.y(2,:) は角速度 $d θ / d t$ を表します。

閉じた軌道の解をプロットします。

figure;

yyaxis left;
plot(sols.x, sols.y(1,:), '-o');
ylabel('\theta (rad)');

yyaxis right;
plot(sols.x, sols.y(2,:), '-o');
ylabel('\theta_t (rad/s)');

grid on;
title('Closed Path in Phase Space');
xlabel('t (s)');

Figure contains an axes object. The axes object with title Closed Path in Phase Space, xlabel t (s), ylabel theta indexOf t baseline blank (rad/s) contains 2 objects of type line.

振子の運動を可視化します。

x_pos = @(t) sin(deval(sols,t,1));
y_pos = @(t) -cos(deval(sols,t,1));
figure;
fanimator(@(t) plot(x_pos(t),y_pos(t),'ko','MarkerFaceColor','k'));
hold on;
fanimator(@(t) plot([0 x_pos(t)],[0 y_pos(t)],'k-'));
fanimator(@(t) text(-0.3,1.5,"Timer: "+num2str(t,2)+" s"));

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

コマンド playAnimation を入力して振子運動のアニメーションを再生します。

手順 9:開いたエネルギー等高線の解

ode45 を使用して、開いたエネルギー等高線の ODE を解きます。エネルギー等高線図から、開いた等高線は条件 $θ_{0} = 0$ 、 $θ_{t 0} / 2 ω_{0} > 1$ を満たします。

x0 = [0; 2.01*omega_0Value];
sols = ode45(f, [tInit, tFinal], x0);

開いたエネルギー等高線の解をプロットします。

figure;

yyaxis left;
plot(sols.x, sols.y(1,:), '-o');
ylabel('\theta (rad)');

yyaxis right;
plot(sols.x, sols.y(2,:), '-o');
ylabel('\theta_t (rad/s)');

grid on;
title('Open Path in Phase Space');
xlabel('t (s)');

Figure contains an axes object. The axes object with title Open Path in Phase Space, xlabel t (s), ylabel theta indexOf t baseline blank (rad/s) contains 2 objects of type line.

振子の運動を可視化します。

x_pos = @(t) sin(deval(sols,t,1));
y_pos = @(t) -cos(deval(sols,t,1));
figure;
fanimator(@(t) plot(x_pos(t),y_pos(t),'ko','MarkerFaceColor','k'));
hold on;
fanimator(@(t) plot([0 x_pos(t)],[0 y_pos(t)],'k-'));
fanimator(@(t) text(-0.3,1.5,"Timer: "+num2str(t,2)+" s"));

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