マスワークスのC MEX S-Functionのチュートリアル(2009年)でバネ・マス・ダンパ系が取り上げられていました。

 

マス、バネ、ダンパ系が例題にあるけど、これのシミュレーションをpython scriptで書ける？

Claude Opus 5:

書けます。上のスクリプトは運動方程式 m x” + c x’ + k x = F(t) を状態方程式（位置・速度の2状態）に直して解いています。関数を S-Function のコールバックに対応させてあり、derivatives() が mdlDerivatives、output() が mdlOutputs の役割です。numpy、scipy、matplotlib があれば動きます。

解き方は2通り入れています。SciPy の solve_ivp（RK45、Simulink の ode45 に相当する可変ステップ法）と、中身が見えるように自前で書いた固定ステップ RK4 です。両者の結果は図のとおりほぼ完全に重なります。m=1, c=0.5, k=4、t=1秒で1Nのステップ入力の設定では、減衰比 ζ=0.125 の不足減衰となって振動しながら、理論値の F/k=0.25 m に収束していきます。

パラメータはスクリプト冒頭の m, c, k を、入力は force(t) を書き換えれば、正弦波加振なども試せます。

理屈の部分を知りたいんだけど。なんか、この手元のチュートリアルには、行列とベクトルが出てくるけど、どういうこと？

Claude Opus 5:

チュートリアルに出てくる行列とベクトルは、**状態空間表現（state-space representation）**と呼ばれる書き方です。MathWorks の C MEX S-Function の代表的なサンプル csfunc.c もまさにこの形で、行列 A, B, C, D を使って書かれています。順を追って説明します。

1. 出発点：2階の微分方程式

マス・バネ・ダンパ系の運動方程式はニュートンの法則から

m x” + c x’ + k x = F

です。x” が入っているので「2階の微分方程式」です。ところが、ルンゲ＝クッタ法などの数値ソルバー（Simulink の ode45 も含む）は、基本的に 1階の微分方程式しか扱えません。そこで2階の式を1階の式に分解します。

2. 1階の連立方程式に分解する

新しい変数を2つ置きます。

x₁ = x（位置）
x₂ = x’（速度）

すると、

x₁’ = x₂ （位置の変化率は速度、という定義そのもの）
x₂’ = (F − c x₂ − k x₁) / m （運動方程式を x” について解いたもの）

となり、2本の1階微分方程式になりました。前回のスクリプトの derivatives() はまさにこの2行です。

3. 行列とベクトルでまとめる

この2本の式は、どちらも x₁, x₂, F の一次式（足し算と定数倍だけ）になっています。一次式の連立はそのまま行列で書けます。

d  [x₁]   [  0      1  ] [x₁]   [ 0  ]
-- [x₂] = [ -k/m  -c/m ] [x₂] + [1/m ] F
dt

記号でまとめると

x’ = A x + B u

で、

• x（状態ベクトル）＝ [x₁, x₂] = [位置, 速度]
• u（入力）＝ F
• A（システム行列）＝ [[0, 1], [−k/m, −c/m]]
• B（入力行列）＝ [0, 1/m]ᵀ

です。行列の積を展開すれば、上の2本の式に戻ることを確認できます。たとえば2行目は (−k/m)·x₁ + (−c/m)·x₂ + (1/m)·F で、x₂’ の式と一致します。

さらに「何を観測・出力するか」を

y = C x + D u

で表します。位置だけを出力するなら C = [1, 0]、D = 0 です。速度を見たければ C = [0, 1] にするだけです。

4. なぜこの書き方が便利なのか

一つ目は汎用性です。どんな線形システムも、行列の中身が違うだけで同じ形 x’ = Ax + Bu になります。モーター、電気回路（RLC回路は数学的にマス・バネ・ダンパと同じ形）、多自由度の振動系でも、S-Function のコードは「A と B を掛けて足す」だけで済みます。csfunc.c がパラメータとして行列を受け取る設計になっているのはこのためで、mdlDerivatives の中身は dx = A·x + B·u を計算するループ、mdlOutputs は y = C·x + D·u を計算するループです。

二つ目は解析ができることです。行列 A の固有値を求めると、系の性質がわかります。A の固有方程式は

λ² + (c/m) λ + k/m = 0

となり、これは運動方程式の特性方程式と同じものです。ここから固有角振動数 ωₙ = √(k/m)、減衰比 ζ = c / (2√(mk)) が出てきます。前回のスクリプトが表示した ωₙ = 2.0、ζ = 0.125 はこれです。固有値の実部がすべて負なら系は安定、虚部があれば振動する、というように、シミュレーションせずに挙動を予測できます。

三つ目は制御設計の基盤になることです。状態フィードバック、極配置、LQR、カルマンフィルタといった現代制御理論は、すべてこの状態空間表現の上に組み立てられています。チュートリアルがこの形を採用しているのは、そうした応用への入り口を意識しているからだと思います。

5. 伝達関数との関係

制御工学では同じ系を伝達関数

G(s) = 1 / (m s² + c s + k)

で表すこともあります。これは入力から出力への関係だけを見る「外から見た」表現で、状態空間表現は内部の状態（位置と速度）まで持つ「中から見た」表現です。1入力1出力の線形系なら両者は相互に変換でき、Python なら scipy.signal.ss2tf や tf2ss で行き来できます。

実際に試したい場合は、前回のスクリプトの derivatives() を

A = np.array([[0, 1], [-k/m, -c/m]])
B = np.array([0, 1/m])
def derivatives(t, x):
    return A @ x + B * force(t)

に置き換えても、まったく同じ結果になります。np.linalg.eigvals(A) で固有値を確認してみると、上の理屈がそのまま数値で見えるはずです。

今のシミュレーション。最初の動きだしは、何がきっかけ？

Claude Opus 5:

きっかけは、t = 1 秒で加わるステップ入力の力 F = 1 Nです。

スクリプトでは初期状態を x0 = [0, 0]（位置ゼロ・速度ゼロ）にしています。この状態で力もゼロなら、x₁’ = x₂ = 0、x₂’ = (0 − 0 − 0)/m = 0 となり、状態はまったく変化しません。つまり系は平衡点で静止していて、放っておけば永遠にそのままです。グラフの 0〜1 秒が平らなのはこのためです。

そこへ force(t) が t = 1 秒で 0 から 1 N に切り替わります。この瞬間に起きることを状態方程式で見ると、x₂’ = (F − c x₂ − k x₁)/m = 1/1 = 1 m/s² となり、加速度だけが突然 1 m/s² に跳ね上がります。一方、位置と速度は積分で決まる量なので、瞬間的には跳べず、その時点ではまだゼロのままです。行列で言えば、入力 F は B = [0, 1/m]ᵀ を通って速度の式（2行目）にだけ入り、位置には速度を経由して間接的に効いてくる、という構造になっています。

その結果、動き出しの直後は速度がほぼ直線的に増え、位置はおおよそ x ≈ (F/2m)(t − 1)² の放物線で滑らかに立ち上がります。グラフで t = 1 秒のところが角ばらず、なだらかに持ち上がっているのはこのためです。そこから変位が増えるにつれてバネの復元力 −k x₁ が効いてきて行き過ぎ（オーバーシュート）、戻り、ダンパで少しずつエネルギーを失いながら、力とバネが釣り合う F/k = 0.25 m に落ち着きます。

なお、動き出しのきっかけは外力だけではありません。力をゼロのまま x0 = [0.1, 0]（最初に 0.1 m 引っ張って手を離す）のように初期条件を与えても系は動き出し、今度は 0 m に向かって減衰振動する「自由振動」になります。外力による応答と初期条件による応答は、線形系では足し合わせで表せる（重ね合わせの原理）ので、両方を同時に与えるとそれぞれの応答の和になります。

マス、バネ、ダンパ系　って、自然界に実在する？例挙げて。特に、生体に関するものある？

Claude Opus 5:

純粋な「質点・理想バネ・理想ダンパ」がそのまま存在するわけではありませんが、「慣性（質量）・弾性（復元力）・粘性（摩擦的な抵抗）」の3要素が揃って振動したり減衰したりする系は自然界にいくらでもあり、それを近似するモデルとしてマス・バネ・ダンパ系が非常によく使われています。生体の例も豊富です。

身近な物理の例

風に揺れる木の枝や草は、幹の質量、材の弾性、空気抵抗と内部摩擦でまさにこの系になり、揺れが収まっていく様子は減衰振動そのものです。地震で揺れる建物も同様で、免震・制振構造の設計はこのモデルが出発点です。分子のレベルでも、化学結合の伸縮振動は原子の質量と結合のバネとしてモデル化され、赤外分光で見ている吸収は結合の固有振動数に対応しています。

生体の例

内耳（蝸牛）の基底膜は代表例です。基底膜は根元が硬く軽く、先端が柔らかく重いため、場所ごとに固有振動数が違う共振器が並んだ構造になっており、音の周波数ごとに異なる場所が大きく揺れます。Békésy のノーベル賞研究の中身がこれです。さらに有毛細胞の感覚毛の束（hair bundle）は、それ自体がバネとダンパを持つ微小な振動子で、細胞が能動的にエネルギーを注入して実効的に減衰をマイナスにすることで微弱な音を増幅している、という理論（Hopf 分岐モデル）が知られています。前回の式で言えば c を負に近づける操作に相当し、とても面白いところです。

筋肉と腱も典型で、Hill の筋モデルは収縮要素・直列弾性要素・並列弾性要素・粘性要素の組み合わせで、マス・バネ・ダンパの発展形そのものです。走ったり跳ねたりする動物の脚は「バネ付きの質点」（SLIP モデル）としてよく記述され、カンガルーやヒトのアキレス腱はバネとしてエネルギーを蓄えて再利用しています。昆虫の飛翔、特にハエやハチの間接飛翔筋は、胸郭の弾性と翅の慣性でできた共振系を固有振動数で駆動することで、神経発火の頻度をはるかに超える毎秒数百回の羽ばたきを実現しています。

眼球運動もきれいな例です。眼球と外眼筋・周囲組織からなる「オクロモーター・プラント」は、Robinson のモデル以来、2次系（マス・バネ・ダンパ）として扱われ、サッケードを生成する神経指令（パルスとステップ）はこの系の特性を打ち消すように設計されている、という説明が定番です。

神経細胞の膜にも類似の現象があります。膜そのものは抵抗とコンデンサの1次系ですが、Ih チャネルや一部の K⁺ チャネルがあると、ゆっくりした負帰還が電気的なインダクタのように振る舞い、膜全体が RLC 回路、つまり数学的にはマス・バネ・ダンパと同じ2次系になります。その結果、特定の周波数の入力に強く応答する「閾値下膜共鳴」が生じ、海馬や嗅内皮質のシータリズムとの関係が議論されています（Hutcheon & Yarom 2000 のレビューが有名です）。

ほかにも、声帯の振動、大動脈の弾性と血液の慣性による脈波、頭部衝撃時の頭蓋内での脳の揺れ（脳震盪のモデル）など、例は尽きません。

注意点：小さい世界では「マス」が消える

細胞や分子のスケールでは話が変わります。サイズが小さくなると慣性に比べて粘性が圧倒的に強くなる（低レイノルズ数の世界）ため、m x” の項は事実上ゼロになり、系は振動せずにじわっと動く過減衰の1次系になります。細胞の粘弾性を表す Kelvin–Voigt モデルや Maxwell モデル、光ピンセットに捕まったビーズの運動などは、マスのない「バネ＋ダンパ」だけのモデルで記述されます。繊毛や精子の鞭毛が慣性に頼らず、常に力を出し続けて泳いでいるのもこのためです。

つまり生体では、眼球や筋肉・脚のようなある程度大きな構造では2次系（振動しうる）、細胞以下では1次系（振動しない）、そして内耳の有毛細胞や神経膜のように能動的な仕組みで振動や共鳴を作り出している系、という3つの見方で整理すると見通しがよくなると思います。