双调和方程的有限积分方法
李书伟;徐定华;余跃
【摘 要】利用有限积分法求解平面矩形区域双调和方程边值问题.首先,对双调和方程以及边界条件分别进行积分,得到一带有任意函数的线性常微分方程组;其次,将积分产生的任意函数分别进行插值估计,进而转化成为一可求解的线性代数方程组;最后,利用正则化方法求解奇异线性方程组,获得近似解误差估计.通过Matlab进行数值模拟实验获得数值结果,并进行误差分析.数值结果表明,与有限差分法、有限元法以及广义有限差分法相比较,有限积分法具有更高精度.
【期刊名称】《浙江理工大学学报》
【年(卷),期】2016(035)001
【总页数】7页(P133-139)
【关键词】双调和方程;有限积分法;正则化方法;误差估计;数值解
【作 者】李书伟;徐定华;余跃
【作者单位】浙江理工大学理学院,杭州310018;浙江理工大学理学院,杭州310018;浙江理工大学理学院,杭州310018
【正文语种】中 文
【中图分类】O242.2;O241.8;O343.1
双调和方程主要来源于流体力学和弹性力学，在工程技术方面应用极为广泛，相关理论和方法的研究一直是人们关注的热点，有限差分法[1-2]和有限元法[3]是常用的求解双调和方程的两种数值解法。有限差分法的优点是方法简单、计算量小，但是对边界区域规则性的要求非常高，难以处理复杂网格。有限元法虽然可以针对任何边界区域，但对于区域网格的划分有严格的要求。有限体积法[4]和边界元法[5]也是求解双调和方程非常重要的数值解法。有限体积法虽然自身包含几何信息、易处理复杂网格，但是计算比较复杂、不易提高精度。边界元法最大的优点是，可以使问题的空间维数降低，从而使计算工作量及其所需计算机容量大大减小，但是该方法需要求出问题的基本解，推导过程比较复杂。无网格法[6-8]是近年来新兴的一种数值方法，它是针对任何边界区域，并且只需在计算区域上选取若干节点及其径向基函数，不需要划分网格。
目前，偏微分方程边值问题的求解方法大多为有限差分法、有限元法、有限体积法和边界元法。迄今为止，对应用有限积分法求解偏微分方程的研究工作还不多，在该方法理论上的研究很少。本文研究了双调和方程边值问题的有限积分求解方法、数值算法、数值模拟以及该方法的收敛性分析，并通过模拟平面薄板受一均匀载荷时，研究静态平衡下薄板与载荷接触面上的位移分布与误差分析。本文针对研究双调和方程边值问题的求解方法，拟采用Wen等 [9]提出的有限积分法对该问题进行数值求解，考虑了该方法的理论基础——积分的离散与函数的插值表示，结合文献[9]给出的利用有限积分法求解一维偏微分方程边值问题的解题思路，以求出双调和方程的数值解，并从理论上对该方法的收敛性做了分析。
1.1　积分离散形式
下面考虑一维情形[9]。考虑定义在[a,b]上的函数u(x)，定义域离散为a=x0<x1<…<xN=b，其中：xk=a+kh，h=(b-a)/N,k=0,1,…,N。设u(x)关于基函数{φi(x),i=0,1,…,N}的插值格式为：
则一次积分可离散为：
U(1)(xk),
k=0,1,…,N
式(2)用矩阵表示为：
其中：
A(1),
u=[u0,u1,…，uN]T.
同理，二次积分的离散形式为：
更一般地，m次积分的离散结果为:
下面考虑二维情形[10]。考虑函数u(x,y)，其定义域为，用M=(N1+1)×(N2+1)个均匀网格点离散，即：
为方便表述，把网格节点标记为：
称这样的顺序为标准顺序。
记Ux(x,y)和Uy(x,y)分别为对x和y的积分，则：
,
s=0,1,2,…,N2，
用矩阵表示为：
具体可写为：
，
或简记为：
其中A(1)是由方程(3)求得的(N1+1)×(N1+1)维积分矩阵。
关于x的m次积分的离散形式为：
同理，关于y的积分的离散形式为：
为了将其转化为标准顺序，可进行一个初等行交换，其对应初等矩阵T=(Tmn)，满足：
.
其中：i=0,1,…,N1;j=0,1,…,N2
1.2　函数的插值表示
下面考虑线性插值[9]。对于1.1节给出的函数u(x)进行线性插值，即使用梯形规则，得：
则：
U(1)(xk)=
,
用矩阵表示为：
对于二阶线性插值有：
对于m阶线性插值有：
正则化定义
下面考虑径向基插值[9-10]。考虑函数u(x)，其定义域为[a,b]，其上任一点的函数值都可以由一些随机分布的节点T，i=0,1,2,…,N,的节点值T插值得到。径向基函数的插值形式为：
u
=R(x)α+P(x)β
其中：Ri(x)是径向基函数，αi是Ri(x)的系数，Pj(x)是多项式基函数，βj是Pj(x)的系数，径向插值的一个附加条件为：
联立方程(5)和(6)得：
则方程(7)的一阶积分插值为：
，
因此，其一阶的有限积分矩阵为：
m阶积分的有限积分矩阵为：
考虑如下双调和方程边值问题：
其中，Ω=(a,b)×(c,d)为平面矩形区域，Γ为其边界，n为Γ的单位外法向向量，由文献[11]可知，在适当条件上述定解问题至多只有一个解。
对式(8)中方程的x和y分别进行四次积分有：
x2f1(y)+xf2(y)+f3(y)+y3g0(x)+
y2g1(x)+yg2(x)+g3(x)=
dξ4dξ3dξ2dξ1dη4dη3dη2dη1
其中fi(y),gi(x),i=0,1,2,3是由积分产生的任意函数。为了处理这些任意函数，我们对其进行插值估计。设为它们的插值节点，{φj(x)}为插值函数，有：
,
,
,
).
根据有限积分，方程(9)可离散为：
+X3Φyρ(0)+X2Φyρ(1)+XΦyρ(2)+Φyρ(3)
+Y3ΦxΛ(0)+Y2ΦxΛ(1)+YΦxΛ(2)+ΦxΛ(3)
其中：,,
,
,
,
,
,
,
,
.
则方程(10)可以根据边界条件确定唯一解，并且其导数边界条件可以表示为：
考虑方程(10)，令：
,
,
,
,
则得到如下线性方程组：
简记为
其中。
假设上式右端有一输入误差δ，使得≤δ,则有：
为避免因系数矩阵K的奇异性而导致错误的数值解，采用正则化方法解方程组(13)，并对有限积分法计算结果做出误差估计。
引理1　设，假设k>1且,则存在常数C=C(k)，使得：
证明可参见文献[12]中定理3.1。
定理1　考虑函数u(x,y)，其定义域为[0,1]×[0,1]，用M=(N1+1)×(N2+1)个均匀网格点离散，假设且误差δ>0，取正则化参数，则误差估计为：
.
特别的，当δ→0时，则有如下误差估计：
这里‖·‖表示L2-范数。
证明：用Tikhonov正则化方法求解线性方程组(13)，出其极小元，并将其代入方程2中，则有：
令，则有：).
对K进行奇异值分解得：
则，有正则解：
令Rα=(αI+KTK)-1KT，由文献[13]可知：
由此可得：
.
根据引理1得：
.
又因为：
,
,
则：
.
令，则误差估计为：
.
特别地，令δ→0且α=hμ-3/E(μ>3)，则误差估计为：
由定理1可知对于固定的σn，当δ→0时，→0；特别的，当δ→0时，对于足够光滑的函数，其光滑性μ也随之增大，则→0。
本节用有限积分法求解一个双调和方程，其数值结果与用有限差分法和有限元方法得出的数值结果相比较，通过实例验证方法是否有效、适用。
一般来说，数值结果的精确性是由其数值结果的误差大小决定。通常数值误差可描述为：绝对误差；相对误差
均方根误差Ermse
。
算例1　考虑定解问题(8)，令：
h(x,y)=8[3x2(1-x)2+3y2(1-y)2+(6x2-6x+1)(6y2-6y+1)],f(x,y)=g(x,y)=0，本文用有限积分法解该问题，并将计算得的数值结果与文献[1]中的十三点格式的有限差分法和广义有限差分法以及文献[3]中的有限元方法得出的数值结果进行对比。数值结果及精确解见表1，数值误差见表2。
由表1、表2易见该方法具有如下优点：该方法是一种边界型无网格法，具有不需要划分网格的优势；与其它数值解法相比较，该方法具有更高的计算精度；该方法的精度与节点数量的多少无关。因此有限积分法是处理双调和边值问题的一种有效数值解法。
算例2　在弹性力学[14]里，对于一受一般载荷的平面薄板，由薄板的小挠度弯曲理论可得该薄板弹性曲面的微分方程为：

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。