減衰振動のデータ解析

01何を求めるか

揺れる場合( $\gamma < \omega_0$ )の動きは $x = A e^{-\gamma t}\cos(\omega_d t + \phi)$ です(減衰振動の定式化)。データから直接求めるのは、減衰定数 $\gamma$ と揺れの角振動数 $\omega_d$ の 2 つで、残りはそこから計算できます。

抵抗がないときの角振動数	$\omega_0 = \sqrt{\omega_d^2 + \gamma^2}$
Q 値	$Q = \dfrac{\omega_0}{2\gamma}$
ばね定数 (質量 m が要る)	$k = m\,\omega_0^2$
抵抗の係数 (質量 m が要る)	$b = 2m\,\gamma$

02方法 1:山の時刻と高さを使う

まず手で計算できる方法から。山は、おもりが同じ側にいちばん大きく振れた点( $x$ の極大)です。続けて来る山の時刻 $t_n$ と高さ $x_n$ (つり合いの位置から測る)を、記録から書き出します。隣の山の間隔を $T_d$ とすると、隣の山の比は近似なしで $e^{-\gamma T_d}$ です。比の対数 $\delta = \ln(x_n/x_{n+1}) = \gamma T_d$ を対数減衰率と呼びます(対数減衰率と Q 値)。

減衰振動の x(t) と、山の点 x0, x1, x2, x3、隣の山の間隔 T_d — 山の高さ *x_n* と、山の間隔 *T_d*

例題 1質量 $m = 0.200\ \mathrm{kg}$ のおもりで、山が次のようになった。 $\gamma, \omega_0, Q, k, b$ を求める。

$n$	0	1	2	3	4	5
時刻 [s]	0.00	0.80	1.60	2.40	3.21	4.01
高さ [cm]	5.00	4.09	3.35	2.74	2.24	1.84

$T_d$	$\dfrac{t_5 - t_0}{5} = 0.802\ \mathrm{s}$
$\delta$	$\dfrac{1}{5}\ln\dfrac{x_0}{x_5} = 0.200$
$\gamma$	$\delta / T_d = 0.249\ \mathrm{s^{-1}}$
$\omega_d$	$2\pi / T_d = 7.834\ \mathrm{s^{-1}}$

$\omega_0$	$\sqrt{7.834^2 + 0.249^2} = 7.838\ \mathrm{s^{-1}}$
$Q$	$7.838 / (2 \times 0.249) = 15.7$
$k$	$0.200 \times 7.838^2 = 12.3\ \mathrm{N/m}$
$b$	$2 \times 0.200 \times 0.249 = 0.0997\ \mathrm{kg/s}$

$\omega_0$ と $\omega_d$ の差は 0.05% です。 $Q$ が 10 を超えるくらいなら、周期は抵抗がないときとほぼ同じです。

山がたくさんあるなら、 $n$ に対する $t_n$ と $\ln x_n$ に直線を合わせると(最小二乗法)、傾きが $T_d$ と $-\delta$ です。 $\ln x_n$ の誤差は山が小さいほど大きいので、きちんと合わせるなら点ごとの誤差で割ります(最小二乗法)。大まかな値なら、そのままでも十分です。 $\ln x_n$ が直線に乗らないときは、まずつり合いの位置のずれを疑います。それでも曲がるなら抵抗が速さに比例していません(摩擦なら $x_n$ 、空気抵抗の 2 乗の項なら $1/x_n$ がおよそ $n$ の 1 次式)。

演習 1

質量 $m = 0.50\ \mathrm{kg}$ のおもりで、山の時刻と高さが次のようになった。 $\gamma,\ \omega_0,\ Q,\ k,\ b$ を求めよう。

n	0	1	2	3
時刻 [s]	0.00	1.26	2.52	3.78
高さ [cm]	8.00	7.20	6.48	5.83

答えを見る

$T_d = \frac{3.78}{3} = 1.26\ \mathrm{s},\qquad \delta = \frac{1}{3}\ln\frac{8.00}{5.83} = 0.1055$ $\begin{gathered}T_d = \frac{3.78}{3} = 1.26\ \mathrm{s},\\[4pt]\delta = \frac{1}{3}\ln\frac{8.00}{5.83} = 0.1055\end{gathered}$

$\gamma = \frac{0.1055}{1.26} = 0.0837\ \mathrm{s^{-1}},\qquad \omega_d = \frac{2\pi}{1.26} = 4.987\ \mathrm{s^{-1}}$ $\begin{gathered}\gamma = \frac{0.1055}{1.26} = 0.0837\ \mathrm{s^{-1}},\\[4pt]\omega_d = \frac{2\pi}{1.26} = 4.987\ \mathrm{s^{-1}}\end{gathered}$

$\omega_0 = \sqrt{\omega_d^2 + \gamma^2} = 4.987\ \mathrm{s^{-1}}$

$Q = \frac{4.987}{2 \times 0.0837} = 29.8$

$k = 0.50 \times 4.987^2 = 12.4\ \mathrm{N/m},\qquad b = 2 \times 0.50 \times 0.0837 = 0.084\ \mathrm{kg/s}$ $\begin{gathered}k = 0.50 \times 4.987^2 = 12.4\ \mathrm{N/m},\\[4pt]b = 2 \times 0.50 \times 0.0837 = 0.084\ \mathrm{kg/s}\end{gathered}$

$\omega_0$ と $\omega_d$ の差は 0.01% ほどで、この有効数字では見分けがつかない。

03方法 2:式を合わせる

方法 1 は山の数点しか使わず、山がはっきりしないときや、つり合いの位置がずれているときに弱くなります。そこで、記録した点すべてに次の式を最小二乗法で合わせます。 $c$ はつり合いの位置のずれです。

x(t) = A\,e^{-\gamma t}\cos(\omega_d t + \phi) + c

\omega_d = \sqrt{\omega_0^2 - \gamma^2}

反復の初期値には、方法 1 で大まかに求めた $\gamma, \omega_d$ を使います(04 節のコードは FFT で自動に決めます)。例(0.05 s ごとに 6 s 分、測定の誤差 0.08 cm):

減衰振動を記録した点と、それに合わせた式の曲線 — 紺の点 = 記録した点、薄い線 = 合わせた式

	求めた値	データを作った値
$\gamma\ [\mathrm{s^{-1}}]$	$0.2512 \pm 0.0025$	$0.2500$
$\omega_d\ [\mathrm{s^{-1}}]$	$7.8476 \pm 0.0026$	$7.8460$
$A\ [\mathrm{cm}]$	$4.994 \pm 0.028$	$5.000$
$c\ [\mathrm{cm}]$	$0.289 \pm 0.007$	$0.300$

「データを作った値」は、この例の点を作るのに使った値です(実際の実験では分からない)。± は測定の誤差 0.08 cm から出した、求めた値の誤差です(近似直線の誤差)。表の 4 つはどれも、誤差の 2 倍以内で元の値と合っています。

04コードと、AI への頼み方

方法 2 の計算はコンピュータでします。そのまま動かせる Python のコードです。data.csv に記録を入れ、sigma を自分の測定の誤差に変えて使います。

Python

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# data.csv: 1 行目は見出し、1 列目が時刻 [s]、2 列目が位置 [cm]
t, x = np.loadtxt("data.csv", delimiter=",", skiprows=1, unpack=True)
sigma = 0.08  # 1 点の測定の誤差 [cm]。自分の値に変える

def model(t, A, g, wd, phi, c):
    return A * np.exp(-g * t) * np.cos(wd * t + phi) + c

# 初期値: 角振動数は FFT のいちばん強い成分から。位相は 4 通り試して、残差がいちばん小さいものを使う
f = np.fft.rfftfreq(t.size, t[1] - t[0])
spec = np.abs(np.fft.rfft(x - x.mean()))
wd0 = 2 * np.pi * f[spec[1:].argmax() + 1]
A0 = (x.max() - x.min()) / 2
bounds = ([0, 0, 0, -np.inf, -np.inf], [np.inf, np.inf, np.inf, np.inf, np.inf])   # A > 0, γ ≥ 0, ω_d > 0
best = None
for phi0 in [0, np.pi / 2, np.pi, -np.pi / 2]:
    try:
        p_, cov_ = curve_fit(model, t, x, p0=[A0, 1 / (t[-1] - t[0]), wd0, phi0, x.mean()],
                             sigma=np.full(t.size, sigma), absolute_sigma=True, bounds=bounds)
    except RuntimeError:
        continue
    s_ = np.sum((x - model(t, *p_)) ** 2)
    if best is None or s_ < best[0]:
        best = (s_, p_, cov_)
_, p, cov = best
p[3] = (p[3] + np.pi) % (2 * np.pi) - np.pi   # φ を -π から π に
err = np.sqrt(np.diag(cov))
A, g, wd, phi, c = p
w0 = np.sqrt(wd**2 + g**2)
print(f"gamma = {g:.4f} +- {err[1]:.4f} 1/s")
print(f"omega_d = {wd:.4f} +- {err[2]:.4f} 1/s")
# ω0 = √(ω_d² + γ²)、Q = ω0/(2γ) の誤差を、共分散行列で伝える (誤差伝播)
Q = w0 / (2 * g)
J0 = np.array([0, g / w0, wd / w0, 0, 0])                        # ∂ω0/∂(A, γ, ω_d, φ, c)
JQ = np.array([0, (g / w0) / (2 * g) - w0 / (2 * g**2), (wd / w0) / (2 * g), 0, 0])
w0_err, Q_err = np.sqrt(J0 @ cov @ J0), np.sqrt(JQ @ cov @ JQ)
print(f"omega_0 = {w0:.4f} +- {w0_err:.4f} 1/s")
print(f"Q = {Q:.2f} +- {Q_err:.2f}")
print(f"c (zero offset) = {c:.3f} +- {err[4]:.3f} cm")

fig, (a1, a2) = plt.subplots(2, 1, sharex=True)
a1.plot(t, x, ".", ms=3); a1.plot(t, model(t, *p)); a1.set_ylabel("x [cm]")
a2.plot(t, x - model(t, *p), ".", ms=3); a2.set_ylabel("residual [cm]"); a2.set_xlabel("t [s]")
plt.show()

AI にコードを書いてもらうなら、何をしたいかを 1 つずつ書いて渡します。

AI に渡す文

減衰振動の実験データ data.csv があります。1 行目は見出し、1 列目が時刻 [s] (等間隔)、2 列目が位置 [cm] です。
Python (numpy, scipy, matplotlib) で、次のことをするコードを書いてください。

1. x(t) = A e^(-γt) cos(ω_d t + φ) + c を、scipy.optimize.curve_fit で最小二乗法により合わせる。c はつり合いの位置のずれ。A > 0、γ > 0、ω_d > 0 に制限し、φ は -π から π の範囲で表示する。
2. 初期値はデータから自動で決める。角振動数は FFT のいちばん強い成分、振幅は最大と最小の差の半分、c は平均、γ は 1/(測定時間)。φ は 0, π/2, π, 3π/2 の 4 通りから始めて、うまくいった中で残差の 2 乗和がいちばん小さい結果を使う。
3. 1 点の測定の誤差を 0.08 cm として、absolute_sigma=True でパラメータの誤差も出す。
4. γ、ω_d、c と、ω_0 = √(ω_d² + γ²)、Q = ω_0 / (2γ) を誤差つきで表示する。ω_0 と Q の誤差は共分散行列を使って誤差伝播で出す。
5. データと合わせた曲線のグラフと、残差 (データ − 式) のグラフを上下に並べて描く。

出てきた値は、方法 1 で出した値とおよそ合うか、残差のグラフがでたらめに散らばっているかを必ず確かめます(残差を見る)。

05どちらを使うか

山の時刻と高さを使う:手で計算できる。山の時刻と高さの誤差と、つり合いの位置のずれが、そのまま結果の誤差になる。
式を合わせる:コンピュータが要る。全部の点を使い、つり合いの位置のずれ $c$ も一緒に求まり、値の誤差も出る。