Main Content

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

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

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

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

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

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

a=rd2θdt2.

subsa の代わりに使用します。

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

mr2t2 θ(t)=-gmsin(θ(t))

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

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

2t2 θ(t)=-gsin(θ(t))r

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

ω0= gr

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

2t2 θ(t)=-ω02sin(θ(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 = 

2t2 θ(t)=-ω02θ(t)

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

方程式 eqnLineardsolve を使用して解きます。初期条件を 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) = 

θ0cos(ω0t)+θt0sin(ω0t)ω0

手順 4: ω0 の物理的意味

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

T=2πω0=2πrg.

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

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

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

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

  • 重力加速度 g=9.81m/s2

  • 振子の長さ r=1m

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 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=12mr2(dθdt)2+mgr(1-cosθ)

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

Emr2=12[(dθdt)2+(2ω0sinθ2)2]

エネルギーは保存されるため、振子の運動は位相空間内の定エネルギー軌道によって記述できます。位相空間は座標 θ および dθ/dt をもつ抽象空間です。これらの軌道を 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 C o n s t a n t blank E n e r g y blank C o n t o u r s blank i n blank P h a s e blank S p a c e blank ( blank theta blank v s . blank theta indexOf t baseline blank ) contains an object of type functioncontour.

一定エネルギーの等高線は θ 軸および dθ/dt 軸に対して対称であり、θ 軸に沿って周期的です。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) = 

(t θ(t)=θt(t)t θt(t)=-ω02sin(θ(t)))

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

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

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

(1001)

F = 

(θt(t)-981sin(θ(t))100)

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

M(t,x(t))dxdt=F(t,x(t)).

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

f = M\F
f = 

(θt(t)-981sin(θ(t))100)

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θt0/2ω01 を満たします。θ および dθ/dt の初期条件を変数 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 ... ]
          y: [2x45 double]
      stats: [1x1 struct]
      idata: [1x1 struct]

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

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

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 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.

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

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

ode45 を使用して、開いたエネルギー等高線の ODE を解きます。エネルギー等高線図から、開いた等高線は条件 θ0=0θt0/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 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.

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