數(shù)值求解波動(dòng)方程 [1]

2022-08-16 16:13 作者:nyasyamorina 0人讀過 | 我要投稿

????????關(guān)于這個(gè)系列的專欄,? 與其說是教程,? 不如說更像是筆記,? 因?yàn)閿?shù)值求解的內(nèi)容實(shí)在在太多太多了,? 我也只是學(xué)到什么寫什么而已,? 希望可以幫到有需要的人?

波動(dòng)方程

????????波動(dòng)方程是一個(gè)用于描述像弦振動(dòng),? 薄膜振動(dòng),? 聲波傳播和地震波等現(xiàn)象的方程,? 甚至麥克斯韋方程組也可以退化為波動(dòng)方程,? 從而使用波動(dòng)方程去研究電磁行為.

????????一般地,? 波動(dòng)方程的形式為? $%5Cfrac%7B%5Cpartial%5E2%20u%7D%7B%5Cpartial%20t%5E2%7D%3Dc%5E2%5CDelta%20u$ ,? 其中 $%5CDelta%3D%5Cfrac%7B%5Cpartial%5E2%7D%7B%5Cpartial%20x_1%5E2%7D%2B%5Cfrac%7B%5Cpartial%5E2%7D%7B%5Cpartial%20x_2%5E2%7D%2B%5Ccdots$ ?為拉普拉斯算子,? $c$ 為波在介質(zhì)里傳播的速度,? 又叫相速度.? 那么求解函數(shù)? $u(x_1%2Cx_2%2C%5Ccdots%2Ct)$ ?即為所求波的解.? 關(guān)于波動(dòng)方程的推導(dǎo),? 數(shù)學(xué)解等,? 在《數(shù)理方程》里就有,? 而且也不是這個(gè)專欄的重點(diǎn).? 因?yàn)椴▌?dòng)方程是二階微分方程,? 所以在給出初值條件 (一般來說是 $t%3D0$ ?的 $u$ ?和 $%5Cpartial%20u%2F%5Cpartial%20t$ ) 時(shí)可以求出相應(yīng)的特解.? 另外波動(dòng)方程具有很多種變形:? 邊值條件,? 強(qiáng)迫振動(dòng),? 非均勻介質(zhì)等.

????????數(shù)值求解波動(dòng)方程就是在已知初值條件時(shí)使用計(jì)算機(jī)求得相應(yīng)的特解,? 主要實(shí)現(xiàn)方法有有限元法和有限差分法等.? 相較于有限元法,? 有限差分法更容易,? 所以這系列的專欄是使用有限差分法進(jìn)行波動(dòng)方程的數(shù)值求解.

有限差分法

????????由于計(jì)算機(jī)無法處理無限的數(shù)據(jù),? 所以當(dāng)進(jìn)行數(shù)值模擬時(shí)必須對(duì)模型進(jìn)行離散化 (避免無限小的間隔) 和截?cái)?(避免無限大的區(qū)間).

????????由導(dǎo)數(shù)的定義?? $%5Cfrac%7Bdf%7D%7Bdx%7D(x)%3D%5Clim_%7Bdx%5Crightarrow0%7D%5Cfrac%7Bf(x%2Bdx)-f(x)%7D%7Bdx%7D$ ?可以得出有限差分法的離散化為 $%5Cfrac%7Bf(x%2B%5CDelta%20x)-f(x)%7D%7B%5CDelta%20x%7D%5Capprox%5Cfrac%7Bdf%7D%7Bdx%7D(x)%2C%5C%2C%5CDelta%20x%3E0$ ,? 但 $%5Clim_%7B%5CDelta%20x%5Crightarrow0%7D$ ?暗含著條件 $%5Clim_%7Bdx%5Crightarrow0%5E-%7D%3D%5Clim_%7Bdx%5Crightarrow0%5E%2B%7D$ ,? 所以有 $%5Cfrac%7Bf(x)-f(x-%5CDelta%20x)%7D%7B%5CDelta%20x%7D%5Capprox%5Cfrac%7Bf(x%2B%5CDelta%20x)-f(x)%7D%7B%5CDelta%20x%7D$ ,? 于是更準(zhǔn)確的離散化應(yīng)為 $%5Cfrac%7Bdf%7D%7Bdx%7D(x)%5Capprox%5Cfrac%7B1%7D%7B2%7D%5Cleft(%5Cfrac%7Bf(x%2B%5CDelta%20x)-f(x)%7D%7B%5CDelta%20x%7D%2B%5Cfrac%7Bf(x)-f(x-%5CDelta%20x)%7D%7B%5CDelta%20x%7D%5Cright)%3D%5Cfrac%7Bf(x%2B%5CDelta%20x)-f(x-%5CDelta%20x)%7D%7B2%5CDelta%20x%7D$ .? 類似地可以得出二階導(dǎo)數(shù)的離散化:?? $%5Cfrac%7Bd%5E2f%7D%7Bdx%5E2%7D(x)%5Capprox%5Cfrac%7Bf(x-%5CDelta%20x)-2f(x)%2Bf(x%2B%5CDelta%20x)%7D%7B%5CDelta%20x%5E2%7D$ .? 雖然可以繼續(xù)推到更高階導(dǎo)數(shù)的離散化,? 但是波動(dòng)方程只需二階就足夠了.

????????函數(shù)在離散化后可以使用數(shù)組儲(chǔ)存,? 假設(shè)? $f(x)$ ?在數(shù)組表示為 $f_x$ ,? 其中下標(biāo)為索引,? 那么有 $f(x%2B%5CDelta%20x)$ ?在數(shù)組中為 $f_%7Bx%2B1%7D$ .? 所以導(dǎo)數(shù)的離散化可以寫為 $%7B%5Cfrac%7Bdf%7D%7Bdx%7D%7D_x%3D%5Cfrac%7Bf_%7Bx%2B1%7D-f_%7Bx-1%7D%7D%7B2%5CDelta%20x%7D$ ,? 二階導(dǎo)數(shù)為 $%7B%5Cfrac%7Bd%5E2f%7D%7Bdx%5E2%7D%7D_x%3D%5Cfrac%7Bf_%7Bx-1%7D%2Bf_%7Bx%2B1%7D-2f_x%7D%7B%5CDelta%20x%5E2%7D$ .

????????有限差分法的區(qū)間截?cái)嗟葍r(jià)于微分問題里的邊值條件,? 通過控制邊界的行為去避免讀取區(qū)間外的數(shù)據(jù),? 從而使計(jì)算范圍限制在區(qū)間內(nèi).? 常用的邊值條件有 Dirichlet 條件和?Neumann 條件,? 假設(shè)? $x_1$ ?為邊界,? 那么 Dirichlet 條件表示為 $f(x_1)%3Dk$ ,? Neumann 條件表示為 $%5Cfrac%7Bdf%7D%7Bdx%7D(x_1)%3Dk$ ,? 其中 $k$ ?為常數(shù).

????????在使用數(shù)值求解時(shí)需要注意一個(gè)嚴(yán)重的問題:? 數(shù)值不穩(wěn)定.? 它會(huì)造成求解誤差較大或者直接發(fā)散出現(xiàn) inf, nan.? 由于我自己并沒有太深究這個(gè)問題,? 所以并沒有具體分析過數(shù)值不穩(wěn)定,? 但從個(gè)人經(jīng)驗(yàn)來說,? 數(shù)值不穩(wěn)定與參數(shù)選擇和模型的具體實(shí)現(xiàn),? 甚至初值有關(guān).

實(shí)現(xiàn)

????????下面以一維波動(dòng)方程為例,? 一維波動(dòng)方程為? $%5Cfrac%7B%5Cpartial%5E2u%7D%7B%5Cpartial%20t%5E2%7D%3Dc%5E2%5Cfrac%7B%5Cpartial%5E2u%7D%7B%5Cpartial%20x%5E2%7D$ ,? 其中函數(shù) $u(x%2Ct)$ ?即為所求的解.? 方程離散化得? $%5Cfrac%7B1%7D%7B%5CDelta%20t%5E2%7D(u_%7Bx%2Ct-1%7D%2Bu_%7Bx%2Ct%2B1%7D-2u_%7Bx%2Ct%7D)%3D%5Cfrac%7Bc%5E2%7D%7B%5CDelta%20x%5E2%7D(u_%7Bx-1%2Ct%7D%2Bu_%7Bx%2B1%2Ct%7D-2u_%7Bx%2Ct%7D)$ ,? 移項(xiàng)得? $u_%7Bx%2Ct%2B1%7D%3D%5Cfrac%7Bc%5E2%5CDelta%20t%5E2%7D%7B%5CDelta%20x%5E2%7D(u_%7Bx-1%2Ct%7D%2Bu_%7Bx%2B1%2Ct%7D)-u_%7Bx%2Ct-1%7D%2B2%5Cleft(1-%5Cfrac%7Bc%5E2%5CDelta%20t%5E2%7D%7B%5CDelta%20x%5E2%7D%5Cright)u_%7Bx%2Ct%7D$ .? 其中兩塊系數(shù)可以提前計(jì)算以提升效率.

????????然后需要設(shè)置邊值條件,? 這里在左邊界使用?Dirichlet 條件 $u(%5Ccdot%2Ct)%3D0$ ,? 右邊界使用?Neumann 條件? $%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20x%7D(%5Ccdot%2Ct)%3D0$ ,? 假設(shè)離散化后的數(shù)組索引范圍為 $%5B1%2Cl%5D$ ?(大部分編程語言是 $%5B0%2Cl-1%5D$ ,? 我都老 julia 人了),? 那么根據(jù)邊值條件有? $u_%7B0%2Ct%7D%3D0$ ?和 $u_%7Bl%2B1%2Ct%7D%3Du_%7Bl-1%2Ct%7D$ ?(對(duì)于其他語言是 $u_%7B-1%2Ct%7D%3D0$ 和 $u_%7Bl%2Ct%7D%3Du_%7Bl-2%2Ct%7D$ ).

????????最后因?yàn)椴▌?dòng)方程是二階微分方程,? 所以至少需要提供兩個(gè)初值條件,? 這里以? $u(%5Ccdot%2C0)$ 和? $%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D(%5Ccdot%2C0)$ ?為例,? 那么有? $u_%7B%5Ccdot%2C0%7D%3Du(%5Ccdot%2C0)$ ?和 $u_%7B%5Ccdot%2C-1%7D%3Du_%7B%5Ccdot%2C0%7D-c%5CDelta%20t%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D(%5Ccdot%2C0)$ .

????????在這個(gè)特定的實(shí)現(xiàn)里,? 當(dāng)滿足? $k%3D%5Cfrac%7Bc%5E2%5CDelta%20t%5E2%7D%7B%5CDelta%20x%5E2%7D%3E1$ ?時(shí),? 會(huì)造成空間數(shù)值不穩(wěn)定,? 值會(huì)以指數(shù)速度發(fā)散;? 當(dāng) $k%3C1$ 時(shí)為時(shí)間數(shù)值不穩(wěn)定,? 時(shí)間不穩(wěn)定不會(huì)造成值發(fā)散,? 但會(huì)產(chǎn)生類似色散的現(xiàn)象;? 只有恰好 $k%3D1$ 時(shí)才是數(shù)值穩(wěn)定的.? 如下圖所示