應用有限差分法於兩條 ODE 聯立之 耦合系統數值求解
一 一以梁受載重而撓曲之分析為例
郭子瑜 · 林珈汝 · 吳南靖
一、 前言
許多在工程上考量的物理現象係以微分方程式 (Differential Equations) 描述之, 故在 工學院必修課程之工程數學中, 微分方程式佔很大比重。 在工程實務上, 因為要分析的問題通常 較為複雜, 或者計算範圍的形狀具不規則性, 或者外力條件無法用解析函數描述, 故多數問題很 難用解析的方式求解, 而必需仰賴數值方法求其近似解。 由於大學畢業學分數之門檻較過去為 低, 所以許多往昔重要的課程都已變成選修, 又因用來做力學分析之電腦軟體的功能越做越強 大, 故在力學分析與設計等應用類型課程的教學內容安排上, 有些偏重於應付國家考試的代公 式技巧, 而有些則著重在培養企業所需具戰力之軟體操作, 卻較少有能夠幫學生建立數學與力 學之間觀念連結的課程。 因此, 現今工程相關科系的學生對於數學之重要性普遍認知不足。
常微分方程式 (Ordinary Differential Equation, 簡稱 ODE) 是工程數學的第一個授課 單元, 學生們一開始接觸到工程數學時, 就被常微分方程式中各種生硬的解題技巧所困惑, 到學 習矩陣相關章節時, 或因高中就有學過矩陣、 或因不知為何而學, 而常有學習意願低落的情形, 以致後續的數值方法、 數值分析等與數學相關的應用課程往往乏人問津。 除非學生真的是很有 心, 否則就此放棄學習數學及其在工程上之相關應用。 有些甚至念到土木、 機械研究所的碩士班 了, 問他什麼是有限元素法, 他卻只能回答就是 ANSYS, 一種電腦程式, 然後就講不出其它東 西。
筆者找大學部專題生來學習並嘗試做相關的數值計算, 而後撰寫此文, 乃與陳正宗教授撰 寫 「工程數學教學拾趣 (2007)」 時所持之理念相同, 就是想激發學生的興趣與潛力。 本文以梁 受載重而產生撓曲之分析為例, 說明如何用數值方法來做常微分方程式聯立之耦合系統數值求 解, 藉以展示工程數學、 材料力學、 結構學、 數值方法、 程式設計等各課程之間的相互關聯, 希
望能讓學生在學習相關課程時, 有明確的學習目標, 減少困惑, 並且能產生學習數學的興趣, 進 而培養出自我成長的能力。
在數值方法的教科書中 (林聰悟、 林佳慧, 1997; Burden et al., 2015), 大多以單一變數 之邊界值問題來做為例說, 只解一條微分方程式。 然而在實務上, 許多問題的變數不只一個。 除 非有師徒傳承, 否則一般讀者, 即使看懂了書中所寫的數值方法, 甚至能照著書中所講的步驟求 出三維 Laplace 方程式的數值解, 要進一步應用在數條方程式聯立的多變數耦合系統, 還是有 一道待跨越的鴻溝。 本文所選之問題, 雖然自變數只有一個, 控制方程式僅為 ODE, 還沒複雜 到偏微分方程式 (Partial Differential Equation, 簡稱 PDE) 那種程度, 但由於應變數有兩 個, 是一個雙變數耦合之聯立系統, 還是比只有一條 ODE 的邊界值問題略為複雜, 恰好可做為 進階的入門範例。 而且, 要求解的變數各具物理意義, 分別是梁的內部彎矩 M(x) 和垂直方向 位移 v(x), 也可再次強調工程數學的應用性。
在眾多的數值方法中, 以有限差分法 (Finite Difference Method) 最為直接易學, 故選 之為本文用來做展示的數值方法。 希望藉由本文之例說, 能讓讀者對數值方法及其應用有進一 步的了解, 進而將來能自行應用在所需分析的問題上。
二、 有限差分法概要說明
函數 y(x) 在 x = xi 附近之泰勒級數展開式為
y(x) = y(xi) + y′(xi)(x − xi) + y′′(xi)(x − xi)2/2! + 高階項 (1) 令 yi = y(xi), yi+1= y(xi+ ∆x), yi−1= y(xi− ∆x), 則
yi+1= yi+ yi′∆x + yi′′(∆x)2/2 + O(∆x)3, (2) yi−1= yi− y′i∆x + yi′′(∆x)2/2 + O(∆x)3. (3) 若 ∆x 是一個微小量, 則 (∆x)2 會更小, (∆x)3 及更高次方項便可忽略。 把式 (2) 及式 (3) 做相加或相減, 可整理得:
y′i ≈ yi+1− yi−1
2∆x , (4)
yi′′ ≈ yi+1− 2yi+ yi−1
(∆x)2 . (5)
若令 yi+2 = y(xi+ 2∆x) 及 yi−2 = y(xi− 2∆x), 則利用類似的方法, 可得到 yi′ ≈ −yi+2+ 4yi+1− 3yi
2∆x , (6)
yi′′≈ 3yi− 4yi−1+ yy−2
2∆x . (7)
在應用時, 式 (4) 至式 (7) 中的 「近似」, 都要當成 「等於」 來看。 運用相同的原理, 可得更高階 的差分式如下:
yi(4) ≈ yi+2− 4yi+1+ 6yi− 4yi−1+ yi−2
(∆x)4 , (8)
yi′′ ≈ 2yi− 5yi+1+ 4yi+2− yi+3
(∆x)2 , (9)
yi′′ ≈ 2yi− 5yi−1+ 4yi−2− yi−3
(∆x)2 . (10)
讀者可自行參閱維基百科 (https://zh.wikipedia.org/wiki/有限差分法;
https://zh.wikipedia.org/wiki/有限差分係數) 找到更多差分公式。
三、 梁受載重撓曲分析之控制方程式與邊界條件
在材料力學的教科書中 (Hibbler, 2016; Goodno and Gere, 2009), 必能找到以下兩條 常微分方程式, 用以描述梁受載重時 (如圖 1 示意), 在彈性範圍內之彎矩與位移之分佈。
d2M
dx2 = −q, d2v
dx2 = M EI,
(11)
其中, M = M(x), 為梁之彎矩, 定義梁頂受壓為正, 受拉為負; x 為水平方向之坐標; q = q(x) 為梁所受之分佈載重, 定義向下為正, 向上為負; v = v(x) 為梁受力變形後之垂直方向位移, 定 義向上為正, 向下為負; EI = EI(x) 則為材料性質與斷面特性。 若梁為均質材料構成, 且全梁 斷面形狀一致, 則 EI為一個定值。
圖 1: 梁受載重之示意圖
常見的邊界條件, 有表 1 所列之幾種型式。 表中所列為右端的邊界條件, 這些邊界條件同 樣適用於左端。
式 (11) 可再合併成為一條四階 ODE, 如下:
d4v
dx4 = − q
EI. (12)
表 1: 梁受載重撓曲分析之常見支承型式及其邊界條件與物理意義
支承型式 邊界條件 物理意義
鉸支承 M = 0
v = 0 不受彎矩, 可旋轉, 無位移。
滾支承 M = 0 不受彎矩, 可旋轉, 無垂直位移,
v = 0 但可有水平位移。
自由端 M = 0 不受彎矩, 可旋轉, 可位移, 當端點受
M′ = given 集中載重時, M′ 與該載重大小有關。
固定端 v = 0
v′ = 0 可受彎矩, 不可旋轉, 無位移。
若是要用解析的方法求解此類問題, 只解一條帶有一個未知函數的 ODE, 即使微分有到四階, 還是比解一個兩條 ODE 聯立的雙變數耦合系統還容易。 然而, 就如前面所言, 實務上要分析 的問題通常複雜到無法用解析的方式求解, 故數值求解還是有其重要性。 而在數值求解上, 則剛 好反過來, 解高階 ODE, 所得到的結果在準確度上往往反而較差。
四、 利用有限差分法做此類問題數值求解之步驟
利用有限差分法做此類問題數值求解之步驟如下:
(1) 首先, 我們用 N 個結點, 把計算域 0 ≤ x ≤ L 分成 (N − 1) 等分, 每一等分之長度為
∆x。 第一個結點, 為 x = 0; 最後一個結點, 為 x = L。
(2) 在內部結點, 即 i = 2, . . . , N − 1, 控制方程式可離散並整理成
Mi−1− 2Mi+ Mi+1= −qi(∆x)2,
−(∆x)2 (EI)i
Mi+ vi−1− 2vi+ vi+1= 0, (13) 其中, qi 和 (EI)i 分別為 q 和 EI 在 x = xi 處的值。
此兩條差分式的等號左邊均為未知, 而等號右邊則均為已知。
(3) 在邊界點, 即 i = 1 及 i = N, 代邊界條件。
(4) 由步驟 2 和 3, 可得到一個共有 2N 條 2N 元一次聯立方程式的線性系統。 解聯立方程組, 可得 M1 ∼ MN 及 v1 ∼ vN 之近似值。
(5) 利用求得的數值解展示 M 和 v 的分佈。
若是要用有限差分法解一條四階 ODE, 則需要用式 (8) 至式 (10) 所列之差分式去形成 一個 N 元一次聯立方程式的線性系統。 求解得到 v1 ∼ vN 近似值之後, 再利用差分公式求出 M1 ∼ MN 之近似值。 由於步驟相當類似, 故在此不再贅述。
五、 計算範例
(1) 靜定問題
首先, 我們取一個簡支梁受均佈載重而撓曲之靜定問題來做展示, 如圖 2。 取左端支承點 當原點, 則這個問題的計算域為 0 ≤ x ≤ 5, x 的單位為公尺。 兩條控制方程式為
d2M
dx2 = −4, d2v
dx2 = M 500.
(14)
在左端為鉸支承, 可以有轉動, 但不能有位移, 也不能受彎矩, 故其邊界條件為 M(0) = 0 及 v(0) = 0。 在右端為滾支承, 可以有水平位移, 但不能有垂直位移, 也一樣是可以有轉動, 但不 能受彎矩, 故其邊界條件為 M(5) = 0 及 v(5) = 0。
圖 2: 簡支梁受均佈載重之撓曲分析靜定問題
直接做積分, 並代邊界條件, 可得到這個問題的正解
M(x) = −2x2+ 10x, (15)
及
v(x) = 8 × 10−3×
− x4
24+ 5x3
12 − 125x 24
, (16)
其中, 彎矩 的單位是 kN-m, 而垂直方向位移 v 的單位則為 m。
我們把計算域切成 5 個等分, 間隔 ∆x 為 1m, 共有 6 個結點, 分別為 x1 = 0、 x2 = 1、. . .、 x6 = 5。 然後依據前述之邊界條件, 可得以下 4 條方程式。
M1 = 0、 v1 = 0、 M6 = 0、 v6 = 0, (17) 再依前面所述之差分公式, 可得另外 8 條方程式。
M1 − 2M2+ M3= −4、 − 1
500M2+ v1− 2v2+ v3 = 0、
M2 − 2M3+ M4= −4、 − 1
500M3+ v2− 2v3+ v4 = 0、
M3 − 2M4+ M5= −4、 − 1
500M4+ v3− 2v4+ v5 = 0、
M4 − 2M5+ M6= −4 及 − 1
500M5+ v4− 2v5+ v6 = 0。
(18)
至此, 讀者可能覺得, ∆x = 1, (∆x)2 和 (∆x)3 也都是 1, 怎麼會說 (∆x)3 以上的高次 方項可以忽略? 在此, 我們要強調這是相對的概念, 雖然 ∆x、 (∆x)2 和 (∆x)3 都是 1, 但相對 之下, ∆x 是計算域的長度 L 的 1/5, (∆x)2 則是計算域的長度平方 1/25, (∆x)3 是計算域 的長度三次方的 1/125, 高次方項還是可以忽略, 即使只把計算域切成 5 等分, 計算結果還是不 會太差。 這些在接下來的計算結果會得到驗證。
式 (17) 及式 (18) 總共列了 12 條一次方程式, 共有 12 個未知數, 將之聯立並求得數值 解, 和正解比較如下:
i x Mi vi
數值解 正解 數值解 正解
1 0 0.000000 0.000000 0.000000 0.000000 2 1 8.000000 8.000000 -0.040000 -0.038667 3 2 12.000000 12.000000 -0.064000 -0.062000 4 3 12.000000 12.000000 -0.064000 -0.062000 5 4 8.000000 8.000000 -0.040000 -0.038667 6 5 0.000000 0.000000 0.000000 0.000000
以上只顯示到小數點以下六位。 讀者可發現, 在彎矩 M 的部分, 數值解和正解完全一致, 而在 垂直方向位移 v 的部分, 數值解和正解之間有一些小小誤差, 這是因為我們只把計算域切成 5 等分而已。 相信一般人不會有耐心去手算解一個 12 元一次方程式的聯立系統, 一定是把線性 聯立方程組變成一個矩陣式, 然後再看是要用軟體去求解, 或是另外寫程式用高斯消去法或其 他方法把這 12 個未知數求出來。 在工程數學的教科書裡 (Kreyszig, 2018; O’neil, 2018), 有 說明線性聯立方程組的係數各別在矩陣裡的哪個位置, 也有說明如何利用數值方法解線性聯立 方程組, 但其實這些在高中的數學課程裡就有了。 係數矩陣裡共有 144 個數, 很少人有耐心把 數字一個一個打進去, 最好也是能寫一個程式去產生線性聯立方程組的係數矩陣及矩陣式等號 右邊的向量。
這個計算例, 由於兩端邊界條件都有 M = 0, 故其實可先解第一條 ODE 得 M 的分佈, 再解第二條 ODE 得 v 的分佈。 不過, 這種解法並不具一般性, 若要分析的是靜不定梁的撓曲 問題, 則不適用。 如果把計算域切成 10 等分, 即 ∆x = 0.5, 共有 11 個結點, 則在形成的 22 條一次方程式的聯立系統中, 中間結點的差分式會變成
Mi−1− 2Mi+ Mi+1= −4 × 0.52,
−0.52
500Mi+ vi−1− 2vi+ vi+1= 0,
i = 2, . . . , 10, (19)
結點數變得更多, 連手寫列出聯立方程組都不可能, 何況手算去解聯立方程組? 一定要寫程式。
如果要以數值方法求解兩條二階 ODE 合併而成的四階 ODE, 則只會形成一個 11 條一 次方程式的聯立系統。 由邊界條件得到的方程式如下:
v1 = 0、 2v1− 5v2+ 4v3− v4 = 0、 − v8+ 4v9 − 5v10+ 2v11= 0、 v11= 0, (20) 而中間結點由差分式得到的方程式如下:
vi−2− 4vi−1+ 6vi − 4vi+1+ vi+2= −4 × 0.54/500, i = 3, . . . , 9, (21) 求出各結點 v 的近似值之後, 由 M = EId2v
dx2, 在邊界點利用式 (9) 或式 (10) 的其中一條, 而在內部點則利用式 (5), 可得到 M1 ∼ M11 之近似值。 求解得到的結果和正解比較如下:
i x Mi vi
正解 數值解 1 數值解 2 正解 數值解1 數值解 2 1 0.0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2 0.5 4.500000 4.500000 3.500000 -0.020438 -0.020625 -0.018375 3 1.0 8.000000 8.000000 7.000000 -0.038667 -0.039000 -0.035000 4 1.5 10.500000 10.500000 9.500000 -0.052938 -0.053375 -0.048125 5 2.0 12.000000 12.000000 11.000000 -0.062000 -0.062500 -0.056500 6 2.5 12.500000 12.500000 11.500000 -0.065104 -0.065625 -0.059375 7 3.0 12.000000 12.000000 11.000000 -0.062000 -0.062500 -0.056500 8 3.5 10.500000 10.500000 9.500000 -0.052938 -0.053375 -0.048125 9 4.0 8.000000 8.000000 7.000000 -0.038667 -0.039000 -0.035000 10 4.5 4.500000 4.500000 3.500000 -0.020438 -0.020625 -0.018375 11 5.0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 其中, 數值解 1 即為求解兩條二階 ODE 聯立的結果; 而數值解 2, 則為把兩條二階 ODE 合 併成一條四階 ODE, 再做數值求解的結果。 由兩種數值解的比較, 可知解一條四階 ODE 的結 果顯然劣於解兩條二階 ODE 聯立, 甚至比前面只有用 6 個結點的計算結果還差。 用了較多結 點之後, 垂直方向位移 v 的數值解和正解之間的差異有變得更小。
(2) 靜不定問題
接下來, 我們把前一例題的載重情形, 改成只局部作用在 2.5 ≤ x ≤ 5, 並且把左端的鉸 支承改成固定端, 如圖 3 所示, 則變成靜不定問題。 所謂的靜不定, 是指無法直接由力平衡方程 式求出支承點的反作用力, 而必須加入材料變形的特性才能進行求解並分析的問題。 這個問題 中的兩條控制方程式為
d2M
dx2 = −4u(x − 2.5), d2v
dx2 = M 500,
(22)
其中,
u(x − 2.5) =
(0, 0 ≤ x < 2.5,
1, 2.5 < x ≤ 5, (23) 稱為單位階梯函數, 在 x = 2.5 處, 為不連續點, 實務上在該處設定其函數值為 0.5。 在邊界條
件方面, 左端為固定端, 不能有轉動, 也不能有位移, 彎矩為未知, 故其邊界條件為 v(0) = 0 及 v′(0) = 0。 在右端的邊界條件, 則和前一例相同, 為 M(5) = 0 及 v(5) = 0。
圖 3: 懸臂梁右端加額外支承之局部受均佈載重撓曲分析靜不定問題
這個問題的正解是
M(x) =
4×
−7×2.52
32 +23×2.5 64 ×x
, 0 ≤ x ≤ 2.5
4× 9
64×2.52+23
64×2.5×(x−2.5)−1
2×(x−2.5)2
, 2.5 ≤ x ≤ 5 (24)
及
v(x) =
4×
−7×2.52×x2
64 +23×2.5×x3 384
, 0 ≤ x ≤ 2.5
4×19 × 2.54
384 −5×2.53×(x−2.5)
128 +9×2.52×(x−2.5)2 128
+23×2.5×(x−2.5)3
384 −(x−2.5)4 24
,
2.5 ≤ x ≤ 5 (25)
要求出這個問題的正解, 程序相當地麻煩, 在此略過, 有興趣的讀者可自行參考結構學教科書 (Hibbeler, 2019) 裡的最小功法、 單位力法或傾角變位法, 先求出兩個邊界點的反作用力, 再 利用力平衡算出彎矩分佈, 然再用積分把位移分佈求出來。 這個問題也無法分兩步先解 M 再 解 v, 不是 M 和 v 要一起解, 就是要解四階 ODE 得 v。 由於前面的計算例已證明把兩條二 階 ODE 合成一條四階 ODE 再進行求解並不佳, 故在此不再多做比較。
和前例一樣, 我們把計算域切成 10 等分, 即 ∆x = 0.5, 共有 11 個結點, 則在左邊界, 可 列出以下兩條方程式
v1 = 0 及 − 3v1 + 4v2− v3 = 0, (26) 在右邊界, 可列出以下兩條方程式
M11 = 0 及 v11 = 0. (27)
在內部點, 可列出方程式如下
Mi−1− 2Mi+ Mi+1= 0
−0.52
500Mi+ vi−1− 2vi+ vi+1= 0 , i = 2, . . . , 5, (28)
Mi−1− 2Mi+ Mi+1= −4 × 0.5 × 0.52
−0.52
500Mi+ vi−1− 2vi+ vi+1= 0 , i = 6, (29)
Mi−1− 2Mi+ Mi+1= −4 × 0.52
−0.52
500Mi+ vi−1− 2vi+ vi+1= 0 , i = 7, . . . , 10. (30) 如此總共列出 22 條一次方程式, 共有 22 個未知數, 然後聯立求解而得到數值解。 與正解比較 如下:
i x Mi vi
數值解 正解 數值解 正解
1 0.0 -5.681818 -5.468750 0.000000 0.000000 2 0.5 -3.863636 -3.671875 -0.000966 -0.001217 3 1.0 -2.045455 -1.875000 -0.003864 -0.004271 4 1.5 -0.227273 -0.078125 -0.007784 -0.008262 5 2.0 1.590909 1.718750 -0.011818 -0.012292 6 2.5 3.409091 3.515625 -0.015057 -0.015462 7 3.0 4.727273 4.812500 -0.016591 -0.016896 8 3.5 5.045455 5.109375 -0.015761 -0.015965 9 4.0 4.363636 4.406250 -0.012409 -0.012521 10 4.5 2.681818 2.703125 -0.006875 -0.006915 11 5.0 0.000000 0.000000 0.000000 0.000000
以上所列, 數值解與正解之間有一些小小誤差, 但若再將計算域切得更細, 用更多結點來描述之, 則得到的數值解會更接近正解。 不過, 因數字太多, 故在此略過, 僅以繪圖展示計算結果與正解 之比較, 如圖 4 及圖 5 所示。 由比較結果可知, 用 21 個結點將計算域切成 20 等分, 所得到的 數值解與正解近乎完全吻合。 一般實例應用的問題, 往往比這問題複雜許多, 比如載重的分佈不 具規則性, 或比如梁的斷面形狀或尺寸隨位置而變, EI 不為一個定值, 則無法求出正解, 唯有 用數值方法才能得到分析結果, 由此亦可知數值方法的價值所在。
圖 4: 靜不定撓曲分析之彎矩分佈數值解與正解比較
圖 5: 靜不定撓曲分析之位移分佈數值解與正解比較
六、 結論與建議
本文選取梁受載重而撓曲之分析, 來說明兩條 ODE 聯立之耦合系統數值求解之步驟。 所 選的例題, 較一般數值方法的教科書中的單一變數之邊界值問題略為複雜, 但又沒複雜到 PDE 那種程度, 恰好可做為介紹數值方法的進階範例。 而且, 要求解的變數各具物理意義, 同時也強 調了工程數學的應用性。 因有限差分法為最直接易學的數值方法, 故選之為本文用來做展示的 數值方法。
本文先對有限差分法做概述, 然後介紹梁受載重撓曲分析問題要求解的變數、 控制方程式, 以及邊界條件, 接下來說明利用有限差分法做此類問題數值求解之步驟, 最後分別以靜定及靜
不定的計算例做展示。 在兩個計算例中, 都顯示若把計算域切得越細, 計算結果越接近正解。
本文所選問題中的兩條二階 ODE 可合併為一條四階 ODE。 直覺上只解這條合併得到的 四階 ODE 應該比解兩條 ODE 聯立的雙變數耦合系統還容易, 然經實例證明, 所得計算結果 之準確度較差, 故較不建議採用。
由於求解的步驟裡, 最後都會形成一個含有許多未知數的線性聯立方程組, 必須寫成矩陣 式, 再給電腦去求解, 才會有效率。 由此也可看出培養程式撰寫能力的重要性。 第二個計算例具 較高的複雜度, 若手算必須花費極大心力, 但若採用數值方法, 則在第一個計算例所形成的矩陣 式中, 只要改變若干位置的數值, 再給電腦解聯立方程組, 即可很快算出極為接近正解的答案。
工程實務上的問題, 又常比本文所舉的計算例更為複雜, 往往只能求近似解, 由此亦可顯示出數 值方法的價值所在。
