二维泊松方程非均匀网格上的高精度紧致差分格式

二维波动方程的高精度交替方向隐式方法马月珍;李小纲;葛永斌【摘要】基于二阶微商的四阶紧致差商逼近公式及加权平均思想,提出了数值求解二维波动方程的2种精度分别为O(τ~2+h~4)和O(τ~4+h~4)的交替方向隐式(ADI)格式,以及与其相匹配的第一个时间层的同阶离散格式,并且通过Fourier方法分析了格式的稳定性.该方法在沿每个空间方向上只涉及3个网格基架点,因此可以重复采用TDMA算法,从而大大节省计算时间.数值实验验证了所用方法的精确性和可靠性.【期刊名称】《四川师范大学学报（自然科学版）》【年(卷),期】2010(033)002【总页数】5页(P179-183)【关键词】波动方程;高阶紧致格式;交替方向隐式方法;稳定性【作者】马月珍;李小纲;葛永斌【作者单位】宁夏大学,应用数学与力学研究所,宁夏,银川,750021;宁夏大学,应用数学与力学研究所,宁夏,银川,750021;宁夏大学,应用数学与力学研究所,宁夏,银川,750021【正文语种】中文【中图分类】O241.82有限差分法[1-3]是数值求解偏微分方程的常用方法之一.在航空、气象、海洋、水利等许多流体力学的问题中,常常遇到双曲型偏微分方程,针对传统的差分离散格式,普遍有着精度低且受到很强的稳定性条件限制的缺陷,因此发展其高精度且稳定性好的差分离散方法具有十分重要的意义[4-7].交替方向隐式(AD I)方法是数值求解该类问题的一种非常有效的方法,它将高维问题转化为若干一维问题进行求解,而一维问题又可采用高效的 TDMA算法,从而可以大大提高计算效率,节省存储空间.传统的AD I方法如 Peaceman-Rachford(P-R)格式和Beam-War ming[8]格式都是二阶精度的.另一方面,文献[9]提出了求解二维非定常对流扩散方程的高精度AD I格式,文献[10]提出了求解高维热传导方程的高精度AD I格式.本文依然根据AD I方法的基本思想,提出两种求解二维波动方程的高精度紧致AD I差分格式,为此考虑初边值问题：其中Ω ={(x,y):0≤x,y≤1},Γ为Ω的边界,u(x, y,t)为待求未知量,f(x,y,t)为源项,φ(x,y),ψ(x,y)和 g(x,y,t)均为已知函数且具有充分的光滑性.1 高阶紧致AD I格式用τ表示时间步长,空间取等间距网格,步长用h表示.网格点为(xi,yj,tn),xi=ih,yj=jh,tn = nτ,i,j=0,1,…,N,h=1/N,n≥0.1.1 AD I(2,4)格式考虑(1)式在 n时刻值,对时间导数项采用中心差分,空间导数项采用Kreiss[11]提出的四阶紧致差分公式:则有其中在空间方向以和的算术平均值代替可得对(5)式进行整理且略去高阶项可得为了构造AD I差分格式,采用与文献[9-10]类似的技巧,在(6)式左端加上可得显然,(7)式与(6)式的截断误差同阶,利用可将(7)式写为如下形式引入一个过渡变量则可将 (8)式写为对于过渡变量的边界条件,可以由下式给出(9)式即为求解二维波动方程的高精度紧致 AD I格式,其精度为O(τ2+h4),记为AD I(2,4).1.2 AD I(4,4)格式对时间和空间导数项均采用四阶紧致差分公式,考虑(1)式在 n时刻值,可得对上式进行整理且略去高阶项可得可将(10)式写为如下形式为了构造AD I差分格式,采用与文献[9-10]类似的技巧,在 (11)式左端加上并进行因式分解,可得显然,(11)式与(12)式的截断误差同阶,引入一个过渡变量则 (12)式可写为对于过渡变量的边界条件,可以由下式给出(13)式即为求解二维波动方程的高精度紧致 AD I格式,其精度为O(τ4+h4),记为AD I(4,4).1.3 初始条件的离散因为格式是 3层的.即每一次时间推进都需要知道前两个时间步的值,初始时刻有(2)式精确给出,第一个时间步的值由 (3)式给出,因此,须对 (3)式进行离散,下面推出与(9)和 (13)式相匹配的第一个时间步的离散格式.利用 Taylor展开式将在处展开可得利用(1)～(3)式,且略去高阶项,即可得与(9)式相匹配的第一个时间步的离散格式与(13)式相匹配的第一个时间步的离散格式2 稳定性分析下面,采用 Fourier分析方法对格式进行稳定性分析.引理[12] 实系数二次方程λ2-bλ-c=0的根按模不大于 1的充要条件为|b|≤1-c≤2.2.1 AD I(2,4)格式的稳定性定理 1 格式(9)是无条件稳定的.证明显然,格式(9)可写成用表示采用上述格式进行计算产生的误差,设源项 f无误差存在,则格式的误差项满足格式相应的齐次方程,即记其特征项表示为其中为虚数单位,ηn为第 n个时间层上的波幅,σ1、σ2为波数,令λ=τ/h,则可得误差的传播矩阵为其中矩阵的特征方程为所以有由于Px≥0、Py≥0、Qx≥0、Qy≥0,故|b|≤2,由引理可得格式(9)是无条件稳定的.2.2 AD I(4,4)格式的稳定性定理 2 格式(13)是条件稳定的,其稳定性条件为证明由相同的分析过程可得格式(13)误差的传播矩阵为矩阵的特征方程为由引理得即上式右端显然成立,考察左端可得解之可得|a|λ≤1/3,即格式(13)是条件稳定的,其稳定性条件为|a|λ≤1/3.3 数值验证对于问题(1)～(4),令问题的精确解为计算是用 Fortran77语言进行编程且在 Pentium IV/2.4G PC机上双精度制下进行的.由于2种格式所得线性方程组均为三对角线型,所以对每一步可以采用TDMA 算法.表 1给出了不同网格步长下,问题在 t=0.125时刻,当τ=h时二阶AD I格式[8]、FULL(4,4)格式[7]和τ=h2时AD I(2,4)格式的数值计算结果的最大误差 E和收敛阶 rate=ln(E1/E2)/ln2(E1和E2分别为粗网格及相邻的细网格上的最大误差).结果表明,3种格式均达到了各自的精度,并且AD I(2,4)格式的计算结果要比其它两种格式精确得多.表 2和表 3给出了不同网格步长下,问题在 t =0.25时刻,当τ =h/2时FULL(4,4)格式[7]、AD I(4,4)格式和τ=h2时 HWAL(2,4)[7]格式的数值计算结果的最大误差 E、收敛阶 rate和 CPU时间.结果表明,3种格式均达到了各自的精度,并且AD I(4,4)格式的计算结果要比其它两种格式精确得多,并且计算时间最少.表 4给出了 h=1/64时,问题在 t=3.2时刻,对不同网格比λ(此时τ=λh不同),3种AD I 格式的收敛性的比较.结果表明,当时,AD I(4,4)格式是发散的,即条件稳定的,其它格式仍然是收敛的,即无条件稳定的.这与本文的理论分析结果一致.本文的方法可推广到三维波动方程的数值求解中去,将另文报道.表 1 t=0.125时刻不同空间步长 h的最大误差和收敛阶Table 1 Max imum error and convergencerate att=0.125 for differenthh 二阶ADI格式[8]FULL(4,4)格式[7]ADI(2,4)格式E rate E rate E rate 1/8 2.81e-2 1.29e-34.81e-4 1/16 7.57e-3 1.89 8.20e-5 3.98 3.01e-5 3.97 1/32 1.93e-3 1.975.15e-6 3.99 1.88e-6 4.00 1/64 4.83e-4 1.99 3.22e-7 4.00 1.18e-7 4.001/128 1.22e-4 2.00 1.82e-8 4.14 7.35e-9 4.00表 2 t=0.25时刻不同空间步长 h的最大误差和收敛阶Table 2 Max imum error and convergencerate att=0.25 for differenthh HWAL(2,4)格式[7]FULL(4,4)格式[7]ADI(4,4)格式E rate E rate E rate 1/8 1.54e-2 2.47e-4 2.87e-5 1/162.46e-5 9.29 1.56e-53.98 1.73e-64.05 1/32 1.53e-6 4.00 9.75e-7 4.001.07e-7 4.01 1/64 8.96e-8 4.09 5.99e-8 4.02 6.68e-9 4.00 1/128 5.60e-9 4.002.73e-9 4.45 4.17e-10 4.00表 3 t=0.25时刻不同空间步长 h的 CPU时间Table 3 CPU t ime att=0.25 for differenthλ HWAL(2,4)格式[7]FULL(4,4)格式[7] ADI(4,4)格式1/8 ＜10e-8 ＜10e-8 ＜10e-8 1/16 0.125 0.015 0.015 1/32 1.985 0.125 0.031 1/64 57.50 1.391 0.219 1/128 2830.51 28.39 1.687表 4 当 h=1/64,t=3.2时刻,对不同λ的最大误差Table 4 Max imum erroratt=3.2 andh=1/64 for differentλλ 二阶格式[8] ADI(2,4)格式 ADI(4,4)格式0.4 6.78e-4 4.68e-4 1.93e-8 0.8 2.04e-3 1.84e-3 1.49e+51 1.6 7.03e-3 6.87e-3 3.09e+65 3.2 1.96e-2 1.95e-2 2.80e+39参考文献[1]李胜坤,冯民富,李珊.Benjamin-Bona-Mahony方程的有限差分近似解[J].四川师范大学学报：自然科学版,2003, 26(4)：363-365.[2]吕胜关.一类高阶方程组的差分方法[J].郑州大学学报：理学版,1999,32(1)：12-19.[3]罗明英,舒国皓,王殿志.RLW方程的有限差分逼近[J].四川师范大学学报：自然科学版,2001,24(2)：138-143.[4]Visher J,Wandzura S,White A.Stable,high-order discretization forevolution of the wave equation in 1+1 dimensions[J].J ComputPhys,2004,194：395-408.[5]Wandzura S.Stable,high-order discretization for evolution of the wave equation in 2+1 dimensions[J].J Comput Phys,2004, 199：763-775.[6]葛永斌,吴文权,田振夫.二维波动方程加权平均隐格式及其多重网格方法[J].上海理工大学学报,2002,24(3)：205-208.[7]葛永斌,吴文权,田振夫.二维波动方程高精度隐式格式及其多重网格方法 [J].厦门大学学报：自然科学版,2003, 42(6)：691-696.[8]孙志忠.偏微分方程数值解法[M].北京：科学出版社,2005.[9]Karaa S,Zhang J.High orderADImethod for solving unsteady convection-diffusion problem[J].J Comput Phys,2004,108：1-9.[10]葛永斌,田振夫,吴文权.高维热传导方程的高精度交替方向隐式方法[J].上海理工大学学报,2007,29(1)：55-58.[11]Hirsh R S.Higher order accurate difference solutions of fluid mechanics problem by a compact differencing technique[J].J Comput Phys,1975,19：90-109.[12]陆金甫,关冶.偏微分方程数值解法[M].北京：清华大学出版社,1987.。

二维有限差分矩阵-概述说明以及解释1.引言1.1 概述引言部分会总体介绍本文将要讨论的主题——二维有限差分矩阵。

本文将首先简要概述二维有限差分方法的基本原理和应用领域，然后详细介绍二维有限差分矩阵的构建方法。

通过本文的阐述，读者将了解到二维有限差分矩阵在数值计算、物理仿真、图像处理等领域的广泛应用，并获得一定的实践指导和理论支持。

二维有限差分方法是一种常用的数值计算技术，广泛应用于解决二维偏微分方程及相关问题。

通过将连续问题离散化为离散点之间的差分，可以利用计算机进行高效且准确的计算。

而二维有限差分矩阵则是二维有限差分方法中的关键组成部分，用于描述问题的离散化形式。

本文着重介绍二维有限差分矩阵的构建方法。

首先，将介绍二维有限差分方法的基本原理，包括空间离散化和时间离散化。

然后，将详细介绍如何根据实际问题的边界条件和离散化精度构建二维有限差分矩阵。

通过合理选择差分格式和边界条件，可以得到满足精度要求的二维有限差分矩阵。

需要注意的是，二维有限差分方法和二维有限差分矩阵的适用范围广泛，不仅仅局限于数值计算领域。

它还可以应用于物理仿真领域，如电磁场模拟和流体动力学分析；以及图像处理领域，如边缘检测和图像恢复等。

通过本文的学习，读者将能够掌握二维有限差分方法的基本原理，了解二维有限差分矩阵的构建方法，并在实际应用中灵活运用。

1.2 文章结构本文共分为引言、正文和结论三个部分。

在引言部分，首先对二维有限差分方法做了一个概述，介绍了其在科学计算和工程领域中的重要性和广泛应用。

接着对文章的结构进行了说明，明确了各个部分的内容和安排。

最后，明确了本文的目的，即探讨二维有限差分矩阵的构建方法。

正文部分主要包括两个部分：二维有限差分方法简介和二维有限差分矩阵的构建。

在第2.1节中，我们将对二维有限差分方法进行简要介绍，包括其基本原理和步骤。

我们将详细解释如何将连续的偏微分方程转化为离散的代数方程，并介绍如何选择合适的差分格式和网格划分方法。

定常对流扩散反应方程非均匀网格上高精度紧致差分格式田芳;田振夫
【期刊名称】《工程数学学报》
【年(卷),期】2009(026)002
【摘要】本文构造了非均匀网格上求解定常对流扩散反应方程的高精度紧致差分格式.我们首先基于非均匀网格上函数的泰勒级数展开,给出了一阶导数和二阶导数的高阶近似表达式;然后将模型方程变形,借助于对流扩散方程高精度紧致格式构造的方法,结合原模型方程,得到定常对流扩散反应方程的高精度紧致差分格式;最后给出的数值算例验证了本文格式高精度和高分辨率的优点.
【总页数】7页(P219-225)
【作者】田芳;田振夫
【作者单位】宁夏大学数学计算机学院,银川,750021;复旦大学数学科学学院,上海,200433
【正文语种】中文
【中图分类】O241.82
【相关文献】
1.一维非定常对流扩散方程非均匀网格上的高精度紧致差分格式 [J], 黄雪芳;郭锐;葛永斌
2.一维定常对流扩散方程非均匀网格上的高阶紧致差分格式 [J], 薛文强;兰斌;葛永斌
3.一维非定常对流扩散反应方程的高精度紧致差分格式 [J], 杨晓佳;田芳
4.一维定常对流扩散反应方程的高精度紧致差分格式 [J], 祁应楠;武莉莉
5.一维非定常对流扩散方程非均匀网格上的高阶紧致差分格式 [J], 赵飞;陈建华;葛永斌
因版权原因，仅展示原文概要，查看原文内容请购买。

poisson方程三维有限差分格式三维Poisson方程有限差分格式主要应用于求解三维空间中的Poisson方程。

与二维情况类似，我们需要将三维空间划分为网格，然后对网格节点上的函数值进行差分。

以下是一个基本的三维有限差分格式求解过程：1. 网格划分：首先对三维求解区域进行网格划分。

网格划分的方向可以采用均匀网格或非均匀网格，取决于问题的特性。

通常，在边界附近的网格节点密度会较大，以更好地捕捉边界附近的梯度变化。

2. 建立差分方程：根据五点差分格式，我们可以得到三维Poisson方程的差分形式。

在x、y、z方向上，分别对函数u(x, y, z)进行差分，得到如下形式的差分方程：u(x+h, y, z) - u(x-h, y, z) / (2h) = λ* (u(x, y+h, z) - u(x, y-h, z)) / (2h) u(x, y+h, z) - u(x, y-h, z) / (2h) = λ* (u(x, y, z+h) - u(x, y, z-h)) / (2h) u(x, y, z+h) - u(x, y, z-h) / (2h) = λ* (u(x+h, y, z) - u(x-h, y, z)) / (2h)其中，h为网格步长，λ为比例系数，可根据边界条件和初始条件进行调整。

3. 迭代求解：将差分方程组转化为矩阵形式，然后采用迭代方法（如Gauss-Seidel迭代法）求解。

对于每个网格节点，迭代更新u(x, y, z)的值，直到达到预设的迭代次数或满足收敛条件。

4. 后处理：在求解过程中，可以采用一些后处理方法来提高解的质量，如欠松弛技术、人工粘性层等。

5. 验证与分析：将求解得到的结果与理论解析解或实验数据进行比较，分析数值解的准确性和稳定性。

需要注意的是，在实际应用中，根据问题的具体情况，可能需要对上述求解过程进行相应的调整，如采用非均匀网格、多重网格技术、自适应步长等方法。

二维有限差分法
二维有限差分法是一种用于求解二维偏微分方程的数值解法。

它将待求解区域分割成有限个网格点，并利用差分近似方法将偏微分方程转化为代数方程组，然后通过迭代求解这个方程组来获得数值解。

具体来说，二维有限差分法将二维区域 $\Omega$ 划分成
$M$ 个横向离散点和 $N$ 个纵向离散点，得到一个 $M \times N$ 的网格。

偏微分方程在网格上被离散化为一组代数方程，其中每个网格点的解被近似表示为该点以及周围点的函数值。

在二维有限差分法中，常用的差分格式包括中心差分、向前差分和向后差分等。

通过差分近似，偏微分方程中的导数被转化为差分系数的线性组合。

然后，可以得到一个线性方程组，其中每个网格点的系数由该点周围网格点的差分系数决定。

解这个线性方程组可以使用迭代方法，如Jacobi迭代、Gauss-Seidel迭代或SOR（逐次超松弛法）迭代等。

迭代过程一般需要设定迭代停止条件，比如迭代次数的上限、残差的收敛精度等。

通过二维有限差分法，可以求解各种边界条件下的二维偏微分方程，比如泊松方程、热传导方程、扩散方程等。

它是一种经典且简单实用的数值方法，广泛应用于科学计算和工程领域。

五点差分格式求解poisson方程在数学和科学领域中，Poisson方程是一种常见的偏微分方程，通常用于描述在某个区域内的物理现象。

求解Poisson方程在科学计算和数值模拟中具有重要的应用。

本文将介绍一种常用的数值求解方法——五点差分格式，用于求解Poisson方程。

1. 引言Poisson方程是二阶偏微分方程，形式如下：∇²u = f(x, y)其中，u是未知函数，f(x, y)是已知函数，∇²是拉普拉斯算子。

求解Poisson方程的目标是找到满足方程的函数u。

2. 五点差分格式五点差分格式是一种常用的离散化方法，用于将连续的Poisson方程转化为离散的数值方程。

在二维情况下，我们可以将区域划分为若干个格点，然后利用差分近似来求解。

假设我们在二维区域(0, 1) x (0, 1)内选择了一组均匀的格点(xi, yj)，其中i和j分别表示x和y方向上的索引。

则方程中的拉普拉斯算子可以被近似为：(∂²u/∂x²) ≈ (u(i+1,j) - 2u(i,j) + u(i-1,j))/Δx²(∂²u/∂y²) ≈ (u(i,j+1) - 2u(i,j) + u(i,j-1))/Δy²其中Δx和Δy是格点的空间间隔。

将上述近似代入Poisson方程，我们可以得到离散的数值方程：(u(i+1,j) - 2u(i,j) + u(i-1,j))/Δx² + (u(i,j+1) - 2u(i,j) + u(i,j-1))/Δy² = f(i,j)3. 离散方程求解根据五点差分格式，我们可以得到离散的数值方程，进而求解Poisson方程。

为了方便计算，我们可以定义一个N x N的矩阵A，其中N表示在x和y方向上的格点数目。

矩阵A中的每个元素对应于方程中的一个未知数。

我们可以将方程表示为矩阵形式：Au = b其中，u是未知函数的向量，b是已知函数f的向量。

一类二维稳态对流——扩散方程的有限差分法一维稳态扩散方程描述了物质在一维空间中的扩散行为。

然而，在某些情况下，我们需要研究物质在二维平面中的扩散行为，例如热传导、流体传输等。

本文将介绍一类二维稳态对流-扩散方程的有限差分法。

二维稳态对流-扩散方程可以写作：∇·(D∇u) + ∇·(cu) + fu = 0 —— (1)其中，D是扩散系数，c是速度场，u是待求解的物理量，f是源项。

在这个方程中，第一项表示物质的扩散项，第二项表示对流项，第三项表示源项。

我们需要求解方程(1)，找到u的分布。

为了应用有限差分法来求解二维稳态对流-扩散方程，需要将二维空间离散化为一个网格。

假设我们将x方向离散为Nx个等距的节点，y方向离散为Ny个等距的节点，那么我们可以得到一个(Nx+1)×(Ny+1)的网格。

我们在网格节点上定义未知量u，然后将方程(1)对节点处的u进行离散化。

首先，我们对方程(1)的扩散项进行离散化。

我们使用五点差分格式来近似二维Laplace算符∇·(D∇u)。

对于网格节点(x,y)，我们可以得到以下差分格式：(Dij(xi+1,yj)ui+1,j + Dij(xi-1,yj)ui-1,j +Dij(xi,yj+1)ui,j+1 + Dij(xi,yj-1)ui,j-1 -4Dij(xi,yj)ui,j) / ∆x^2 + (Dij(xi,yj)ui,j) / ∆y^2其中，∆x和∆y是网格步长，Dij是扩散系数。

接下来，我们对方程(1)的对流项进行离散化。

我们使用中心差分格式来近似二维梯度算符∇·(cu)。

对于网格节点(x,y)，我们可以得到以下差分格式：(cxi+1/2,yj(ui+1,j - ui,j)) / ∆x + (cxi-1/2,yj(ui,j - ui-1,j)) / ∆x + (cyi,j+1/2(ui,j+1 - ui,j)) / ∆y + (cyi,j-1/2(ui,j - ui,j-1)) / ∆y其中，cxi+1/2,yj、cxi-1/2,yj、cyi,j+1/2和cyi,j-1/2是速度场在节点(x,y)处的中心点处的x和y分量。

方程utt-uxx-uxxtt=f(u)的紧致差分格式紧致差分格式是一种在数值计算中比较常用的方法，用于解决求解常微分方程的问题。

本文将讨论如何使用紧致差分格式来求解ut-uxx-uxxtt=f(u)的常微分方程。

首先，我们来看看ut-uxx-uxxtt=f(u)的常微分方程的几何意义。

这个方程的左边表示的是一个变量u的二阶时间导数，其中u的二阶空间导数也参与其中。

右边的f(u)表示的是一个函数，我们可以认为它是外部的一个影响因素。

接下来，我们要使用紧致差分格式来求解ut-uxx-uxxtt=f(u)的常微分方程。

首先，我们将方程进行分析，可以得出：$$u_{tt}-u_{xx}-u_{xxtt}=f(u)$$可以将上述方程分解为：$$u_{tt}-u_{xx}=g(u)$$$$g(u)=-u_{xxtt}+f(u)$$此时，我们可以使用紧致差分格式来求解上述方程，即：$$\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{\Delta x^2}-\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{\Deltax^2}=g(u_{i,j})$$其中，$u_{i,j}$表示的是时间和空间上的网格点的u值，$\Delta x$表示的是网格的间距，$g(u_{i,j})$表示的是外部影响因素f(u)在网格点$u_{i,j}$上的值。

最后，我们可以使用上述紧致差分格式来求解ut-uxx-uxxtt=f(u)的常微分方程，其中$u_{i,j}$表示时间和空间上的网格点的u值，$\Delta x$表示的是网格的间距，$g(u_{i,j})$表示的是外部影响因素f(u)在网格点$u_{i,j}$上的值。

使用紧致差分格式，可以很容易地求解出ut-uxx-uxxtt=f(u)的常微分方程的解。

总的来说，紧致差分格式是一种比较常用的数值计算方法，可以用来求解ut-uxx-uxxtt=f(u)的常微分方程。

大气数值模式及模拟（数值天气预报）习题第一章大气数值模式概论1．试述原始方程组、全球模式、区域模式和非静力模式之间的区别。

2．试述天气模式、气候模式的主要区别？3．区域气候模式、大气环流模式、中尺度模式、陆面模式、边界层模式各有什么特点？第二章 大气运动方程组1. 试证明球坐标系中单位矢量i 的个别变化率为(sin cos )cos di u j k dt r ϕϕϕ=- 2.试说明局地直角坐标系（即z 坐标系）中的运动方程与球坐标系中的运动方程有何异同？3．用球坐标导出下面两个方程：(sin cos )cos d i u j k dt r ϕϕϕ=- tan d j u v i k dt r rϕ=-- 4.由热力学方程v dT d C p Q dt dtα+=推导出如下方程： p dT C Q dt αω-= ()dp dtω= 式中v dT C dt为单位质量理想空气内能的变化率，v C 为空气的定容比热，d p dtα为可逆过程中单位质量非粘性气体在单位时间里膨胀所作的功。

Q 为外界对单位质量空气的加热率。

第三章 数值计算方案1. 什么是差分格式的收敛性和稳定性？二者之间有何关系？2. 试证明一阶偏微商u x ∂∂的三点差商近似式：3(,)(,)213(,)4(,)(2,)22u u x x t u x t x x u x t u x x t u x x t x ∂+∆-⎡⎤=⎢⎥∂∆⎣⎦-++∆-+∆⎡⎤-⎢⎥∆⎣⎦的截断误差为2()O x ∆。

3. 用中央差分将涡度方程()()()l l u u u v l t x y x y∂Ω∂Ω+∂Ω+∂∂++=-+∂∂∂∂∂ 写成有限差形式。

设(,)l l x y =，并取水平坐标步长为s δ，时间步长为t δ。

4. 分别对x 轴上的i+1和i+3格点，以d 和2d 为步长，写出一阶微商dF dx的前差、后差和中央差的差分近似式，以及二阶微商22d F dx 的二阶中央差分近似式。

二维泊松方程边值问题有限差分方程的病态结构和最优预条件子张衡;郑汉垣【摘要】基于结构分析的思想,讨论大规模病态稀疏线性方程的病态机理和预处理原理,定义该方程组的病态结构、病态因子、去病因子.针对病态结构,设计去病因子,以去病因子为预条件子,并对预条件子的性能进行定量分析,结果表明去病因子是最优预条件子,该预条件子的使用,几乎不增加迭代的计算量,预处理后方程组的主体保持正定对称,条件数接近常数.【期刊名称】《福建师大福清分校学报》【年(卷),期】2018(000)005【总页数】6页(P1-6)【关键词】病态机理;病态结构;病态因子;去病因子;预处理【作者】张衡;郑汉垣【作者单位】福建师范大学福清分校无损检测技术福建省高校重点实验室,福建福清350300;福建师范大学福清分校电子与信息工程学院,福建福清350300;龙岩学院传播与设计学院,福建龙岩364012【正文语种】中文偏微分方程大规模数值求解问题，通常转化为大规模病态（高条件数）稀疏线性方程组的求解 [1].该方程组的条件数经常随着问题规模的增加而增加[2]，成为影响求解效率和精度的瓶颈因素.因此，在迭代求解之前，使用预处理方法来减少方程组的病态，成为提高求解效率和精度的必要措施.所谓“预处理技术”是指在求解方程组时，构造简单的可逆矩阵 M 1 , M 2，使得Cond( ) < <Cond(A )（Cond(A)是指矩阵A的条件数），, 容易计算，从而方程组（1）化成等价易解的方程组其中x ′ =M2 x，矩阵 M 1 M 2称为预条件子 [3-4]．如果 C ond( 1)与A的阶数无关，则称M1 M 2为最优预条件子[5]205-209.病态方程组的成功求解，常常以适当的预条件子作为前提.由于缺乏对病态机理和预处理原理的研究，目前关于预处理问题的理论仍然不完善.一是缺乏一般的预处理方法，没有通用的预条件子，只能针对具体问题，根据预条件子的基本要求，设计具体的预条件子[6-8]；二是对预处理的效果（条件数下降的程度，预处理后的条件数，对计算量的影响等）缺少科学的定量分析，多用实验结果说明[9-15].目前关于病态机理和预处理原理的研究鲜见有成果发表.本文讨论病态机理和预处理原理.定义方程组的病态结构、病态因子、去病因子，讨论他们的性质和作用. 使用有限差分方法，大规模求解二维泊松方程边值问题时，基于非均匀网格形成的有限差分方程，是稀疏病态方程组.基于结构分析的思想[16-19]，针对该方程，研究病态结构、病态因子、去病因子.将病态分为由病态因子表达的原发性病态和由具体问题和方法引起的继发性病态，更准确地说明了不同的因素对病态的影响；以去病因子为预条件子，对预处理的效果（条件数下降的程度，对计算量的影响）定量分析，结果说明去病因子是最优预条件子[5]205-209；该预条件子的使用，几乎不增加求解的计算量，预处理后方程组的主体保持正定对称，条件数接近常数（与问题规模无关）.1 矩阵的病态结构与去病因子1.1 矩阵的病态结构定义1 可逆矩阵的如下结构称为病态结构：其中是整数，，Y1C Y2 ,Y1Y2 均可逆，Cond(C )很小，C o nd(Y 1C Y2 )和 C ond(Y 1Y2)很大为A的病态主体，分别称 Y1 , Y2 为左、右病态因子.显然，在定义1中，C o nd(A) ≈Cond(Y 1C Y2)1.2 去病因子如果A有(3)式的病态结构，为减少或者消除病态，可针对病态因子，设计预条件子M1 , M 2，使得) < <Cond(Y 1C Y2 )，从而减少或者消除病态因子 1 2,Y Y 的致病作用.因为 1 2,Y Y 不是方阵，所以即使找到病态因子，也难以用逆矩阵的方法设计预条件子.特别地，有定义2 如果可逆矩阵 1 2,M M 满足：，或者，则称这样的 1 2,M M 为 1 2,Y Y 的去病因子.定义2中可认为去病因子 1 2,M M 消除了病态因子 1 2,Y Y 的致病作用.通常考虑寻找满足=I的可逆矩阵 1 2,M M.特别地，对于正定对称矩阵的情形，有命题 1：设C分别是阶正定对称矩阵，则有2)如果有α阶矩阵M满足则有证明：1)利用Rayleigh-Ritz定理[20]所以，即(4)成立.2)根据条件有 M -1 YYTT M-TT =，I是α α阶单位矩阵，利用结论1)有Cond(M -1 YCYTT M-TT )≤≤Cond(C ) Cond(M -1 YYTT M -TT )=Cond(C)，即(6)成立.根据命题1，如果 Y , 是病态因子，则满足(5)式的矩阵 M , 可以消除 Y , 的作用，因此有定义3 称满足(5)式的矩阵M为属于Y的去病因子.1.3 一个特别的病态因子及其去病因子阶正弦变换矩阵，有命题2：Z Z TT =HHTT命题3：3)对于任意 m ( n + 1 )+ n ( m + 1 )阶正定对称矩阵P，有证明：1)根据命题2 ，Z Z TT 的特征值为，1 ≤ i≤ n，1 ≤ j≤ m.所以即(7)成立.2)记，则容易验证：1 + 2min(x, y) ≥ g ( x, y) ≥ m in(x, y)，所以3)根据(4)可得.证毕.根据命题3和定义1，Z , Z T 是 Z PZT或者 Z Z T 的病态因子；根据命题2和定义3，矩阵H是属于Z的去病因子.2 二维泊松方程边值问题非均匀网格的有限差分方程使用有限差分方法，求解二维泊松方程边值问题[4]100-101其中f =f( x, y) ,(x, y) ∈ D = [ a, b] × [ c, d ] ⊂R2对D做非均匀网格剖分[5]59：在(xi, y j )，使用δ x x u i, j , δ y y u i, j , fi, j 分别代替u x x , u y y,f，得到非均匀差分格式其中0≤≤i≤n+1，0≤≤j≤≤ m +1.记则有命题4：问题(9)有如下的有限差分网格方程证明：容易验证根据上述记号的定义，有根据差分格式(10)，有 - δxx u -δyyu=f，所以(11)式成立.证毕.3 有限差分方程的病态结构与最优预条件子及其预处理效果显然，当 h i , x , h j , y都很小时，根据 P , Q, Z的定义，命题4的(11)式中，ZPZ T >>Q，在系数矩阵A中，Z P Z T是A的大范数部分，即主要部分，Co nd(A) ≈ Cond(Z PZ T )；根据定义1和命题3，有限差分方程(11)的系数A矩阵具有病态结构，Z P Z T是A的病态主体，Z, Z T是病态因子；根据命题2，矩阵H是属于Z的去病因子.根据命题3，A的病态大部分是由病态因子Z表达的，少部分是P表达的，Q对A的病态影响很小.病态因子Z表达的病态来自微分算子，是本质的，离散精度越高，A的阶数越大，Z表达的病态越严重；P表达的病态来自网格大小，几乎不受A的阶数影响，是非本质的，可以随着应用问题的不同而不同，在应用中可以调整. 定义4 称病态因子Z表达的病态为原生病态，P表达的病态为继发病态使用H作为预条件子，则方程(11)化成容易验证：Cond()P 几乎不受 ,m n的影响.根据命题1的结论2),有因此，H H T 是最优预条件子[5]205-209.H, H - 1都是离散傅里叶变换矩阵的直积与对角矩阵的积，有简单、确定的结构，因此预条件子的构造和预处理计算不需要增加大量成本，其中预处理需要的计算是H , H -1与向量的乘积，每次需要的计算操作数是O( m n l og2(n m )) ，即有限次离散快速傅里叶变换的计算量，所以每个迭代步的计算量没有显著增加.预处理后，方程(12)的系数矩阵的主体仍然是正定对称矩阵，所以仍然可以使用经典Krylov子空间方法—共轭梯度法求解，即为预处理共轭梯度法[4]139-151.4 结论1）二维泊松方程边值问题的有限差分网格方程，是稀疏病态方程组，它的结构中，隐藏着包含病态因子的病态结构，这是病态产生的根本原因.2）有限差分方程组的病态分为由病态因子表达的原发性病态和由具体方法引起的继发性病态，原发性病态是主要的、本质的.3）针对病态因子，可找到去病因子，将去病因子作为预条件子，可消除病态因子的致病作用，即消除原发性病态.预处理后，系数矩阵主体（大范数部分）的条件数降为接近常数，几乎不受方程阶数的影响，因此，去病因子是最优预条件子[5]205-209.4）将去病因子作为预条件子，预处理过程基本不增加计算操作数，并且预处理后矩阵主体保持正定对称.5）与经典的不完全LU分解法比较，本文的预条件子针对病态因子设计，有明显的针对性和标准，得到的预条件子与具体的问题的继发性病态没有关系，有一定的通用性.【相关文献】[1] JIA Z X. The convergence of krylov subspace methods for large nonsymmertric linear systems [J]. Acta Mathematical Sinica, 1988, 14(4): 507-518.[2] 吴勃英,王德明,丁效华,等.数值分析原理[M].北京：科学出版社，2003:71-101.[3] 张科.两类线性方程组的预处理技术及数值求解方法[D].上海：上海大学，2014:7-9.[4] 李荣华 . 偏微分方程数值解法 [M]. 2 版 .北京：高等教育出版社 , 2010.[5] 张文生. 科学计算中的偏微分方程有限差分法[M].北京：高等教育出版社，2006.[6] BAI Z Z, LI G Q. Restrictively preconditioned conjugate gradient methods for systemsof linear equations [J]. IMA Journal of Numerical Analysis, 2003 (23): 561–580.[7] WANG Z Q. Restrictively preconditioned chebyshev method for solving systems of linear equations [J]. Eng Math, 2015 (93):61-76.[8] BRU R, MARW J, MAS J, et al. Preconditioned iterative methods for solving linear least squares problems [J].SIAM J Sci Comput, 2014, 36(4):2002-2022.[9] 于春肖 , 苑润浩 , 穆运峰 . 新预处理 ILUCG 法求解稀疏病态线性方程组[J]. 数值计算与计算机应用,2014, 35(1): 21-27.[10] XUE Q F, GAO X B, LIU X G，Comparison theorems for a class of preconditioned AOR iterative methods[J].Journal of Mathematics, 2014, 34(3): 448-460[11] 潘春平,马成荣 ,曹文方,等 . 一类预条件AOR 迭代法的收敛性分析 [J]. 数学杂志 , 2013, 33(3): 479-484.[12] 罗芳.L-矩阵的多参数预条件AOR迭代法[J]. 数学的实践与认识 , 2013, 43(15): 277-282.[13] 吴建平,赵军,马怀发,等.一般稀疏线性方程组的因子组合型并行预条件研究[J]. 计算机应用与软件 , 2012, 29(5): 6-9+108.[14] 李继成,蒋耀林. 预条件IMGS迭代方法的比较定理 [J].数学物理学报 , 2011, 31A(4): 880-886.[15] PAN C P. A new effective preconditioned Gauss-Seidel iteration method[J]. Journal on Numerical Methods and Computer Applications, 2011, 32(4): 267-273[16] 张衡. 一维问题一次形函数有限元方程的条件数与预处理[J]. 福建师大福清分校学报,2016(2):1-3.[17] 张衡. 一维问题m次lagrange形函数有限元方程的条件数与预处理[J]. 石河子大学学报（自然科学版）, 2017, 35(4):518-521.[18] 张衡. 两点边值问题2次lagrange形函数有限元方程的条件数和预处理[J]. 福州大学学报（自然科学版）, 2017,45(5): 617-622[19] 张衡. 两点边值问题3次lagrange形函数有限元方程的条件数和预处理[J]. 计算力学学报, 2017,34(5): 672-676。

§10. 高阶紧致差分格式先考虑导数的差分近似。

若某一差分近似的精度是 p 阶的，则近似的误差就是 ()p h O 。

要想进一步提高精度，通常有两种途径：减小 h （h -version ）或是提高 p （p -version ）。

但由于计算机资源的限制，h 不可能无限地减小，因此在需要高精度流场计算的情形（如，粘性边界层、湍流等），就要考虑采用高阶格式。

通常情形，构造高阶格式需要更多的点。

例如：两点差分近似()()()f x h f x f x h+-¢»只有一阶精度。

而使用三个点，就可以构造出二阶近似()()()()2432f x h f x h f x f x h-+++-¢»精度越高，需要的点就更多。

对于中心差分近似也有类似的结果。

但是这种高阶近似用在差分格式中，除了计算公式更加复杂，计算量增加之外，还会造成其他困难。

例1：以一个简单的常微分方程初值问题为例。

设 0a > 。

0duau dx+= （01x < ） ， ()0u =α 取 M 个网格，空间步长 1h M= ，网格点记作 j x jh=（0,1,2,,j M =L ），网格点上的近似解记作 ()j j u u x » 。

因 0a > ，导数采用向后差分近似，就有10j j j u u au h--+= （1,2,3,,j M =L ）实际的计算方案为0u =α ， 111j j u u ha-=+ （1,2,3,,j M =L ）上述格式用到两个点，但只有一阶精度。

如果采用二阶差分近似，则成为12340j j j j u u u au h---++= （2,3,,j M =L ）这个格式具有二阶精度。

可是由于涉及三个点，所以只能从 2j = 开始计算。

而初始条件只提供了 0u =α 。

因此 1u 的计算就需要补充另外的等式。

对于更为复杂的流动控制方程以及更复杂、精度更高的数值格式，这种问题就更加严重。