「え、バネの硬さって、伸ばし方によって変わるの？」
そうなんです。学校で習うフックの法則「$F = kx$」は、あくまで理想的な線形バネの話。現実の世界では、ゴムが伸びきると急に硬くなったり、ある種の金属は大きな変形で柔らかくなったりします。この「バネの硬さが変わる振動」を非線形振動と呼び、自動車のエンジンマウント、超高層ビルの免震装置、さらにはスマホの中のMEMSセンサーまで、あらゆる工学分野で登場します。
そしてこの非線形振動、「周波数をゆっくり変えるだけで振幅が突然ジャンプする」 という、線形の常識を覆す現象が起こります。今回は、その謎を数式→実装→シミュレーションで完全に理解しましょう。
ざっくり本質：バネが「気まぐれ」になると何が起きる？
「非線形振動とは、振幅が大きくなると固有振動数が変化する振動である」
線形バネでは、振幅がどんなに大きくても固有振動数は一定。だから共振はいつも同じ周波数で起こります。
ところが非線形バネでは：
- ハードニング（硬化）型：振幅が大きいほどバネが硬くなる → 共振周波数が高くなる
- ソフトニング（軟化）型：振幅が大きいほどバネが柔らかくなる → 共振周波数が低くなる
この「振幅に応じて共振点が動く」という性質が、ジャンプ現象や複数解の存在など、線形では絶対に起こらない面白い挙動を生み出します。
数式で理解する：Duffing方程式の核心
非線形振動の代表格がDuffing振動子。その運動方程式は：
\ddot{x} + 2\zeta\omega_0\dot{x} + \omega_0^2 x + \varepsilon x^3 = F\cos(\Omega t)
各項の意味：
- $\ddot{x}$：加速度（慣性力）
- $2\zeta\omega_0\dot{x}$：減衰項（$\zeta$は減衰比）
- $\omega_0^2 x$：線形剛性（$\omega_0$は線形固有角周波数）
- $\varepsilon x^3$：非線形項（これが全ての鍵！$\varepsilon > 0$でハードニング、$\varepsilon < 0$でソフトニング）
- $F\cos(\Omega t)$：外部からの強制振動
この非線形項があるせいで、解析的に厳密解を求めるのは不可能。そこで登場するのが調和バランス法です。
調和バランス法の一次近似
「応答は励振と同じ周波数$\Omega$の単一余弦波で近似できる」と仮定します：
x(t) \approx A\cos(\Omega t)
これを運動方程式に代入し、$\cos(\Omega t)$の係数を比較すると、以下の振幅-周波数関係式が得られます：
\left[(\omega_0^2 + \tfrac{3}{4}\varepsilon A^2 - \Omega^2)^2 + (2\zeta\omega_0\Omega)^2\right]A^2 = F^2
この式が意味すること：
- 左辺の $(\omega_0^2 + \tfrac{3}{4}\varepsilon A^2 - \Omega^2)^2$ の部分で、非線形項が振幅$A$に依存した固有振動数を作り出している
- $\omega_0^2$ が $\omega_0^2 + \tfrac{3}{4}\varepsilon A^2$ に置き換わったと解釈できる
- この式を$A$について解くと、ある周波数$\Omega$に対して複数の振幅解が現れる（これがジャンプ現象の正体！）
骨格曲線（スケルトンカーブ）
減衰と励振がない場合（$\zeta=0, F=0$）の自由振動における振幅と周波数の関係：
\Omega_{bb}= \omega_0\sqrt{1 + \frac{3\varepsilon A^2}{4\omega_0^2}}
これが骨格曲線。振幅が大きくなるにつれて、ハードニングでは高周波側、ソフトニングでは低周波側に曲がっていく様子を示します。
コードで実装する：振幅-周波数応答を計算する
ここが一番大事。上の代数方程式を数値的に解いて、周波数応答曲線を描くPythonコードです。
import numpy as np
def duffing_frequency_response(omega0=1.0, zeta=0.02, epsilon=0.5, F=0.2, omega_min=0.5, omega_max=2.0, n_points=500):
"""
Duffing振動子の周波数応答（振幅 vs 周波数）を計算する
パラメータ:
omega0: 線形固有角周波数
zeta: 減衰比
epsilon: 非線形係数（正:ハードニング, 負:ソフトニング）
F: 励振力振幅
omega_min, omega_max: 周波数掃引範囲
n_points: 周波数分割数
戻り値:
omega_list: 周波数配列
A_list: 対応する振幅配列（安定な実数解のみ）
"""
# 解くべき方程式: [(ω0^2 + 3/4 ε A^2 - Ω^2)^2 + (2ζω0Ω)^2] A^2 = F^2
# これをAについて解く → 両辺の差が0になるAを探索
omega_list = np.linspace(omega_min, omega_max, n_points)
A_list = []
for Omega in omega_list:
# 振幅Aを0.001から5.0まで細かくスキャン
A_candidates = np.linspace(0.001, 5.0, 2000)
# 方程式の左辺 - 右辺（誤差）
left = ((omega0**2 + 0.75 * epsilon * A_candidates**2 - Omega**2)**2
+ (2 * zeta * omega0 * Omega)**2) * A_candidates**2
error = left - F**2
# 符号が変化する点が解（振幅が存在する周波数）
# 誤差が最小の点を解として採用（複数解がある場合、最も近いものを選ぶ）
sign_changes = np.where(np.diff(np.sign(error)))[0]
if len(sign_changes) > 0:
# 各符号変化点の前後で線形補間して精密な解を求める
for idx in sign_changes:
A1, A2 = A_candidates[idx], A_candidates[idx+1]
e1, e2 = error[idx], error[idx+1]
# 線形補間
A_solution = A1 - e1 * (A2 - A1) / (e2 - e1)
if 0 < A_solution < 5.0:
A_list.append(A_solution)
else:
# 解が見つからない場合はNaN（プロットで無視される）
A_list.append(np.nan)
return omega_list, np.array(A_list)
# 実際に計算してみる
omega0 = 1.0
zeta = 0.02
epsilon = 0.5 # ハードニング
F = 0.2
omega_vals, A_vals = duffing_frequency_response(
omega0=omega0, zeta=zeta, epsilon=epsilon, F=F,
omega_min=0.6, omega_max=1.8, n_points=1000
)
# 結果表示（最初の10点）
print("Ω/ω0, 振幅A")
for i in range(min(10, len(omega_vals))):
if not np.isnan(A_vals[i]):
print(f"{omega_vals[i]:.3f}, {A_vals[i]:.4f}")
# 骨格曲線も計算
A_skel = np.linspace(0.01, 2.0, 100)
omega_skel = omega0 * np.sqrt(1 + 0.75 * epsilon * A_skel**2 / omega0**2)
このコードのポイント：
- 振幅$A$を細かくスキャンして方程式の誤差がゼロになる点を探す
- 複数の解が見つかる場合がある（これがジャンプ現象の原因！）
- 骨格曲線は別途計算
数値例で確かめる：具体的な値を入れてみる
条件設定：
- $\omega_0 = 1.0$ [rad/s]
- $\zeta = 0.02$（2%の減衰）
- $\varepsilon = 0.5$（ハードニング）
- $F = 0.2$ [N/kg]
周波数$\Omega = 1.2$ [rad/s] における振幅を計算：
方程式：
[(1.0^2 + \tfrac{3}{4} \times 0.5 \times A^2 - 1.2^2)^2 + (2 \times 0.02 \times 1.0 \times 1.2)^2] A^2 = 0.2^2
計算過程：
- まず線形項：$1.0^2 = 1.0$、$\Omega^2 = 1.44$
- 減衰項：$(2 \times 0.02 \times 1.0 \times 1.2)^2 = (0.048)^2 = 0.002304$
- 非線形項を含む：$(1.0 + 0.375A^2 - 1.44)^2 = (-0.44 + 0.375A^2)^2$
したがって：
[(-0.44 + 0.375A^2)^2 + 0.002304] A^2 = 0.04
この方程式を解くと、3つの実数解が現れます：
- $A \approx 0.18$（小振幅、安定）
- $A \approx 0.85$（中振幅、不安定）
- $A \approx 1.52$（大振幅、安定）
これが非線形振動の特徴！同じ周波数1.2[rad/s]でも、履歴によって振幅が変わる（ジャンプ現象）ことを示しています。
コードで確認：上のコードを実行すると、$\Omega/\omega_0 = 1.2$付近で振幅が3つ存在することが確認できます。
シミュレーターで遊ぶ：実際に動かして体感
▶ 非線形振動シミュレーター を開いて、以下の実験をやってみましょう。
実験1：ハードニングのジャンプを観察
- パラメータ設定: $\varepsilon = 0.5$, $\zeta = 0.02$, $F = 0.2$
- 周波数比 $\Omega/\omega_0$ を0.8からゆっくり増加させる
- $\Omega/\omega_0 \approx 1.15$ で、振幅が突然下の枝から上の枝にジャンプ
- 今度は1.4からゆっくり減少させる
- $\Omega/\omega_0 \approx 1.05$ で、振幅が上の枝から下の枝にジャンプ
なぜ？：周波数応答曲線が「S字」に折れ曲がっているため。増加時と減少時で異なる経路をたどるヒステリシス現象が起こります。
実験2：ソフトニングの挙動
- パラメータ設定: $\varepsilon = -0.5$, $\zeta = 0.02$, $F = 0.2$
- 周波数比を1.2から減少させていく
- 今度は低周波側でジャンプが発生！
ハードニングとソフトニングでは、ジャンプの方向が完全に逆になることを確認しましょう。
実験3：減衰の影響
- $\varepsilon = 0.5$, $F = 0.2$ で固定
- $\zeta = 0.01$（非常に小さい減衰）→ ジャンプが顕著、広い周波数範囲で複数解
- $\zeta = 0.05$（やや大きい減衰）→ ジャンプが消え、曲線がなだらかに
実務での教訓：減衰が小さいほど非線形効果が強く現れる。防振設計では、減衰材の選定がジャンプ現象の抑制に直結します。
現場でハマるポイント：適用限界と落とし穴
1. 「$\varepsilon$が大きいほど危険」は誤解
確かに$\varepsilon$が大きいと共振点が高周波側にシフトし、想定外の回転数で大振幅が発生するリスクがあります。しかし、この特性を逆手に取れば、ある周波数帯域の振動を意図的に抑制する設計も可能。
例えば、$\varepsilon=0.5$, $F=0.3$の設定では、共振ピークが右にずれて幅が狭まる。これは「振動しにくい領域」が広がったとも解釈できるんです。
2. 減衰比$\zeta$を軽視しない
線形振動では$\zeta$が小さいと共振ピークが鋭くなるだけ。しかし非線形では、ジャンプ現象の発生範囲と安定解の存在周波数幅に直接影響します。
$\zeta=0.01$から$0.05$に少し増やすだけで：
- ジャンプが起こる周波数領域が大きく変化
- 場合によってはジャンプそのものが消失
実務では材料や構造で減衰がばらつくため、シミュレーション結果を過信せず、安全マージンを十分に取ることが重要です。
3. 骨格曲線と周波数応答曲線を混同しない
骨格曲線（$\omega_0^2 + \tfrac{3}{4}\varepsilon A^2$で表されるピークの軌跡）は、減衰ゼロの理想的な共振周波数の変化を示します。
実際の応答曲線は：
- この骨格曲線を中心に、減衰によって幅を持つ
- ジャンプ現象で「折れ曲がった」形になる
設計では「まず骨格曲線の傾向を掴む → その後に減衰を加えた実際の応答を評価する」という二段階の思考が役立ちます。
4. 調和バランス法の限界
今回の計算は一次近似（応答が単一の余弦波）を仮定しています。非線形性が強い場合や共振付近では、高調波成分（$3\Omega$, $5\Omega$など）が無視できなくなり、誤差が生じることがあります。
実務では：
- 弱非線形（$\varepsilon$が小さい）→ 一次近似で十分
- 強非線形（ゴムの大変形など）→ 数値積分（Runge-Kutta法など）で直接解く必要あり
まとめ：非線形振動を「見える化」する力
非線形振動は一見複雑ですが、Duffing方程式と調和バランス法というシンプルな道具で、その本質を捉えることができます：
- 非線形項$\varepsilon x^3$が、振幅に依存した固有振動数を生み出す → ハードニング・ソフトニングの原因
- 振幅-周波数関係式が複数の解を持つ → ジャンプ現象・ヒステリシスの正体
- 骨格曲線が「理想的な共振点の変化」を示し、減衰が「実際の応答の幅と安定性」を決める
今回紹介したシミュレーターを使えば、これらの現象をブラウザ上で直感的に体験できます。パラメータを動かしながら、数式とグラフの対応を確認してみてください。
▶ 非線形振動シミュレーター — ブラウザで即動作、登録不要
NovaSolverでは1,600以上の工学シミュレーターを無料公開中 👉 一覧はこちら
「バネの硬さが変わる」という一見奇妙な現象が、実は最先端のエンジニアリングで活用されている。その面白さと奥深さを、ぜひ自分の手で確かめてみてください。