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

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

Plot showing the amplitude of a damped harmonic oscillator over time, for underdamped, overdamped, and critically damped conditions

内容

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

1.運動方程式の導出

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

Mass-spring-damper system

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

$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("\omega_0 t");
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 omega indexOf 0 baseline t, 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 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 = 1, xlabel omega indexOf 0 baseline 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 t, ylabel amplitude contains 6 objects of type functionline, text. These objects represent \pi, 1, 1/\pi.$

参考

関数

diff | subs | dsolve | fplot | collect | expand | rewrite