二維泊松方程求解--點迭代法

1. 問題描述

本算例來自B站Up主“Red-Green鯉魚”的系列教程佩谣。本文主要介紹計算代數(shù)方程組的三種點迭代方法。

1.1. 泊松方程

含有二階偏導(dǎo)數(shù)的偏微分方程：
$\frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial y^2}=f(x,y)$

當 $f=0$ 時，上述方程被稱為拉普拉斯方程窄绒。許多物理過程都可以用泊松方程來描述，如熱傳導(dǎo)方程
$\frac{\partial \phi}{\partial t}=\frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial y^2}$

在求解不可壓縮流動的NS方程時鹃愤，通常將已知壓力場代入動量方程來預(yù)估速度場抄瑟，然后將預(yù)估的速度場代入連續(xù)性方程中，由于預(yù)估的速度場中包含壓力梯度項藏鹊，因此代入連續(xù)性方程后會得到 $\nabla \cdot(\nabla P)$ ，即壓力泊松方程

1.2. 算例

令上述泊松方程的源項為如下形式：
$f(x,y)=-4\cdot sin(x-y)\cdot e^{(x-y)}$

其解析解為
$\phi(x,y)=cos(x-y)\cdot e^{(x-y)}$

比較數(shù)值解與解析解在定義域 $x\in [-1,1], y\in [-1,1]$ 上的差別转锈。邊界條件為Dirichlet盘寡，邊界上的值由解析解求出。

2. 區(qū)域離散和方程離散

$x$ 方向設(shè)置 $N$ 個結(jié)點撮慨，編號從 $1-N$ 竿痰， $y$ 方向上設(shè)置 $M$ 個結(jié)點脆粥，編號從 $1-M$ ，結(jié)點間距分別為 $\Delta x$ 和 $\Delta y$ 影涉，將 $(i,j)$ 號結(jié)點記為 $P$ 变隔，該結(jié)點上下左右四個結(jié)點分別記為 $N,S,W,E$

采用有限差分法計算泊松方程，二階偏導(dǎo)數(shù)項采用二階中心差分離散蟹倾，
$\frac{\partial^2 \phi}{\partial x^2}\bigg|_P=\frac{\phi_W+\phi_E-2\phi_P}{\Delta x^2}$

$\frac{\partial^2 \phi}{\partial y^2}\bigg|_P=\frac{\phi_N+\phi_S-2\phi_P}{\Delta y^2}$

$\frac{\phi_W+\phi_E-2\phi_P}{\Delta x^2}+\frac{\phi_N+\phi_S-2\phi_P}{\Delta y^2}=f_P$

$\phi_P=\frac{(\phi_W+\phi_E)\Delta y^2+(\phi_N+\phi_S)\Delta x^2-f_P \Delta x^2 \Delta y^2}{2(\Delta x^2+\Delta y^2)}$

令 $\Delta x=\Delta y=h$ 匣缘，則上式化為
$\phi_P=\frac{1}{4}(\phi_W+\phi_E+\phi_N+\phi_S-h^2f_P)$

2.1. 邊界條件

左邊界
$\phi_{1,j}=cos(-1-y_j)\cdot e^{(-1-y_j)}$

右邊界
$\phi_{N,j}=cos(1-y_j)\cdot e^{(1-y_j)}$

下邊界
$\phi_{i,1}=cos(x_i+1)\cdot e^{(x_i+1)}$

上邊界
$\phi_{i,N}=cos(x_i-1)\cdot e^{(x_i-1)}$

3. 代數(shù)方程組求解

上一篇文章介紹了代數(shù)方程組求解方法中的直接解法TDMA，本文介紹另一大類求解方法：迭代法喊式。迭代法的思想也可以概況為“預(yù)測-校正”孵户，給出初始值，通過不斷迭代逐步改進岔留，直到達到一定精度要求為止夏哭。
該方法需要首先構(gòu)造迭代方式；其次是所構(gòu)造的迭代序列是否收斂献联，如果收斂則要進一步提高收斂速度竖配。

迭代法可分為點迭代、塊迭代里逆、交替方向迭代法以及強隱迭代法进胯。在點迭代法中，每一步計算只能改進求解區(qū)域中一個結(jié)點的值原押，且該值是由一個顯函數(shù)形式由其余各點的已知值來確定胁镐，因而點迭代法又稱為顯式迭代法。

下面將討論點迭代法的三種實施方式诸衔。

3.1. 雅可比迭代

$\phi_P^{n+1}=\frac{1}{4}(\phi_W^n+\phi_E^n+\phi_N^n+\phi_S^n-h^2f_P)$
上式中上標 $n$ 為當前預(yù)測值盯漂， $n+1$ 為代入迭代方程后的校正值

3.2. 高斯-賽德爾迭代

$\phi_P^{n+1}=\frac{1}{4}(\phi_W^{n+1}+\phi_E^n+\phi_N^n+\phi_S^{n+1}-h^2f_P)$

在逐點計算過程中， $W$ 和 $S$ 點的值在本次迭代過程中已知笨农，因此將已知值代入迭代方程中

3.3. SOR迭代

Successive Over Relaxation就缆，逐次超松弛。SOR迭代法收斂的充要條件是松弛因子 $0<\beta <2$ 谒亦，當 $\beta>1$ 時能夠起到加速收斂的效果竭宰。
$\phi_P^{n+1}=\frac{\beta}{4}(\phi_W^{n+1}+\phi_E^n+\phi_N^n+\phi_S^{n+1}-h^2f_P)+(1-\beta)\phi_P^n$
當 $\beta=1$ ，SOR迭代退化為Gauss-Seidel. 在《Computational Methods for Fluid Dynamics》第四版5.3.3節(jié)給出了一些關(guān)于SOR的討論份招。

上述方程變換形式可得

$\phi^{n+1}=\phi^n+\beta(\phi_{GS}^{n+1}-\phi^n)$
其中 $\phi_{GS}$ 為Gauss-Seidel法求出的值切揭。進一步移項， $\beta$ 用其他幾項表示锁摔，可以得出 $\beta>1$ 能夠加速迭代伴箩，即縮小初始值變化到最終值所需時間。

3.4. 迭代收斂標準

本文使用如下標準來定義收斂標準
$Residual=\frac{1}{N}\sum_{i=1}^{N}\frac{|\phi_i^{n+1}-\phi_i^{n}|}{|\phi_i^{n+1}+\varepsilon|}<10^{-5}$
此處 $N$ 為全局結(jié)點個數(shù)鄙漏，上式表示結(jié)點平均相對偏差小于 $10^{-5}$ 時嗤谚，認為達到收斂。 $\varepsilon$ 為極小值怔蚌，防止分母為0

除此之外巩步，《數(shù)值傳熱學(xué)》還介紹了其他收斂標準。

3.5. 迭代法收斂的分析

《數(shù)值傳熱學(xué)》第二版7.4節(jié)寫道桦踊，對于如下形式的方程：
$a_PT_P=\sum a_{nb}T_{nb} +b$

Jacobi與Gauss-Seidel迭代法收斂的一個充分條件是：系數(shù)矩陣不可約且按行或按列弱對角占優(yōu)椅野。其中“弱對角占優(yōu)”需滿足：
$\frac{\sum |a_{nb}|}{|a_P|}\leq 1 \quad or \quad |a_P| \ge \sum |a_{nb}|$
對各行成立，且其中至少對一行不等號成立籍胯。在本文的算例中竟闪，系數(shù)矩陣的第一行和最后一行對應(yīng)為不等號，其他各行均是等號成立杖狼，即 $|a_P| = \sum |a_{nb}|$

3.6. 上述迭代方法的計算結(jié)果

兩個方向各自的結(jié)點數(shù) $N=101$

Method	Iteration number
Jacobi	9141
GS	5121
SOR, $\beta=1.5$	2065
SOR, $\beta=1.9$	403

Jacobi

jacobi.png

Gauss-Seidel

gs.png

SOR, $\beta=1.5$

sor-1.5.png

SOR, $\beta=1.9$

sor-1.9.png

Analytical solution

analytical.png

4. 代碼

代碼鏈接

最后編輯于：2022.03.29 13:38:54

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者

人面猴
序言：七十年代末炼蛤，一起剝皮案震驚了整個濱河市，隨后出現(xiàn)的幾起案子蝶涩，更是在濱河造成了極大的恐慌理朋，老刑警劉巖，帶你破解...
沈念sama閱讀 218,204評論 6贊 506
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件绿聘，死亡現(xiàn)場離奇詭異嗽上，居然都是意外死亡，警方通過查閱死者的電腦和手機熄攘，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 93,091評論 3贊 395
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進店門兽愤，熙熙樓的掌柜王于貴愁眉苦臉地迎上來，“玉大人挪圾，你說我怎么就攤上這事浅萧。” “怎么了洛史？”我有些...
開封第一講書人閱讀 164,548評論 0贊 354
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵惯殊，是天一觀的道長。經(jīng)常有香客問我也殖，道長土思，這世上最難降的妖魔是什么？我笑而不...
開封第一講書人閱讀 58,657評論 1贊 293
?港島之戀（遺憾婚禮）
正文為了忘掉前任忆嗜，我火速辦了婚禮己儒，結(jié)果婚禮上，老公的妹妹穿的比我還像新娘捆毫。我一直安慰自己闪湾，他們只是感情好，可當我...
茶點故事閱讀 67,689評論 6贊 392
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布绩卤。她就那樣靜靜地躺著途样，像睡著了一般江醇。火紅的嫁衣襯著肌膚如雪。梳的紋絲不亂的頭發(fā)上何暇，一...
開封第一講書人閱讀 51,554評論 1贊 305
城市分裂傳說
那天陶夜，我揣著相機與錄音，去河邊找鬼裆站。笑死条辟，一個胖子當著我的面吹牛，可吹牛的內(nèi)容都是我干的宏胯。我是一名探鬼主播羽嫡，決...
沈念sama閱讀 40,302評論 3贊 418
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼，長吁一口氣：“原來是場噩夢啊……” “哼肩袍！你這毒婦竟也來了杭棵？” 一聲冷哼從身側(cè)響起，我...
開封第一講書人閱讀 39,216評論 0贊 276
萬榮殺人案實錄
序言：老撾萬榮一對情侶失蹤了牛，失蹤者是張志新（化名）和其女友劉穎颜屠，沒想到半個月后，有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體鹰祸，經(jīng)...
沈念sama閱讀 45,661評論 1贊 314
?護林員之死
正文獨居荒郊野嶺守林人離奇死亡甫窟，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點故事閱讀 37,851評論 3贊 336
?白月光啟示錄
正文我和宋清朗相戀三年，在試婚紗的時候發(fā)現(xiàn)自己被綠了蛙婴。大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片粗井。...
茶點故事閱讀 39,977評論 1贊 348
活死人
序言：一個原本活蹦亂跳的男人離奇死亡，死狀恐怖街图，靈堂內(nèi)的尸體忽然破棺而出浇衬，到底是詐尸還是另有隱情，我是刑警寧澤餐济，帶...
沈念sama閱讀 35,697評論 5贊 347
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布耘擂，位于F島的核電站，受9級特大地震影響絮姆，放射性物質(zhì)發(fā)生泄漏醉冤。R本人自食惡果不足惜，卻給世界環(huán)境...
茶點故事閱讀 41,306評論 3贊 330
男人毒藥：我在死后第九天來索命
文/蒙蒙一篙悯、第九天我趴在偏房一處隱蔽的房頂上張望蚁阳。院中可真熱鬧，春花似錦鸽照、人聲如沸螺捐。這莊子的主人今日做“春日...
開封第一講書人閱讀 31,898評論 0贊 22
一樁弒父案，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽定血。三九已至赔癌，卻和暖如春，著一層夾襖步出監(jiān)牢的瞬間澜沟，已是汗流浹背届榄。一陣腳步聲響...
開封第一講書人閱讀 33,019評論 1贊 270
情欲美人皮
我被黑心中介騙來泰國打工，沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留倔喂，地道東北人。一個月前我還...
沈念sama閱讀 48,138評論 3贊 370
代替公主和親
正文我出身青樓靖苇，卻偏偏與公主長得像席噩，于是被迫代替她去往敵國和親。傳聞我的和親對象是個殘疾皇子贤壁，可洞房花燭夜當晚...
茶點故事閱讀 44,927評論 2贊 355