二维泊松方程求解--点迭代法

1. 问题描述

本算例来自B站Up主“Red-Green鲤鱼”的系列教程。本文主要介绍计算代数方程组的三种点迭代方法。

1.1. 泊松方程

含有二阶偏导数的偏微分方程:
\frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial y^2}=f(x,y)

f=0时,上述方程被称为拉普拉斯方程。许多物理过程都可以用泊松方程来描述,如热传导方程
\frac{\partial \phi}{\partial t}=\frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial y^2}

在求解不可压缩流动的NS方程时,通常将已知压力场代入动量方程来预估速度场,然后将预估的速度场代入连续性方程中,由于预估的速度场中包含压力梯度项,因此代入连续性方程后会得到\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)}

比较数值解与解析解在定义域x\in [-1,1], y\in [-1,1]上的差别。边界条件为Dirichlet,边界上的值由解析解求出。

2. 区域离散和方程离散

x方向设置N个结点,编号从1-Ny方向上设置M个结点,编号从1-M,结点间距分别为\Delta x\Delta y,将(i,j)号结点记为P,该结点上下左右四个结点分别记为N,S,W,E

采用有限差分法计算泊松方程,二阶偏导数项采用二阶中心差分离散,
\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. 代数方程组求解

上一篇文章介绍了代数方程组求解方法中的直接解法TDMA,本文介绍另一大类求解方法:迭代法。迭代法的思想也可以概况为“预测-校正”,给出初始值,通过不断迭代逐步改进,直到达到一定精度要求为止。
该方法需要首先构造迭代方式;其次是所构造的迭代序列是否收敛,如果收敛则要进一步提高收敛速度。

迭代法可分为点迭代、块迭代、交替方向迭代法以及强隐迭代法。在点迭代法中,每一步计算只能改进求解区域中一个结点的值,且该值是由一个显函数形式由其余各点的已知值来确定,因而点迭代法又称为显式迭代法。

下面将讨论点迭代法的三种实施方式。

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为当前预测值,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)

在逐点计算过程中,WS点的值在本次迭代过程中已知,因此将已知值代入迭代方程中

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节给出了一些关于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为全局结点个数,上式表示结点平均相对偏差小于10^{-5}时,认为达到收敛。\varepsilon为极小值,防止分母为0

除此之外,《数值传热学》还介绍了其他收敛标准。

3.5. 迭代法收敛的分析

《数值传热学》第二版7.4节写道,对于如下形式的方程:
a_PT_P=\sum a_{nb}T_{nb} +b

Jacobi与Gauss-Seidel迭代法收敛的一个充分条件是:系数矩阵不可约且按行或按列弱对角占优。其中“弱对角占优”需满足:
\frac{\sum |a_{nb}|}{|a_P|}\leq 1 \quad or \quad |a_P| \ge \sum |a_{nb}|
对各行成立,且其中至少对一行不等号成立。在本文的算例中,系数矩阵的第一行和最后一行对应为不等号,其他各行均是等号成立,即|a_P| = \sum |a_{nb}|

3.6. 上述迭代方法的计算结果

两个方向各自的结点数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. 代码

代码链接

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 219,110评论 6 508
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 93,443评论 3 395
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 165,474评论 0 356
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,881评论 1 295
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,902评论 6 392
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,698评论 1 305
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,418评论 3 419
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,332评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,796评论 1 316
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,968评论 3 337
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 40,110评论 1 351
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,792评论 5 346
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,455评论 3 331
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 32,003评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 33,130评论 1 272
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 48,348评论 3 373
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 45,047评论 2 355

推荐阅读更多精彩内容