两点边值问题的三类边值条件的有限元解法实现

卢仁洋 于陆洋 江山

摘 要：文章研究二阶微分方程的两点边值问题，使用有限元方法对三类不同的边值条件具体进行讨论和处理。对于可齐次化的Dirichlet、Neumann边值， 给出相应分析以简化和规范计算步骤。对于Robin边值， 基于之前的分析给出实现技巧以达到有效的数值模拟。

关键词：有限元方法；Dirichlet边值；Neumann边值；Robin边值

中图分类号：O172 文献标志码：A 文章编号：2096-000X（2018）05-0058-03

Abstract： A two-point boundary value problem of the second-order ordinary differential equation is studied in this paper， and a finite element method is used to deal with three kinds of boundary conditions in detail. The corresponding analysis is provided to simplify and standardize the steps for the homogeneous Dirichlet， Neumann boundary values. For the Robin boundary value， we utilize the previous analyses to present a specific strategy for a good performance in numerical simulations.

Keywords： finite element method； Dirichlet boundary； Neumann boundary； Robin boundary

一、本文模型及介绍

二阶微分方程的两点边值问题是科学工程计算中的经典问题，也是微分数值解法的必要基础[1，2]。常见的做法为有限差分法[3，4]、有限元方法[5-7]、有限体积法[8，9]等。

本文考虑两点边值问题（BVP）

其中变系数 。

上述区间I也可以是I=[a，b]，但为研究方便，能通过变换x'=转化为I=[0，1]。

针对上述方程，根据A1，A2，B1，B2取值的不同，分成三类边值条件：

1. Dirichlet边值：A1≠0，B1≠0，A2=0，B2=0，即u（0）=？琢/A1，u（1）=？茁/B1。

2. Neumann边值：A1=0，B1=0，A2≠0，B2≠0，即u'（0）=？琢/A2，u（1）=？茁/B2。

3. Robin边值：A1≠0，B1≠0，A2≠0，B2≠0。

本文采用有限元法（Finite Element Method， FEM）處理模型，首先需确定变分形式。对（1.1）移项后Lu-f=0，两边同乘v并在[0，1]上积分，得：

即有： （1.2）

其中a（u，v）=（p+qv+ruv）dx，（f，v）=fvdx，H1为一阶可导的完全内积空间。设Uh是U的n维子空间，选取Un的一组基底？渍1，？渍2…，？渍n作为基函数，有限元解即为uh=u？渍，取线性元基函数？渍1，？渍2…，？渍n计算。通过||u-uh||∞计算误差、log2（）计算收敛率，其中u真解、uh为有限元解。

二、三类边值问题的分析和算例

根据边值条件是否齐次，又可具体细分。

（一）Dirichlet边值

（1）？琢=？茁=0，即u（0）=0，u（1）=0

（2）？琢=0，？茁≠0，即u（0）=0，u（1）=？茁/B1

处理（2）的做法：令w（x）=u（x）-u（1），则有：

即 ，就转化为（1）的情形。

（3）？琢≠0，？茁=0，即u（0）≠？琢/A1=？琢D，u（1）=0。可以通过令：（x）=u（x）+u（0），后续处理与（2）一致。

（4）？琢≠0，？茁≠0，即u（0）≠？琢/A1=？琢D，u（1）=？茁/B1=？茁D，可齐次化左边界，令w（x）=u（x）-u（0），则有：

即 ，进而转化为（2）。

算例1：式（1.1）中，取p（x）=x q（x）=3x，r（x）=-2x，f（x）=-2e2x-ex。

A1=1，A2=0，？琢=1，B1=1，B2=0，？茁=e2+e精确解为u=e2x+xex。

通过上述分析，令w（x）=u（x）-u（0）-（u（1）-u（0））=u（x）-（e2+e-1）x-1， 则有

经典有限元法求出wh后，得到uh=wh+（e2+e-1）x+1。 计算结果如下图表：

可以看出，在N=32较粗网格即有10-3级误差，随着网格不断剖分、有限元法的逼近效果越来越好，并且收敛率达到稳定的2阶收敛，这是符合预期结果的。

（二）Neumann边值

（1）？琢=？茁=0，即u'（0）=0，u'（1）=0。

（2）？琢=0，？茁≠0，即u'（0）=0，u'（1）=？茁/B2=？茁N。

（3）？琢≠0，？茁=0，即u'（0）≠？琢/A2=？琢N，u'（1）=0。

（4）？琢≠0，？茁≠0，即u'（0）≠？琢/A2=？琢N，u'（1）=？茁/B2=？茁N。

通过与Dirichlet边值类似的处理方法，先齐次化左边界，再齐次化右边界。具体做法：令w（x）=u（x）-u'（0）x-■ （u'（1）-u'（0））x2

即有

算例2：式（1.1）中，取p（x）=x q（x）=3x，r（x）=-2x，f（x）=-2e2x-ex。

A1=0，A2=1，？琢=3，B1=0，B2=1，？茁=2e2+2e该方程的精确解为u=e2x+xex。

由上述分析，令w（x）=u（x）-3x-■（2e2+2e-3）x2，原问题转化为

经典有限元法求出vh后，得到uh=wh+■（2e2+2e-3）x2+3x。

表2 Neumann边值问题的误差与收敛率

（三）Robin边值

u'（0）=（？琢-A1u（0））/A2=？琢R1-？琢R2u（0），u'（1）=（？茁-B1u（1）））/B2=？茁R1-？茁R2u（1）

由于Robin边值特性，导数项与函数项混合，无法使用齐次化方法。这里给出一般化的解决方法：在形成总刚度矩阵以及右端载荷向量时，根据边值条件相对应地修改矩阵和向量，从线性元性质出发，实际需要修改的是矩阵的第1×1项和第（N+1）×（N+1）项以及载荷向量的第1项和第N+1项。

算例3：式（1.1）中，取p（x）=x，b（x）=3x，c（x）=-2x，f（x）=-2e2x-ex。

A1=1，A2=1，？琢=4，B1=-2，B2=1，？茁=0該方程的精确解为 u=e2x+xex。

利用经典有限元法形成总刚度矩阵A以及载荷向量 F，构成了（N+1）×（N+1）阶数，将Robin边值的导数项用函数项表示并代入（1.2）式得：

表3 Robin边值问题的误差与收敛率

通过详细分析和具体验证，我们验证了三类边值中无论哪类模型，使用有限元法和数学技巧处理都得到了有效精度和稳定收敛的2阶结果。本文是作者学习过程中对于有限元知识的一些总结和融合，将三类边值问题做了具体的分析以及数值案例的展示，也对线性有限元的误差进行了实际验证。

参考文献：

[1]李荣华.偏微分方程数值解法[M].北京：高等教育出版社，2012.

[2]陆金甫.偏微分方程数值解法[M].北京：清华大学出版社，2016.

[3]许芝卉，王凤艳.用差分方法求解两点边值问题[J].山西大同大学学报（自然科学版），2009，25（4）：10-11.

[4]赵秉新，郑来运.两点边值问题的差分方法及快速求解[J].华南师范大学学报（自然科学版），2010（4）：27-30.

[5]Huang Z Y. Tailored finite point method for the interface problem[J].Networks and Heterogeneous Media，2009，4（1）：91-106.

[6]付小龙，李名，李小瑞.两点边值问题的有限元算法[C].北京力学会学术年会，2014.

[7]Constantinou P， Xenophontos C. Finite element analysis of an exponentially graded mesh for singularly perturbed problems[J]. Computational Methods in Applied Mathematics，2015，15（2）： 135-143.

[8]Cao W， Zhang Z， Zou Q. Superconvergence of any order finite volume schemes for 1D general elliptic equations[J].Journal of Scientific Computing，2013，56（3）：566-590.

[9]魏保军，张武军，石金娥.两点边值问题有限体积法的一种加权模估计[J].广西师范大学学报（自然科学版），2015，33（3）：75-78.