數(shù)值求解波動方程 [4]
這篇里講述了另外一種求得 PML 的方法以及簡單實現(xiàn)二維波動方程的數(shù)值模擬.
另一個 PML 求解
????????在上一篇求一維 PML 時,? PML 變換為?,? 以及在假設(shè)
?時得到關(guān)系
,? 并且在最后恰好全部虛數(shù)單位
?都被抵消,? 留下一個純實數(shù) (復(fù)數(shù)也可以) 的方程.? 在部分情況下,? 需要施加 PML 變換的波動方程比較復(fù)雜 (比如強(qiáng)迫振動或粘度等),? 造成虛數(shù)單位無法恰好被消除,? 留下一個需要在復(fù)數(shù)域上求解的方程.? 也不是說數(shù)值求解復(fù)數(shù)方程不可以,? 只不過考慮到內(nèi)存和計算開銷,? 復(fù)數(shù)真的不太行.
????????實際上,? 考慮傅里葉變換?,? 無論
?如何都有
,? 這與上面提到的關(guān)系是等價的.? 另外再考慮拉普拉斯變換
,? 除了積分域和定義域不一樣外,? 它跟傅里葉變換是相同的.? 在波動方程里,? 不需要太關(guān)心
?部分 (這里只討論正演而忽略反演),? 那么把
?替換掉所有出現(xiàn)的頻率得到新的 PML 變換
?和關(guān)系
.? 下面繼續(xù)以一維波動方程做例子:
????????對一維波動方程??施加拉普拉斯變換得到
,? 即
?(偏微分算符和拉普拉斯變換的線性),? 這個就是波動方程在 s 域的形式.? 對 s 域的波動方程施加 PML 變換得?
,? 兩邊同乘
?得
.? 引入新量
?使得
,? 即得到
,? 把方程組從 s 域轉(zhuǎn)為時域得
,? 可以看到這與之前得到的一模一樣.
????????盡管兩個方法得到的結(jié)果都是一樣的,? 但在 s 域上處理 PML 變換是非常通用的,? 至少不需要假設(shè)一些條件.
二維波動方程
????????二維波動方程為?,? 其中
?為所求波的解.? 習(xí)慣上使用
作為二維笛卡爾坐標(biāo),? 但是 unicode 沒有下標(biāo) y,? 所以這里使用?
?作為坐標(biāo)了 (方便以后的編程實現(xiàn)).? 跟一維波動方程類似,? 二維波動方程離散化為?
有一說一, 確實有點長了.? 整理后得到
其中?.? 在離散化后還需要截斷計算區(qū)域,? 截斷方式與一維的一樣,? 但只不過現(xiàn)在不管是截斷一個維度,? 而是兩個.? 但是值得留意的是,? 邊界圍成的形狀不一定只能是正方形長方形,? 但奇形怪狀的邊界很多時候需要額外計算邊界處的離散化,? 所以這里只展示正方形的邊界了.
二維實現(xiàn)
????????在模擬二維波動方程里,? 有一個比較大的問題就是計算速度,? 因為網(wǎng)格大小會隨著空間精度和大小的增大呈平方上升,? 所以與其放在 cpu 上面跑,? 不如放在 gpu 上面.
????????julia 的 cuda 支持在?https://cuda.juliagpu.org/stable?有詳細(xì)記述.? 一般來說有顯卡驅(qū)動時直接在 REPL 里 `]add CUDA` 就可以安裝了,? 并且 julia 函數(shù)可以編譯為 cuda 代碼,? 而不用像 python 那樣在字符串里寫 C++ CUDA 函數(shù).? 下面是 Wave2D 相關(guān)的實現(xiàn),? 在代碼里上下邊界使用了 Neumann 邊界,? 而左右邊界則是循環(huán)邊界.
????????不過在渲染 gif 的時候無論如何都不能把文件大小降下來,? 明明隔壁 mp4 都超級小了,? 怎么 gif 就那么大 (mp4: 211K, gif: 13.8M).? 所以下面貼圖片算了,? 但是代碼是直接渲染視頻的.
????????觀察到圖中,? 波包在空間擴(kuò)散時,? 在后面留下了不為零的"尾跡",? 并且尾跡不會消失 (可以參考《數(shù)理方程》)?,? 亦即不符合 Huggens 原理,? 這種現(xiàn)象稱為波的彌散或有后效現(xiàn)象,? 并且這是偶數(shù)維度空間特有的.? 下圖是沿軸方向的切片,? 可以更直觀地看到這個現(xiàn)象:
???????在上面的代碼實現(xiàn)里,? 已經(jīng)包括了不均勻介質(zhì),? 下圖為空間中心有一個半徑為 0.1 的實心玻璃圓.
????????渲染的代碼稍微有點長,? 所以這里單獨拿出來貼一份:
數(shù)值穩(wěn)定性
????????因為二維波動方程有三個維度 (兩個空間一個時間),? 那么應(yīng)該討論三個維度的數(shù)值穩(wěn)定性.? 據(jù)自己的觀察經(jīng)驗,? 當(dāng)??時,?
?會造成空間數(shù)值不穩(wěn)定,? 否則
是空間數(shù)值不穩(wěn)定 (? 我也不知道為何).? 在空間維度上,??
?是衡量時間不穩(wěn)定的指標(biāo),? 與一維波動方程類似,? 這兩個數(shù)值在?
?時會在相應(yīng)維度上產(chǎn)生時間不穩(wěn)定,? 但必須確保全局上不能有空間不穩(wěn)定,? 即
,? 所以任何維度都不能不產(chǎn)生時間不穩(wěn)定.
????????盡管在任何維度上必定是時間不穩(wěn)定的,? 但是當(dāng)??時,? 在空間對角線上是數(shù)值穩(wěn)定的,? 如下圖所示:
????????為了應(yīng)對這種情況,? 可以選擇具有復(fù)雜網(wǎng)格的實現(xiàn)方式,? 比如旋轉(zhuǎn)網(wǎng)格或自適應(yīng)網(wǎng)格等.? 但這些技術(shù)的輪子實在不太好搓,? 這里就 pass 了.
