用解耦(Decoupling)概念求解雙木塊耦合震盪問題(含程式模擬)

這篇文章想要討論一個物理問題， 雙木塊耦合震盪問題(Coupled Oscillator)，想要了解這個問題，需要用常微分方程(Ordinary Differential Equation, ODE)、普通物理、複變，線性代數。

我會從一個有趣的脈絡來切入這個問題，用線性代數中的「Decoupling (解耦)」的概念來解這個問題，也會用Python程式語言來模擬物理現象，希望透過這個問題讓大家了解線性代數理論的巧妙。

• 以下是雙彈簧耦合震動問題的意示圖

單木塊彈簧問題

求解雙木塊耦合震盪問題(Coupled Oscillator)之前，我們先從比較簡單的單木塊彈簧問題下手，由從虎克定律(Hooke's Law)木塊的受力在彈性限度範圍內，受到的力與彈簧的形變成線性關係，再根據牛頓第二運動，單木塊彈簧問題滿足以下ODE：

其中是木塊的重量，表示彈簧的彈性係數，且 是木塊距離靜力平衡點(Static Equilibrium)的位移。

微分方程式描述物理現象為：

木塊受力會與木塊與靜力平衡點的位移成負號的關係，且距離越遠(大)，受到的反方向的力也越大。

• 單木塊彈簧問題意示圖

為了呼應後面的我們的討論，我們令，式改寫為：

回想最簡單的ODE，的解為指數函數 ( 為積分常數)。

想要求解，可「猜」解也是指數的形式，且

將 帶入式：

滿足以上等式有兩個可能， 或 ，

是我們「猜」出來的解，我們希望解不要等於0。事實上對任意的齊次微分方程，所有都是一個簡單解(trivial solution)，這個解太過簡單，我們不在乎這個解，因為 的任意微分都還是，永遠等於齊次微分方程的等式右邊的。

由於上述原因，我們考慮一般的狀況，也就是當的時的狀況

這個方程也稱為微分方程 的特徵方程式 (Characteristic Equation)，透過二次方程式的公式解，得到兩個根，，這裡注意 Hooke's Law 指出彈簧的受力與位移成負號的關係，所以。

且，所以特徵方程式的根是兩個複數根 ，將複數根代回得到的兩個解 。

因為 除上 不是一個常數，這兩個解之間為線性獨立。

注意，指數的虛數次方是透過 Euler Formula 來計算，所以兩個線性獨立解均為複數。

在ODE裡面有個重要的定理，叫做重疊定理(Superposition Theorem)：

對於齊次的常微分方程(ODE)，所有解都可以透過線性獨立解的線性組合所生成。

重疊定理的證明不困難，我們可以簡易說明如下：

齊次ODE的解集合實際上就是一個subspace，因為子空間的定義，子空間中的線性獨立向量的線性組合，還會落在該子空間中。用線性代數的語言來說，齊次ODE解集合的是一個kernel space，kernel space也是線性變換四個基本子空間中的一個。

透過重疊原理，如果我們找出方程式的兩個線性獨立解，可以透過將做線性組合得到式的通解(General Solution)

其中

又因為 的物理意義是木塊距離靜力平衡點(Static Equilibrium)的位移，位移必須實數，且均為複數，所以係數 需要是複數才可能讓是實數，才符合物理量的意義，我們把這個條件稱為「實條件」 (Realily)。

「實條件」(Realily)：如果一個複數，他是實數，則他會與自己的共軛複數相等，也就是

是實數，根據「實條件」須滿足條件，則

其中利用到

由比較係數得到，與

根據「實條件」，我們得到式的通解(General Solution)需滿足以下形式

其中

用Euler Formula 且假設 代入

令，我們進一步將式的通解(General Solution)改寫成:

其中

將位移函數對時間微分，會得到木塊的速度對時間的關係：

其中

因為有 兩個變數，所以需要給兩個初始條件才能計算，因為，則

所以只要給定木塊的初始位置 與初始速度 就能刻畫出單木塊彈簧問題的解：

我們用 Python 程式語言模擬單木塊彈簧的動畫，用到的套件為 Vpython 與 Numpy

給定初始參數為

彈簧模擬動畫與 x-t 圖曲線

程式碼如下: 程式的主要邏輯是用推導出來的解, 來模擬木塊與彈簧的行為

# coding=utf-8
from vpython import *
import numpy as np

"""
 1. 參數設定, 設定變數及初始值
"""
size = 0.1   # 木塊邊長
L = 1       # 地板長度

k = 3
m = 1
x0 = 0  
x0_dot = 0.3


t = 0        # 時間
T = 30     
dt = 0.01    # 時間間隔


"""
 2. 畫面設定
"""
scene = canvas(title=f"單木塊彈簧，假設地版無摩擦 m={m}, k={k}", width=800, height=300, x=0, y=0, center=vec(0, 0.1, 0), background=vec(0, 0.6, 0.6))
floor = box(pos=vec(0, 0, 0), size=vec(L, 0.1*size, 0.5*L), color=color.blue)

wall1 = box(pos=vec(-L/2, 0.3*L/2, 0), size=vec(0.1*size, 0.3*L, 0.5*L), color=color.blue)

cube1 = box(pos=vec(x0, 0.55*size, 0), size=vec(size, size, size), color=color.red, v=vec(0, 0, 0))

spring1 = helix(pos=vec(-L/2,  0.55*size, 0), radius=0.3*size, thickness=0.05*size, color=color.yellow)
spring1.axis = cube1.pos - spring1.pos - vec(size/2, 0, 0)

gd = graph(title="x-t plot", width=600, height=450, x=0, y=600, xtitle="t(s)", ytitle="x(m)", label="cos(x)")
gd2 = graph(title="v-t plot", width=600, height=450, x=0, y=1050, xtitle="t(s)", ytitle="v(m/s)")
xt = gcurve(graph=gd, color=color.red)
vt = gcurve(graph=gd2, color=color.red)


"""
 3. 物體運動部分, 木塊到達地板邊緣時停止執行
"""

while(t<T):
    rate(100)
    cube1.pos.x = - 0.2*L + x0*(np.cos(np.sqrt(k/m)*t)) \
                           + (x0_dot/np.sqrt(k/m))*np.sin(np.sqrt(k/m)*t)

    cube1.v.x = -x0*(np.sqrt(k/m))*np.sin(np.sqrt(k/m)*t) \
                + x0_dot*np.cos(np.sqrt(k/m)*t)
    
    spring1.axis = cube1.pos - spring1.pos - vec(size/2, 0, 0)


    xt.plot(pos=(t, cube1.pos.x))
    vt.plot(pos=(t, cube1.v.x))
    t += dt

print("Done")

雙木塊耦合震盪問題(Coupled Oscillator)

我們有兩個質量為m的木塊，有兩面固定不動的牆，木塊之間，木塊與牆壁之間用彈性係數為的彈簧連接。

雙木塊耦合震盪問題的意示圖如下圖：

因為虎克定律(Hooke's Law)，木塊的受力在彈性限度範圍內，受到的力與彈簧的形變成線性關係，假設向右的力量為正。

假設「左」邊的木塊往右位移了，「右」邊的木塊往右位移了 ，如果的值為負，表示該木塊向左移動，位移量表示整個系統處於靜力平衡的狀態。

假設兩個木塊分別向右移動，木塊移動前的位置用虛線表示，移動後的位置用實線表示，見下圖：

在以上的意示圖中，兩個木塊位移的關係式為

• 觀察「左」木塊:

• 對於「左」木塊的左邊的彈簧，也就是被拉伸了的位移，根據Hooke's Law 會產生向左的力量，因為向左所以記做

• 對於「左」木塊的右邊的彈簧(整個系統中間的彈簧)，可以考慮以下幾個情況

• 如果，也就是把「右」木塊當作一個不動的牆壁，受力為。

• 如果，對於「左」木塊的右邊彈簧，彈簧沒有形變，根據Hooke's Law 不會產生力。

• 如果，對於「左」木塊的右邊彈簧，彈簧被壓縮了，根據Hooke's Law 會產生將「左」木塊往「左」的力量，力量大小為，考慮向右為正，記做。

• 如果，對於「左」木塊的右邊彈簧，彈簧被拉伸了，根據Hooke's Law 會產生將「左」木塊往「右」的力量，力量大小為，考慮向右為正，記做。

總和以上討論

「左」木塊受的力量總合為 ，根據牛頓第二運動定律

• 觀察「右」木塊:

與「左」木塊的討論類似，不過整個系統中間的彈簧受力的方向與「左」木塊的情況相反，這裡留給讀者自己推導

「右」木塊受的力量總合為 ，根據牛頓第二運動定律

雙木塊耦合震盪問題，可以寫成二階齊次ODE方程組：

1. 透過變數變換求解雙木塊耦合震盪問題

這裡我們用一個技巧，先把的兩式加得到

接著再將的上式減下式

令

式 可以改寫成

透過這個變數變換，原本的耦合系統 (Coupled System)，被解耦(Decoupling)了!!

方程組中，兩個式子均滿足我們前面討論的單木塊彈簧問題，只要給定初始條件，我們就可以用刻劃出的解，方程組的解為：

但是我們並不想要解，因為的物理意義並不是我們原先定義的問題中兩個木塊的位移,，我們想要的解是 Coupled Oscillator 的解

由，上式加下式除二，上式減下式除二，我們也能夠將在寫回的函數

把代入

再把代回，且微分是線性，所以

Coupled Oscillator 的解可以寫成：

速度，也就是位移函數的微分，透過對微分得

觀察解，如果再以下兩種特別的初始值，Coupled Oscillator 的解會產生兩個最「簡單」的狀態：

• 當：

Coupled Oscillator的解簡化為

• 當：

Coupled Oscillator的解簡化為