Main Content

このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。

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

この例では減衰調和振動子の物理特性について、以下を行うことにより考察します。

  • 駆動力がない場合の運動方程式の求解

  • 不足減衰、過減衰、および臨界減衰の各ケースについての調査

内容

  1. 運動方程式の導出

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

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

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

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

  6. まとめ

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) = 

m2t2 x(t)+ct x(t)+kx(t)=F(t)m*diff(x(t), t, 2) + c*diff(x(t), t) + k*x(t) == F(t)

c=mγ および k=mω02 を使用して方程式を書き換えます。

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

m2t2 x(t)+γmt x(t)+mx(t)ω02=F(t)m*diff(x(t), t, 2) + gamma*m*diff(x(t), t) + m*x(t)*omega_0^2 == F(t)

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

eq = collect(eq, m)/m
eq(t) = 

2t2 x(t)+γt x(t)+x(t)ω02=F(t)mdiff(x(t), t, 2) + gamma*diff(x(t), t) + x(t)*omega_0^2 == 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 = 

e-tγ2-σ12γ+σ12σ1-e-tγ2+σ12γ-σ12σ1where  σ1=γ-2ω0γ+2ω0(exp((-t*(gamma/2 - sqrt((gamma - 2*omega_0)*(gamma + 2*omega_0))/2)))*(gamma + sqrt((gamma - 2*omega_0)*(gamma + 2*omega_0))))/(2*sqrt((gamma - 2*omega_0)*(gamma + 2*omega_0))) - (exp((-t*(gamma/2 + sqrt((gamma - 2*omega_0)*(gamma + 2*omega_0))/2)))*(gamma - sqrt((gamma - 2*omega_0)*(gamma + 2*omega_0))))/(2*sqrt((gamma - 2*omega_0)*(gamma + 2*omega_0)))

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

sol = expand(sol)
sol = 

σ1σ22+σ1etγ2-4ω0222-γσ1σ22γ2-4ω02+γσ1etγ2-4ω0222γ2-4ω02where  σ1=e-γt2  σ2=e-tγ2-4ω022(exp(-(gamma*t)/2)*exp(-(t*sqrt(gamma^sym(2) - 4*omega_0^2))/2))/2 + (exp(-(gamma*t)/2)*exp(((t*sqrt(gamma^sym(2) - 4*omega_0^2))/2)))/2 - (gamma*exp(-(gamma*t)/2)*exp(-(t*sqrt(gamma^sym(2) - 4*omega_0^2))/2))/(2*sqrt(gamma^sym(2) - 4*omega_0^2)) + (gamma*exp(-(gamma*t)/2)*exp(((t*sqrt(gamma^sym(2) - 4*omega_0^2))/2)))/(2*sqrt(gamma^sym(2) - 4*omega_0^2))

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

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

e-σ12+eσ12-γe-σ12γ2-4ω02+γeσ12γ2-4ω02e-γt2where  σ1=tγ2-4ω022(exp((-(t*sqrt(gamma^sym(2) - 4*omega_0^2))/2))/2 + exp((t*sqrt(gamma^sym(2) - 4*omega_0^2))/2)/2 - (gamma*exp((-(t*sqrt(gamma^sym(2) - 4*omega_0^2))/2)))/(2*sqrt(gamma^sym(2) - 4*omega_0^2)) + (gamma*exp(((t*sqrt(gamma^sym(2) - 4*omega_0^2))/2)))/(2*sqrt(gamma^sym(2) - 4*omega_0^2)))*exp((-(gamma*t)/2))

γ2-4ω02 は、解のさまざまな部分に現れます。減衰比 ζγ2ω0 を導入することで、より単純な形式に書き換えます。

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

γ24ω02=2ω0(γ2ω0)21=2ω0ζ21,

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

e-γt2σ22+σ12+γσ24ω0ζ2-1-γσ14ω0ζ2-1where  σ1=e-ω0tζ2-1  σ2=eω0tζ2-1exp((-(gamma*t)/2))*(exp(omega_0*t*sqrt(zeta^sym(2) - 1))/2 + exp(-omega_0*t*sqrt(zeta^sym(2) - 1))/2 + (gamma*exp(omega_0*t*sqrt(zeta^sym(2) - 1)))/(4*omega_0*sqrt(zeta^sym(2) - 1)) - (gamma*exp(-omega_0*t*sqrt(zeta^sym(2) - 1)))/(4*omega_0*sqrt(zeta^sym(2) - 1)))

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

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

e-ω0tζσ22+σ12+ζσ22ζ2-1-ζσ12ζ2-1where  σ1=e-ω0tζ2-1  σ2=eω0tζ2-1exp((-omega_0*t*zeta))*(exp(omega_0*t*sqrt(zeta^sym(2) - 1))/2 + exp(-omega_0*t*sqrt(zeta^sym(2) - 1))/2 + (zeta*exp(omega_0*t*sqrt(zeta^sym(2) - 1)))/(2*sqrt(zeta^sym(2) - 1)) - (zeta*exp(-omega_0*t*sqrt(zeta^sym(2) - 1)))/(2*sqrt(zeta^sym(2) - 1)))

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

  • 不足減衰 (ζ<1)

  • 過減衰 (ζ>1)、および

  • 臨界減衰 (ζ=1)

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

ζ<1 の場合、ζ2-1=i1-ζ2 は純虚数です。

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

e-ω0tζσ12+σ22+ζσ1i21-ζ2-ζσ2i21-ζ2where  σ1=e-ω0t1-ζ2i  σ2=eω0t1-ζ2iexp((-omega_0*t*zeta))*(exp(-omega_0*t*sqrt(1 - zeta^sym(2))*sym(1i))/2 + exp(omega_0*t*sqrt(1 - zeta^sym(2))*sym(1i))/2 + (zeta*exp(-omega_0*t*sqrt(1 - zeta^sym(2))*sym(1i))*sym(1i))/(2*sqrt(1 - zeta^sym(2))) - (zeta*exp(omega_0*t*sqrt(1 - zeta^sym(2))*sym(1i))*sym(1i))/(2*sqrt(1 - zeta^sym(2))))

上記の方程式内の項 eiω0tζ2-1±e-iω0tζ2-1 に注目すると、eix=cos(x)+isin(x). という恒等式が想起されます。

解を cos で書き換えます。

solUnder = coeffs(solUnder, zeta);
solUnder = solUnder(1);
c = exp(-omega_0 * zeta * t);
solUnder = c * rewrite(solUnder / c, 'cos')
solUnder = e-ω0tζcos(ω0t1-ζ2)exp((-omega_0*t*zeta))*cos(omega_0*t*sqrt(1 - zeta^sym(2)))
solUnder(t, omega_0, zeta) = solUnder
solUnder(t, omega_0, zeta) = e-ω0tζcos(ω0t1-ζ2)exp((-omega_0*t*zeta))*cos(omega_0*t*sqrt(1 - zeta^sym(2)))

系は、ω01-ζ2 の固有振動数で振動し、1/ω0ζ の指数関数的速度で減衰します。

ω0tζ の関数として解を 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. The axes with title Underdamped contains 4 objects of type functionline. These objects represent 0, 1/4, 1/2, 3/4.

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

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

solOver = sol
solOver = 

e-ω0tζσ22+σ12+ζσ22ζ2-1-ζσ12ζ2-1where  σ1=e-ω0tζ2-1  σ2=eω0tζ2-1exp((-omega_0*t*zeta))*(exp(omega_0*t*sqrt(zeta^sym(2) - 1))/2 + exp(-omega_0*t*sqrt(zeta^sym(2) - 1))/2 + (zeta*exp(omega_0*t*sqrt(zeta^sym(2) - 1)))/(2*sqrt(zeta^sym(2) - 1)) - (zeta*exp(-omega_0*t*sqrt(zeta^sym(2) - 1)))/(2*sqrt(zeta^sym(2) - 1)))

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

e-ω0tζeω0tζ2-12+e-ω0tζ2-12exp((-omega_0*t*zeta))*(exp(omega_0*t*sqrt(zeta^sym(2) - 1))/2 + exp((-omega_0*t*sqrt(zeta^sym(2) - 1)))/2)

(eω0tζ2-1+e-ω0tζ2-1)2 に注目すると、cosh(x)=ex+e-x2 という恒等式が想起されます。

式を cosh で書き換えます。

c = exp(-omega_0*t*zeta);
solOver = c*rewrite(solOver / c, 'cosh')
solOver = cosh(ω0tζ2-1)e-ω0tζcosh(omega_0*t*sqrt(zeta^sym(2) - 1))*exp((-omega_0*t*zeta))
solOver(t, omega_0, zeta) = solOver
solOver(t, omega_0, zeta) = cosh(ω0tζ2-1)e-ω0tζcosh(omega_0*t*sqrt(zeta^sym(2) - 1))*exp((-omega_0*t*zeta))

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

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. The axes with title Overdamped 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-ω0tω0t+1exp((-omega_0*t))*(omega_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. The axes with title Critically damped, \zeta = 1 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. The axes with title Damped Harmonic Oscillator contains 6 objects of type functionline, text. These objects represent \pi, 1, 1/\pi.