有限元数值解法在MATLAB中的实现及可视化
摘 要：偏微分方程的数值解法在数值分析中占有很重要的地位，很多科学技术问题的数值计算包括了偏微分方程的数值解问题。在学习初等函数时，总是先画出它们的图形，因为图形能帮助我们了解函数的性质。而对于偏微分方程，画出它们的图形并不容易，尤其是没有解析解的偏微分方程，画图就显得更加不容易了。为了从偏微分方程的数学表达式中看出其所表达的图形、函数值与自变量之间的关系，通过MATLAB编程，用有限元数值解法求解了偏微分方程，并将其结果可视化。
关键词：偏微分方程；MATLAB；有限元法；可视化
1 引言（Introduction）
偏微分方程的数值解法在数值分析中占有很重要的地位，很多科学技术问题的数值计算包括了偏微分方程的数值解问题。近三十多年来，它的理论和方法都有了很大的发展，而且在各个科学技术的领域中应用也愈来愈广泛。例如，核武器的研制要有理论设计和核试验。但核反应和核爆炸的过程是在高温高压的条件下进行的，而且巨大的能量在极短的时间内释放出来，核装置内部的细致反应过程及各个物理量的变化是根本不能用仪器测量出来的，核试验只是提供综合的数据。而描述核反应和爆炸物理过程的数学模型是一个很复杂的非线性偏微分方程组，也根本没有办法得到这个方程组理论上的精确解。所以发展核武器的国家都在计算机上对核反应过程进行数值模拟，这也称为“数值核实验”，它可以大大减少核试验的次数，节约大量的经费，缩短研制的周期[1]。
在学习初等函数时，总是先画出它们的图形，因为图形能帮助我们了解函数的性质。而对于偏微分方程，画出它们的图形并不容易，尤其是没有解析解的偏微分方程，画图就显得更加不容易了。所以本文主要研究如何用MATLAB数值求解偏微分方程，并将其数值解绘制成三维图形的形式，从而可以从复杂的数学表达式中看出其所表达的图像、函数值与自变量之间的关系[2]。
2 有限元法（Finite element method）
2.1 有限元法概述
有限元法的基本思想是将结构离散化，用有限个容易分析的单元来表示复杂的对象，单元之间通过有限个节点相互连接，然后根据变形协调条件综合求解。由于单元的数目是有限的，节点的数目也是有限的，所以称为有限元法。
一般来说，用差分法解偏微分方程，解得的结果就是方程的准确解函数在节点上的近似值。而用变分近似方法求解，是将近似解表示成有限维子空间中基函数的线性组合。在古典变分方法中，这样的基函数一般采用幂函数和三角函数等初等函数，又要求在区域的边界上满足边界条件，如果是二维或三维的不规则区域，这样的基函数往往很难构造出来。所以，古典的变分方法虽然是得到近似的解析解（与差分方法不同），但是对一般的区域，却往往难以实现。有限元方法，也是基于变分原理，由于选择了特殊的基函数，使它能适用于较一般的区域。这种基函数是区域的剖分有关的，近似解u表示为基函数的线性组合，而线性组合中的系数，又是剖分节点上u或其导数的近似解。所以有限元方法既是基于变分原理，又具有差分方法的一些特点，并且适合于较复杂的区域和不同粗细的网格。正是由于这些特点，20世纪60年代以来，有限元方法的理论和应用得到迅速的发展，适用范围也愈来愈广泛。
2.2 有限元法的基本思想
有限元法是用有限个单元将连续体离散化，通过对有限个单元作分片插值求解各种力学、物理问题的一种数值方法。有限元法把连续体离散成有限个单元：杆系结构的单元是每一个杆件；连续体的单元是各种形状（如三角形、四边形、六面体等）的单元体。每个单元的场函数是只包含有限个待定节点参量的简单场函数，这些单元场函数的集合就能近似代表整个连续体的场函数。根据能量方程或加权残量方程可建立有限个待定参量的代数方程组，求解此离散方程组就得到有限元法的数值解。
有限元方法的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解[3]。采用不同的权函数和插值函数形式，便构成不同的有限元方法。
2.3 有限元法的计算格式
根据所采用的权函数和插值函数的不同，有限元方法也分为多种计算格式。从权函数的选择来说，有配置法、矩量法、最小二乘法和伽辽金法；从计算单元网格的形状来划分，有三角形网格、四边形网格和多边形网格；从插值函数的精度来划分，又分为线性插值函数和高次插值函数等。不同的组合同样构成不同的有限元计算格式。插值函数一般由不同次幂的多项式组成，但也有采用三角函数或指数函数组成的乘积表示，但最常用的多项式插值函数。有限元插值函数分为两大类，一类只要求插值多项式本身在插值点取已知值，称为拉格朗日（Lagrange）多项式插值；另一种不仅要求插值多项式本身，还要求它的导数值在插值点取已知值，称为哈密特（Hermite）多项式插值。在二维有限元中，三角形单元应用的最早，近来四边形等参元的应用也越来越广[4]。对于二维三角形和四边形电源单元，常采用的插值函数为有Lagrange插值直角坐标系中的线性插值函数及二阶或更高阶插值函数、面积坐标系中的线性插值函数、二阶或更高阶插值函数等。
2.4 有限元法的解题步骤
对于有限元方法，其基本思路和解题步骤可归纳为：
（1）建立积分方程：根据变分原理或方程余量与权函数正交化原理，建立与微分方程初边值问题等价的积分表达式，这是有限元法的出发点。
（2）区域单元剖分：根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相互连接、不重叠的单元[5]。区域单元划分是采用有限元方法的前期准备工作，这部分工作量比较大，除了给计算单元和节点进行编号和确定相互之间的关系之外，还要表示节点的位置坐标，同时还需要列出自然边界和本质边界的节点序号和相应的边界值。
（3）确定单元基函数：根据单元中节点数目及对近似解精度的要求，选择满足一定插值条件的插值函数作为单元基函数。有限元方法中的基函数是在单元中选取的，由于各单元具有规则的几何形状，在选取基函数时可遵循一定的法则。
（4）单元分析：将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近；再将近似函数代入积分方程，并对单元区域进行积分，可获得含有待定系数（即单元中各节点的参数值）的代数方程组，称为单元有限元方程。
（5）总体合成：在得出单元有限元方程之后，将区域中所有单元有限元方程按一定法则进行累加，形成总体有限元方程。
（6）边界条件的处理：一般边界条件有三种形式，分为本质边界条件（狄里克雷边界条件）、自然边界条件（黎曼边界条件）、混合边界条件（柯西边界条件）。对于自然边界条件，一般在积分表达式中可自动得到满足。对于本质边界条件和混合边界条件，需按一定法则对总体有限元方程进行修正满足。
（7）解有限元方程：根据边界条件修正的总体有限元方程组，是含所有待定未知量的封闭方程组，采用适当的数值计算方法求解，可求得各节点的函数值[6]。
3 有限元法在MATLAB中的实现（Reality of finite
element numerical method in MATLAB）
在用有限元法时，MATLAB编程的方法与其他语言如FORTRAN、C语言相似，但用MATLAB编程会更简单一些。以下详细介绍了用有限元方法计算十字形截面的方形同轴电缆线内的电势，该同轴电缆线内的电势值在每一点都不一样，无法用解析解来描述，所以用了有限元的方法解数值解，最后画出了等势线，对程序稍加修改，也可以用表面图反映数值解的分布。
（1）问题描述
有横截面如图1所示的电缆芯线，电缆的外导体的电势为零，芯线的电势为10V。求该同轴电缆线内的电势[7]。
（2）定解问题
这可以看成是二维问题，由于问题的对称性，只要解第一象限（四分之一区域）就可以，根据问题建立如下的方程组：
（3）区域单元剖分
根据本问题的性质，将区域划分为三角形，网格的划分如图2所示。
4 结论（Conclusion）
科学计算在各门自然科学（物理学、气象学、地质学和生命科学等）和技术科学与工程科学（核技术、石油勘探、航空航天和大型土木工程等）中起着越来越重要的作用，在很多重要领域中成为不可缺少的工具。而科学与工程计算中最重要的内容就是求解在科学研究和工程技术中出现的各种各样的偏微分方程或方程组。解偏微分方程已经成为科学与工程计算的核心内容，包括一些大型的计算和很多已经成为常规的计算。原则上，可以用FORTRAN或C语言来完成这些计算，但很少有人这样去做，原因是成本太高，编程太复杂。而MATLAB是一种用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境[8]。所以使用MATLAB，可以较使用传统的编程语言（如C、C++和Fortran）更快地解决技术计算问题。