減衰調和振動子の物理特性

この例では、駆動力がない場合の運動方程式を解くことにより、減衰調和振動子の物理特性を調べます。この例では、不足減衰、過減衰、および臨界減衰の各ケースについて調査します。

内容

運動方程式の導出
運動方程式 (F = 0) の求解
不足減衰のケース ( $ζ < 1$ )
過減衰のケース ( $ζ > 1$ )
臨界減衰のケース ( $ζ = 1$ )
まとめ

1.運動方程式の導出

以下に示す減衰での強制調和振動子を考えます。振動子が動く速度に比例した抵抗力をモデル化します。

以下の場合の運動方程式を定義します。

$m$ が質量
$c$ が減衰係数
$k$ がばね定数
$F$ が駆動力

syms x(t) m c k F(t)
eq = m*diff(x,t,t) + c*diff(x,t) + k*x == F

eq(t) = 
 $m \frac{\partial^{2}}{\partial t^{2}} x (t) + c \frac{\partial}{\partial t} x (t) + k x (t) = F (t)$

$c = m γ$ および $k = m ω_{0}^{2}$ を使用して方程式を書き換えます。

syms gamma omega_0
eq = subs(eq, [c k], [m*gamma, m*omega_0^2])

eq(t) = 
 $m \frac{\partial^{2}}{\partial t^{2}} x (t) + γ m \frac{\partial}{\partial t} x (t) + m x (t) {ω_{0}}^{2} = F (t)$

質量 $m$ を除算します。これで、分析に使いやすい形式の方程式が得られました。

eq = collect(eq, m)/m

eq(t) = 
 $\frac{\partial^{2}}{\partial t^{2}} x (t) + γ \frac{\partial}{\partial t} x (t) + x (t) {ω_{0}}^{2} = \frac{F (t)}{m}$

2.運動方程式 (F = 0) の求解

$F = 0$ であるような、外力がない場合に、dsolve を使用して運動方程式を解きます。単位変位および速度ゼロの初期条件を使用します。

vel = diff(x,t);
cond = [x(0) == 1, vel(0) == 0];
eq = subs(eq,F,0);
sol = dsolve(eq, cond)

sol = 
 $\begin{array}{l} \frac{e^{- t (\frac{γ}{2} - \frac{σ_{1}}{2})} (γ + σ_{1})}{2 σ_{1}} - \frac{e^{- t (\frac{γ}{2} + \frac{σ_{1}}{2})} (γ - σ_{1})}{2 σ_{1}} \\ where \\ σ_{1} = \sqrt{(γ - 2 ω_{0}) (γ + 2 ω_{0})} \end{array}$

解を展開して単純化する方法を確認します。

sol = expand(sol)

sol = 
 $\begin{array}{l} \frac{σ_{1} σ_{2}}{2} + \frac{σ_{1} e^{\frac{t \sqrt{γ^{2} - 4 {ω_{0}}^{2}}}{2}}}{2} - \frac{γ σ_{1} σ_{2}}{2 \sqrt{γ^{2} - 4 {ω_{0}}^{2}}} + \frac{γ σ_{1} e^{\frac{t \sqrt{γ^{2} - 4 {ω_{0}}^{2}}}{2}}}{2 \sqrt{γ^{2} - 4 {ω_{0}}^{2}}} \\ where \\ σ_{1} = e^{- \frac{γ t}{2}} \\ σ_{2} = e^{- \frac{t \sqrt{γ^{2} - 4 {ω_{0}}^{2}}}{2}} \end{array}$

各項には、 $σ_{1}$ 、つまり $e^{- γ t / 2}$ の因数があることに注意してください。collect を使用してこれらの項を集約します。

sol = collect(sol, exp(-gamma*t/2))

sol = 
 $\begin{array}{l} (\frac{e^{- σ_{1}}}{2} + \frac{e^{σ_{1}}}{2} - \frac{γ e^{- σ_{1}}}{2 \sqrt{γ^{2} - 4 {ω_{0}}^{2}}} + \frac{γ e^{σ_{1}}}{2 \sqrt{γ^{2} - 4 {ω_{0}}^{2}}}) e^{- \frac{γ t}{2}} \\ where \\ σ_{1} = \frac{t \sqrt{γ^{2} - 4 {ω_{0}}^{2}}}{2} \end{array}$

項 $\sqrt{γ^{2} - 4 ω_{0}^{2}}$ は、解のさまざまな部分に現れます。減衰比 $ζ \equiv \frac{γ}{2 ω_{0}}$ を導入することで、より単純な形式に書き換えます。

ζ を上記の項に代入すると、次の結果を得ます。

$\sqrt{γ^{2} - 4 ω_{0}^{2}} = 2 ω_{0} \sqrt{{(\frac{γ}{2 ω_{0}})}^{2} - 1} = 2 ω_{0} \sqrt{ζ^{2} - 1},$

syms zeta;
sol = subs(sol, ...
    sqrt(gamma^2 - 4*omega_0^2), ...
    2*omega_0*sqrt(zeta^2-1))

sol = 
 $\begin{array}{l} e^{- \frac{γ t}{2}} (\frac{σ_{2}}{2} + \frac{σ_{1}}{2} + \frac{γ σ_{2}}{4 ω_{0} \sqrt{ζ^{2} - 1}} - \frac{γ σ_{1}}{4 ω_{0} \sqrt{ζ^{2} - 1}}) \\ where \\ σ_{1} = e^{- ω_{0} t \sqrt{ζ^{2} - 1}} \\ σ_{2} = e^{ω_{0} t \sqrt{ζ^{2} - 1}} \end{array}$

$ω_{0}$ および $ζ$ で $γ$ を置き換えることで、解をさらに単純化します。

sol = subs(sol, gamma, 2*zeta*omega_0)

sol = 
 $\begin{array}{l} e^{- ω_{0} t ζ} (\frac{σ_{2}}{2} + \frac{σ_{1}}{2} + \frac{ζ σ_{2}}{2 \sqrt{ζ^{2} - 1}} - \frac{ζ σ_{1}}{2 \sqrt{ζ^{2} - 1}}) \\ where \\ σ_{1} = e^{- ω_{0} t \sqrt{ζ^{2} - 1}} \\ σ_{2} = e^{ω_{0} t \sqrt{ζ^{2} - 1}} \end{array}$

駆動力のない減衰調和振動子の運動の一般解を得ました。次に、運動がより単純な形をとる、減衰比 $ζ$ の 3 つの特殊なケースについて検討します。これらのケースは以下のように呼ばれます。

不足減衰 $(ζ < 1)$ 、
過減衰 $(ζ > 1)$ 、および
臨界減衰 $(ζ = 1)$ 。

3.不足減衰のケース ( $ζ < 1$ )

$ζ < 1$ の場合、 $\sqrt{ζ^{2} - 1} = i \sqrt{1 - ζ^{2}}$ は純虚数です。

solUnder = subs(sol, sqrt(zeta^2-1), 1i*sqrt(1-zeta^2))

solUnder = 
 $\begin{array}{l} e^{- ω_{0} t ζ} (\frac{σ_{1}}{2} + \frac{σ_{2}}{2} + \frac{ζ σ_{1} i}{2 \sqrt{1 - ζ^{2}}} - \frac{ζ σ_{2} i}{2 \sqrt{1 - ζ^{2}}}) \\ where \\ σ_{1} = e^{- ω_{0} t \sqrt{1 - ζ^{2}} i} \\ σ_{2} = e^{ω_{0} t \sqrt{1 - ζ^{2}} i} \end{array}$

上記の方程式内の項 $e^{i ω_{0} t \sqrt{ζ^{2} - 1}} \pm e^{- {i ω}_{0} t \sqrt{ζ^{2} - 1}}$ に注目すると、 $e^{i x} = \cos (x) + i \sin (x) .$ という恒等式が想起されます。

解を $\cos$ で書き換えます。

solUnder = coeffs(solUnder, zeta);
solUnder = solUnder(1);
c = exp(-omega_0 * zeta * t);
solUnder = c * rewrite(solUnder / c, 'cos')

solUnder =  $e^{- ω_{0} t ζ} \cos (ω_{0} t \sqrt{1 - ζ^{2}})$

solUnder(t, omega_0, zeta) = solUnder

solUnder(t, omega_0, zeta) =  $e^{- ω_{0} t ζ} \cos (ω_{0} t \sqrt{1 - ζ^{2}})$

系は、 $ω_{0} \sqrt{1 - ζ^{2}}$ の固有振動数で振動し、 $1 / ω_{0} ζ$ の指数関数的速度で減衰します。

$ω_{0} t$ と $ζ$ の関数として解を fplot でプロットします。

z = [0 1/4 1/2 3/4];
w = 1;
T = 4*pi;
lineStyle = {'-','--',':k','-.'};

fplot(@(t)solUnder(t, w, z(1)), [0 T], lineStyle{1});

hold on;
for k = 2:numel(z)
    fplot(@(t)solUnder(t, w, z(k)), [0 T], lineStyle{k});
end
hold off;
grid on;
xticks(T*linspace(0,1,5));
xticklabels({'0','\pi','2\pi','3\pi','4\pi'});
xlabel('t / \omega_0');
ylabel('amplitude');
lgd = legend('0','1/4','1/2','3/4');
title(lgd,'\zeta');
title('Underdamped');

Figure contains an axes object. The axes object with title Underdamped, xlabel t / blank omega indexOf 0 baseline, ylabel amplitude contains 4 objects of type functionline. These objects represent 0, 1/4, 1/2, 3/4.

4.過減衰のケース ( $ζ > 1$ )

$ζ > 1$ の場合、 $\sqrt{ζ^{2} - 1}$ は純実数であり、解は次のように書き換えられます。

solOver = sol

solOver = 
 $\begin{array}{l} e^{- ω_{0} t ζ} (\frac{σ_{2}}{2} + \frac{σ_{1}}{2} + \frac{ζ σ_{2}}{2 \sqrt{ζ^{2} - 1}} - \frac{ζ σ_{1}}{2 \sqrt{ζ^{2} - 1}}) \\ where \\ σ_{1} = e^{- ω_{0} t \sqrt{ζ^{2} - 1}} \\ σ_{2} = e^{ω_{0} t \sqrt{ζ^{2} - 1}} \end{array}$

solOver = coeffs(solOver, zeta);
solOver = solOver(1)

solOver = 
 $e^{- ω_{0} t ζ} (\frac{e^{ω_{0} t \sqrt{ζ^{2} - 1}}}{2} + \frac{e^{- ω_{0} t \sqrt{ζ^{2} - 1}}}{2})$

項 $\frac{(e^{ω_{0} t \sqrt{ζ^{2} - 1}} + e^{- ω_{0} t \sqrt{ζ^{2} - 1}})}{2}$ に注目すると、 $\cosh (x) = \frac{e^{x} + e^{- x}}{2}$ という恒等式が想起されます。

式を $\cosh$ で書き換えます。

c = exp(-omega_0*t*zeta);
solOver = c*rewrite(solOver / c, 'cosh')

solOver =  $\cosh (ω_{0} t \sqrt{ζ^{2} - 1}) e^{- ω_{0} t ζ}$

solOver(t, omega_0, zeta) = solOver

solOver(t, omega_0, zeta) =  $\cosh (ω_{0} t \sqrt{ζ^{2} - 1}) e^{- ω_{0} t ζ}$

解をプロットし、振動なしで減衰することを確認します。

z = 1 + [1/4 1/2 3/4 1];
w = 1;
T = 4*pi;
lineStyle = {'-','--',':k','-.'};

fplot(@(t)solOver(t, w, z(1)), [0 T], lineStyle{1});

hold on;
for k = 2:numel(z)
    fplot(@(t)solOver(t, w, z(k)), [0 T], lineStyle{k});
end
hold off;
grid on;
xticks(T*linspace(0,1,5));
xticklabels({'0','\pi','2\pi','3\pi','4\pi'});
xlabel('\omega_0 t');
ylabel('amplitude');
lgd = legend('1+1/4','1+1/2','1+3/4','2');
title(lgd,'\zeta');
title('Overdamped');

Figure contains an axes object. The axes object with title Overdamped, xlabel omega indexOf 0 baseline blank t, ylabel amplitude contains 4 objects of type functionline. These objects represent 1+1/4, 1+1/2, 1+3/4, 2.

5.臨界減衰のケース ( $ζ = 1$ )

$ζ = 1$ の場合、解は次のように単純化されます。

solCritical(t, omega_0) = limit(sol, zeta, 1)

solCritical(t, omega_0) =  $e^{- ω_{0} t} (ω_{0} t + 1)$

臨界減衰のケースの解をプロットします。

w = 1;
T = 4*pi;

fplot(solCritical(t, w), [0 T])
xlabel('\omega_0 t');
ylabel('x');
title('Critically damped, \zeta = 1');
grid on;
xticks(T*linspace(0,1,5));
xticklabels({'0','\pi','2\pi','3\pi','4\pi'});

Figure contains an axes object. The axes object with title Critically damped, zeta blank = blank 1, xlabel omega indexOf 0 baseline blank t, ylabel x contains an object of type functionline.

6.まとめ

減衰比 $ζ$ を用いてその運動を表す ODE を解くことで、調和振動子の異なる減衰状態を調べました。3 つすべてのケースを一緒にプロットし、それらを比較して対比します。

zOver  = pi;
zUnder = 1/zOver;
w = 1;
T = 2*pi;
lineStyle = {'-','--',':k'};

fplot(@(t)solOver(t, w, zOver), [0 T], lineStyle{1},'LineWidth',2);
hold on;
fplot(solCritical(t, w), [0 T], lineStyle{2},'LineWidth',2)
fplot(@(t)solUnder(t, w, zUnder), [0 T], lineStyle{3},'LineWidth',2);
hold off;

textColor = lines(3);
text(3*pi/2, 0.3 , 'over-damped'      ,'Color',textColor(1,:));
text(pi*3/4, 0.05, 'critically-damped','Color',textColor(2,:));
text(pi/8  , -0.1, 'under-damped');

grid on;
xlabel('\omega_0 t');
ylabel('amplitude');
xticks(T*linspace(0,1,5));
xticklabels({'0','\pi/2','\pi','3\pi/2','2\pi'});
yticks((1/exp(1))*[-1 0 1 2 exp(1)]);
yticklabels({'-1/e','0','1/e','2/e','1'});
lgd = legend('\pi','1','1/\pi');
title(lgd,'\zeta');
title('Damped Harmonic Oscillator');

$Figure contains an axes object. The axes object with title Damped Harmonic Oscillator, xlabel omega indexOf 0 baseline blank t, ylabel amplitude contains 6 objects of type functionline, text. These objects represent \pi, 1, 1/\pi.$