PDE有限差分方法(8)

您所在的位置：网站首页 › 积分的周期性公式推导过程图 › PDE有限差分方法(8)

PDE有限差分方法(8)

2024-05-21 05:09| 来源: 网络整理| 查看: 265

今天回顾一维扩散方程（系数充分光滑）的差分格式怎么构造，以及稳定性的分析方法. 这一章的重点是格式的构造方式与变系数格式稳定性分析方法, 关于相容性的分析, 在前面应该很熟悉了, 自证.

参考书：(1) J.W. Thomas - Numerical Partial Differential Equations_ Finite Difference Methods (1995, Springer)(2) 张强《偏微分方程的有限差分方法》科学出版社，2019年1月版。(3) K. W. Morton, D. F. Mayers - Numerical solution of partial differential equations (2005, Cambridge University Press)

今天的内容包括：

下面考虑非均匀介质的热传导方程, 有如下两种形式:

(i)非守恒型扩散方程

$u_t=a(x,t)u_{xx} \\$

(ii)守恒型扩散方程

$u_t=(a(x,t)u_x)_x. \\$

其中, a(x,t) 是扩散系数. 为了保证方程是抛物型的, 这里设 $\inf\limits_{\forall x\forall t}a(x,t)\geq \delta_00,$ 即扩散系数在计算区域内有正的下确界. 如果给出一定的初始条件与边界条件, 可以证明真解存在唯一(回顾PDE). 下面都考虑纯初值问题或者周期边值问题, 并设 a(x,t),u(x,t) 都足够光滑, 时空网格 $T_{\Delta x,\Delta t}$ 等距, 空间步长和时间步长分别为 $\Delta x,\Delta t,$ 网比为 $\mu=\dfrac{\Delta t}{(\Delta x)^2}.$

1 非守恒型扩散方程

用局部冻结技术可以构造非守恒型扩散方程的差分方程.

在离散焦点冻结扩散系数, 可以得到非守恒型扩散方程的全显格式:

$u_j^{n+1}-u_j^n=\mu a_j^n\delta_x^2u_j^n,$

同理有全隐格式

$u_j^{n+1}-u_j^n=\mu a_j^n\delta_x^2u_j^{n+1}.$

这两个格式无条件具有(2,1)阶局部截断误差.

加权平均格式的扩散系数有两种冻结策略:

(1) 多焦点策略

$u_j^{n+1}-u_j^n=\mu\left[\theta a_j^{n+1}\delta_x^2a_j^{n+1}+(1-\theta)a_j^n\delta_x^2u_j^n\right]$

其中 $\theta\in[0,1],$ 且当 $\theta\neq \dfrac{1}{2}$ 时, 它无条件有(2,1)阶局部截断误差. 当 $\theta=\dfrac{1}{2}$ 时, 这是CN格式, 无条件具有(2,2)阶局部截断误差.

(2) 单焦点策略

$u_j^{n+1}-u_j^n=\mu a_j^*\left[\theta \delta_x^2a_j^{n+1}+(1-\theta)\delta_x^2u_j^n\right]$

这里 a_j^* 是扩散系数的局部冻结.

当 $\theta\neq\dfrac{1}{2}时, 可以令a_j^*=a(x_j,t^*),$ 其中 $t^*\in[t^n,t^{n+1}],$ 此时差分方程无条件有(2,1)阶局部截断误差.

当 $\theta=\dfrac{1}{2}$ 时, 这是CN格式. 为了保持(2,2)阶局部截断误差, 需要精心设置 a_j^*, 设置方法有两种:

(a)仿照多焦点策略, 定义

$a_j^*=\dfrac{1}{2}(a_j^n+a_j^{n+1});$

(b)仿照线性常系数CN格式的构造过程, 直接把扩散系数冻结在最佳的离散焦点上:

$a_j^*=a_j^{n+\frac{1}{2}}\equiv a\left(x_j,\dfrac{t^n+t^{n+1}}{2}\right).$

可以发现, (a)是基于时间积分的梯形公式, (b)是基于时间积分的中点矩形公式. 容易证明这两种冻结方式都有(2,2)阶局部截断误差.

注：上述所有格式都与常系数格式的相容阶保持一致.

下面给出相容阶更高的差分格式. 为了简单, 设扩散系数与时间无关, 即 $a(x,t)\equiv a(x),$ 基于加权平均格式, 构造非守恒型扩散方程的整体四阶格式(Douglas格式, 有条件相容).

回顾CN格式, 根据构造过程可知

$\dfrac{[u]_j^{n+1}-[u]_j^n}{a_j\Delta t}- \dfrac{1}{2(\Delta x)^2}\left[\delta_x^2[u]_j^n+\delta_x^2[u]_j^{n+1}\right] =-\dfrac{(\Delta x)^2}{12}[u_{xxxx}]_j^{n+\frac{1}{2}}+O((\Delta x)^4+(\Delta t)^2).$

只需建立 $[u_{xxxx}]_j^{n+\frac{1}{2}}$ 的二阶相容的离散.

把空间导数变成时间导数, 并用相应的中心差商离散, 可得

$\begin{aligned} \, [u_{xxxx}]_j^{n+\frac{1}{2}}&=[(a^{-1}u_t)_{xx}]_j^{n+\frac{1}{2}} =[(a^{-1}u)_{xxt}]_j^{n+\frac{1}{2}} \\ &=\dfrac{\delta_x^2[a^{-1}u]_j^{n+1}-\delta_x^2[a^{-1}u]_j^n}{(\Delta x)^2\Delta t} +O((\Delta x)^2+(\Delta t)^2). \end{aligned}$

联立两式并用数值解代替真解, 可得整体四阶的差分方程:

$\dfrac{\Delta_tu_{j+1}^n}{12a_{j+1}}+\dfrac{5\Delta_tu_j^n}{6a_j^n} +\dfrac{\Delta_tu_{j-1}^n}{12a_{j-1}}=\dfrac{1}{2}\mu(\delta_x^2u_j^{n+1}+\delta_x^2u_j^n).$

若 $a(x)\equiv a,$ 此时它就是常系数的Douglas格式.

注1：随着相容阶升高, 局部冻结会变得更加繁琐.

注2：以偏微分方程作为桥梁, 把时空方向导数进行相互转换, 克服了时空信息分布不均的困难,使得离散模板的空间网格点分布更加紧凑, 从而可以把低阶格式转化为高阶格式.

2 守恒型扩散方程

最容易想到的差分格式设计思想是直接展开然后用差商离散导数.

$u_t=a(x,t)u_{xx}+b(x,t)u_x,\text{其中}b=a_x.$

用中心差商离散空间导数, 向前差商离散时间导数, 并采用冻结系数方法, 可得全显格式

$u_j^{n+1}-u_j^n=\mu a_j^n\delta_x^2u_j^n+\dfrac{\Delta x}{2}\mu b_j^n(u_{j+1}^n-u_{j-1}^n).$

它无条件具有(2,1)阶局部阶段误差. 但是当 a(x,t) 变化剧烈(或者是间断)的时候, 稳定性缺陷很严重. 这是因为, 守恒型扩散方程的热量守恒性质

$\int_{z_1}^{z_2}u_t(x,t)dx=W(z_2,t)-W(z_1,t), \forall z_1,z_2,\forall t,$

没有得到满意的数值保持, (即流入、流出不“相容”). 其中,

$W=au_x \\$

称为热流通量. 本章范围内的都采取这个定义保持不变.

我们需要数值实现热量的局部守恒性质. 常见方法如下:

(1) 积分插值方法: 适合用于散度型导数的离散. 设计思想: 离散对象在某个局部区域的积分近似. 步骤如下:

选取适当的积分区域, 对PDE作积分.用散度定理把高阶导数的高维积分变成低阶导数的低维积分.用恰当的数值积分公式处理所得积分恒等式.离散低阶导数.

下面用积分插值方法建立守恒型扩散方程的全显格式.

(a)首先选取区域 $\left(x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}\right)\times (t^n,t^{n+1}),$ 对PDE作积分, 然后根据Newton-Leibniz公式可以得到如下恒等式:

$\int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}u(x,t^{n+1})dx -\int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}u(x,t^{n})dx =\int_{t^n}^{t^{n+1}}W(x_{j+\frac{1}{2}},t)dt-\int_{t^n}^{t^{n+1}}W(x_{j-\frac{1}{2}},t)dt$

(b)左侧采用中点公式近似, 右边采用左矩形公式近似, 可得

$\left([u]_j^{n+1}-[u]_j^n\right)\Delta x\approx \left([W]_{j+\frac{1}{2}}^n- [W]_{j-\frac{1}{2}}^n\right)\Delta t$

(c)用一阶中心差商技术与冻结系数方法来离散热流通量可得

$[W]_{j+\frac{1}{2}}^n\approx a_{j+\frac{1}{2}}^n\dfrac{[u]_{j+1}^n-[u]_j^n}{\Delta x},$

其中 $a_{j+\frac{1}{2}}^n$ 是扩散系数的局部冻结, 冻结方式可以是

$a_{j+\frac{1}{2}}^n=a(x_{j+\frac{1}{2}},t^n)\text{或}a_{j+\frac{1}{2}}^n =\dfrac{1}{2}(a_j^n+a_{j+1}^n).$

略去无穷小量, 用数值解代替真解, 可得守恒型扩散方程的全显格式:

$u_j^{n+1}-u_j^n=\mu\left[a_{j+\frac{1}{2}}^n(u_{j+1}^n-u_j^n) -a_{j-\frac{1}{2}}^n(u_j^n-u_{j-1}^n)\right].$

或者写

$u_j^{n+1}-u_j^n=\mu\delta_x(a_j^n\delta_xu_j^n)$

它无条件具有(2,1)阶局部截断误差.

注1：如果把步骤(b)的右边改为用右矩形公式近似, 可得全隐格式.

注2：右侧积分在区间 $[t^n,t^n+\theta(t^{n+1}-t^n)]$ 用左矩形公式, 在区间 $[t^n+\theta(t^{n+1}-t^n),t^{n+1}]$ 用右矩形公式, 可得加权平均格式

$u_j^{n+1}=u_j^n+\theta\mu\delta_x(a_j^n\delta_xu_j^n) +(1-\theta)\mu\delta_x(a_j^{n+1}\delta_xu_j^{n+1}).$

局部截断误差阶与常系数的情形保持一致.

(2) 盒子格式: 基本思想是对微分方程进行降阶处理. 离散对象不是热传导方程, 而是与其等价的一阶微分方程组

$\boxed{u_t=v_x, \qquad v=a(x,t)u_x.}$

这里是原始变量, 是辅助变量. 理论上辅助变量可以任意设置, 但它通常具有一定的计算目标或物理意义. 在这里, 辅助变量就是热流通量, 即 W=v=au_x.

具体步骤如下:

(a) 在时空网格中取四个紧凑排列的网格点, 搭建出盒子形状的离散模板, 离散焦点为 $(x_{j+\frac{1}{2}},t^{n+\frac{1}{2}}).$

(b) 用中心差商技术离散第一个方程, 可得

$\dfrac{[u]_{j+\frac{1}{2}}^{n+1}-[u]_{j+\frac{1}{2}}^n}{\Delta t} \approx \dfrac{[v]_{j+1}^{n+\frac{1}{2}}-[v]_j^{n+\frac{1}{2}}}{\Delta x},$

$[v]_{j+\frac{1}{2}}^n\approx a_{j+\frac{1}{2}}^n\dfrac{[u]_{j+1}^n-[u]_j^n}{\Delta x}.$

这里 $a_{j+\frac{1}{2}}^n$ 的冻结方式如前.

(d) 用整点的真解来逼近各边中点的真解, 即

$\begin{aligned} \,[w]_{j+\frac{1}{2}}^n\approx\dfrac{1}{2}\left([w]_j^n+[w]_{j+1}^n\right) =\Pi_x[w]_{j+\frac{1}{2}}^n, \\ [w]_j^{n+\frac{1}{2}}\approx\dfrac{1}{2}\left([w]_j^n+[w]_j^{n+1}\right) =\Pi_t[w]_{j}^{n+\frac{1}{2}}, \end{aligned}$

这里 $w=u\text{或}v, \Pi$ 是算术平均算子.

(e)略去无穷小量, 用数值解代替真解, 可得守恒型扩散方程的盒子格式:

$\begin{aligned} &\dfrac{\Pi_xu_{j+\frac{1}{2}}^{n+1}-\Pi_xu_{j+\frac{1}{2}}^n}{\Delta t} =\dfrac{\Pi_tv_{j+1}^{n+\frac{1}{2}}-\Pi_tv_j^{n+\frac{1}{2}}}{\Delta x}, \\ &\Pi_xv_{j+\frac{1}{2}}^n=a_{j+\frac{1}{2}}^n\dfrac{u_{j+1}^n-u_j^n}{\Delta x}. \end{aligned}$

这个格式无条件具有(2,2)阶局部截断误差(分析相容阶的时候需要消去v).

如果扩散系数恒定, 即 $a(x,t)\equiv a0,$ 则利用Fourier方法可以知道盒子格式无条件具有 L^2 模稳定性.

不是常系数的差分格式不能用Fourier方法! 这就是为什么我们前面一直都没有提到“稳定性”这个概念. 接下来我们来看一下上述变系数差分格式的稳定性分析方法.

3 稳定性分析方法

Lax-Richtmyer等价定理在变系数情况下依旧成立, 即相容差分格式的稳定性与收敛性等价. 所以重点讨论对象依然是相容性与稳定性. 相容性分析方法依然是Taylor展开, 但稳定性的分析方法不能用Fourier方法, 要采取其他方式.

$\mathbf{1\quad 冻结系数方法}$

出发点: 当离散网格变密的时候, 扩散系数可以视为局部恒定, 甚至可以把线性变系数差分格式看作某个线性(或分片)常系数差分格式的微小扰动. 步骤:

把差分系数冻结为常数, 导出相应差分格式.用其他准确的分析技术(Fourier方法、最大模原理等等)给出相应的稳定性结论.考虑所有合理的系数冻结范围, 所有稳定性结论的交集就是冻结系数法给出的结果.例1 用冻结系数法给出非守恒型扩散方程的全显格式的最大模与 L^2

模稳定性条件: $u_j^{n+1}-u_j^n=\mu a_j^n\delta_x^2u_j^n.$

把 a_j^n 冻结为某个常数 , 则全显格式可以变成线性常系数差分格式

$u_j^{n+1}-u_j^n=\mu a\delta_x^2u_j^n.$

根据已知的稳定性结果可知, 最大模稳定性条件与 L^2 模稳定性条件都是 $\mu a \leq\dfrac{1}{2},$

但这个结论只能反映网格点 (x_j,t^n) 附近的情形, 需要综合考虑所有的网格点. 让遍历 a_j^n 的所有可能取值, 相应稳定性结论的交集是

$\max\limits_{\forall x\forall t}a(x,t)\mu\leq\dfrac{1}{2}.$

这就是非守恒型扩散方程的全显格式的两种稳定性条件. \QED

注1：根据离散最大模原理, 这个条件也可以保证全显格式的最大模稳定性. 但是在临界状态下, 全显格式的 L^2 模稳定性无法得到严格的理论验证

注2： 冻结系数法完全忽略了系数变化带来的数值影响, 对于线性常系数问题稳定的某些格式, 有可能对于线性变系数问题出现“数值共振”现象(不满足最大模原理), 使得部分简谐波出现无法控制的增长, 造成线性不稳定现象.

类似可以分析守恒型扩散方程的最大模稳定性条件是

$\max\limits_{\forall x\forall t}a(x,t)\mu\leq\dfrac{1}{2}, \text{且还要满足}\Delta x\leq\inf\limits_{\forall x\forall t}\dfrac{2|a(x,t)|}{|b(x,t)|}.$

若扩散系数 a(x,t) 沿着空间的变化非常剧烈( b(x,t) 很大), 如

$a(x,t)=\sin^2\dfrac{x}{\varepsilon}+1, \varepsilon1,$

则 $\Delta x$ 要足够小, 空间网格要足够密集, 计算量很大, 效率很低. 一般来说, 若导数具有紧凑的散度型结构, 相应的数值格式更具有优势.

$\mathbf{2\quad 能量方法}$

能量方法可以用在线性变系数、非周期边值问题、非等距网格等情况. 过程与PDE里面的能量方法类似. 分析步骤:

选取适当的检验函数, 建立能量范数的递推关系式.指出能量范数与离散 L^2

模的等价关系.给出

模稳定性的充分条件.例2 考虑守恒型扩散方程 u_t=(au_x)_x

的周期边值问题, 扩散系数 a(x,t),u(x,t)

均以1为周期. 用偏微分方程的能量方法可知 $\int_0^1u^2dx\leq \int_0^1u_0^2dx,\text{其中}u_0\text{是初值函数}$ 基于等距时空网格构造周期边值问题的全显格式 $u_j^{n+1}-u_j^n=\mu\left[a_{j+\frac{1}{2}}^n(u_{j+1}^n-u_j^n) -a_{j-\frac{1}{2}}^n(u_j^n-u_{j-1}^n)\right],$ 用能量方法给出 L^2

模稳定的充分条件.

在差分方程两端同乘检验函数 $u_j^{n+1}+u_j^n, j=0:J-1,$ 把个恒等式相加可得

$\begin{aligned} LHS&\triangleq \sum\limits_{j=0}^{J-1}(u_j^{n+1}-u_j^n)(u_j^{n+1}+u_j^n)\Delta x \\ &=\mu\sum\limits_{j=0}^{J-1}\left[a_{j+\frac{1}{2}}^n(u_{j+1}^n-u_j^n) -a_{j-\frac{1}{2}}^n(u_j^n-u_{j-1}^n)\right](u_j^{n+1}+u_j^n)\Delta x\triangleq RHS \end{aligned}$

上式中

$LHS=\sum\limits_{j=0}^{J-1}(u_j^{n+1})^2\Delta x-\sum\limits_{j=0}^{J-1}(u_j^n)^2\Delta x =\|u^{n+1}\|_{2,\Delta x}^2-\|u^n\|_{2,\Delta x}^2.$

为了简单, 我们设 a(x,t)=a(x), 即扩散系数与时间无关(扩散系数与时间有关的时候分析过程是类似的).

为了简洁, 记 $\theta_{j+\frac{1}{2}}^n=a_{j+\frac{1}{2}}^n(u_{j+1}^n-u_j^n), v_j^n=u_j^{n+1}+u_j^n,$ 根据周期边界条件 $u_0^n=u_J^n,u_{J+1}^n=u_1^n$ 以及 a(x) 周期性, 利用分部求和, 可知

$\begin{aligned} RHS&=\sum\limits_{j=0}^{J-1}(\theta_{j+\frac{1}{2}}^n-\theta_{j-\frac{1}{2}}^n) (u_j^{n+1}+u_j^n)\Delta x \\ &=-\mu\sum\limits_{j=0}^{J-1}\theta_{j-\frac{1}{2}}^n(v_j^n-v_{j-1}^n) -\mu\theta_{J-\frac{1}{2}}^nv_j^n+\mu\theta_{-\frac{1}{2}}^nv_0^n \\ &=-\mu\sum\limits_{j=0}^{J-1}\theta_{j-\frac{1}{2}}^n(v_j^n-v_{j-1}^n) \\ &=-\mu\sum\limits_{j=0}^{J-1}a_{j-\frac{1}{2}}^n(u_j^n-u_{j-1}^n) (u_{j}^{n+1}-u_{j-1}^{n+1}+u_j^n-u_{j-1}^n). \end{aligned}$

注意到恒等式 $p(p+q)=\dfrac{1}{2}(p+q)^2+\dfrac{1}{2}(p^2-q^2),$ 这里 $p=u_j^n-u_{j-1}^n,q=u_{j-1}^{n+1}-u_{j}^{n+1},$ 则

$RHS=-\mu\sum\limits_{j=0}^{J-1}a_{j-\frac{1}{2}}\left(\dfrac{1}{2}(p+q)^2 +\dfrac{1}{2}(p^2-q^2)\right).$

定义系统的能量范数

$\mathcal{E}(u^n)=\|u^n\|_{2,\Delta x}^2-\dfrac{\mu}{2} \sum\limits_{j=0}^{J-1}a_{j-\frac{1}{2}}(u_j^n-u_{j-1}^n)^2\Delta x,$

则 $\mathcal{E}(u^{n+1})-\mathcal{E}(u^n) =-\dfrac{1}{2}\mu\sum\limits_{j=0}^{J-1}a_{j-\frac{1}{2}}(p+q)^2\leq 0.$

则 $\mathcal{E}(u^n)$ 是不增的. 由算术平均值不等式,

$\|u^n\|_{2,\Delta x}^2-\dfrac{1}{2}\mu A(2\|u^n\|_{2,\Delta x}^2+2\|u^n\|_{2,\Delta x}^2) \leq\mathcal{E}(u^n)\leq \mathcal{E}(u^0)\leq\|u^0\|_{2,\Delta x}^2.$

其中 $A=\sup\limits_{x\in(0,1)} a(x)0.$ 若存在正常数 $\delta0$ 使得 $1-2\mu A\geq\delta,$ 则 L^2 模稳定性成立, 即

$\|u^n\|_{2,\Delta x}^2\leq \dfrac{1}{\delta}\|u^0\|_{2,\Delta x}^2.$

易知 $\delta1,$ 则稳定性结论弱于偏微分方程定解问题的适定性结论. \QED

注：在临界状态即 $2\mu A=1$ 下, 变系数全显格式的 L^2 模稳定性结论不明确. 如果冻结系数与时间、空间都有关, 且 $|a_t|\leq C,$ 此时的稳定性结论需要用到Gronwall不等式, 用归纳法可以证明.

定理 [Gronwall不等式]设 $\{f_n\}_{n\geq 0},\{g_n\}_{n\geq 0}$ 是非负序列, $\{g_n\}_{n\geq 0}$ 是递增序列, 如果存在

使得 $f_{n+1}\leq C\sum\limits_{m=0}^{n}f_m\Delta t+g_{n+1},\forall n\geq 0,$ 则当 $\Delta t$ 充分小时, $f_n\leq e^{Cn\Delta t}g_n.$

定义

$\mathcal{E}(u^n)=\|u^n\|_{2,\Delta x}^2-\dfrac{\mu}{2} \sum\limits_{j=0}^{J-1}a_{j-\frac{1}{2}}^n(u_j^n-u_{j-1}^n)^2\Delta x,$

此时可以化出

$\mathcal{E}(u^{n+1})-\mathcal{E}(u^n) +\dfrac{\mu}{2}\sum\limits_{j=0}^{J-1}(a_{j-\frac{1}{2}}^{n+1}-a_{j-\frac{1}{2}}^n)(u_{j}^{n+1} -u_{j-1}^{n+1})^2\Delta x\leq 0.$

由微分中值定理, $|a_{j-\frac{1}{2}}^{n+1}-a_{j-\frac{1}{2}}^n|\leq C\Delta t,$ 则

$\mathcal{E}(u^{n+1})-\mathcal{E}(u^n) \leq \dfrac{\mu C\Delta t}{2}\sum\limits_{j=0}^{J-1}(u_{j}^{n+1} -u_{j-1}^{n+1})^2\Delta x \leq 2C\mu\Delta t\|u^{n+1}\|_{2,\Delta x}^2,$

所以

$\mathcal{E}(u^{n+1})-\mathcal{E}(u^0)\leq 2C\mu \Delta t\sum\limits_{k=1}^{n+1}\|u^k\|_{2,\Delta x}^2.$

由于

$(1-2\mu A)\|u^n\|_{2,\Delta x}^2\leq \mathcal{E}(u^n)\leq \|u^n\|_{2,\Delta x}^2,$

则

$(1-2\mu A)\|u^{n+1}\|_{2,\Delta x}^2\leq \|u^0\|_{2,\Delta x}^2 +2C\mu \Delta t\sum\limits_{k=1}^{n+1}\|u^k\|_{2,\Delta x}^2.$

当 $1-2\mu A-2C\mu\Delta t\geq \delta0$ 时，由Gronwall不等式可知

$\|u^n\|_{2,\Delta x}^2\leq e^{\frac{2C\mu n\Delta t}{1-2\mu A-2C\mu\Delta t}} (1-2C\mu\Delta t)\|u^0\|_{2,\Delta x}^2 \leq e^{2CT\delta^{-1}}\|u^0\|_{2,\Delta x}^2,$

显然, 稳定性结论更弱.

【本文地址】

公司简介

联系我们