使用 NAG MATLAB^® 工具箱 - Part 4

NAG MATLAB  工具箱能够让使用者在 MATLAB 中调用所有的 NAG 函数， 本文是这一系列的最后一篇文章，再次演示工具箱中的各种功能 (先前的文章可参考 Part 1、Part 2 与 Part 3)。 我们以前的文章中，都是单独使用工具箱的函数进行说明，在这篇的内容中， 我们会讨论一个较大的应用范例，同时会使用到工具箱中的许多函数，并且会采用不同的领域中的函数； 本应用程序在最新版本的工具箱中也有提供。

背景与简介

我们将介绍一个大量采用 NAG 工具箱函数，在二维空间中求解偏微分方程的应用程序。 偏微分方程通常是用来建构各种不同科学现象的数学模型 – 例如：热或声音的传导现象、电磁场的作用、流力以及固体变形等问题。 我们的应用程序是要用来求解所谓的 Poisson PDE，用来找出对应于给定的电荷分布的电位能； 在我们的程序中，用户可以交互式的定义默认函数的形式；同时也能求解拉普拉斯方程 (Poisson 偏微分方程右侧为 0)；它描述了给定介质中热散布的方式。

我们采用简单的有限元法求解偏微分方程的数值解，要被求解的偏微分方程区域可被分割成网状， 且偏微分方程可被转换成联立线性方程，并利用先大型稀疏系统算法求解。 利用有限元方法的优点是，它非常适合任一形状的区域；我们将在这个范例中说明，另外，我们使用 MATLAB 的 GUI 接口来输入与显示求解的区域。

我们的应用程序中采用了许多 NAG 工具箱的网格生成与求解函数。我们会在以下的部分详细说明。

D06 函数 – 网格生成

NAG 算法库中的 D06 函数可在用户定义的区域中，自动生成网格。它提供区域边界网格的建立， 并可在范围内提供三个不同的网格生成函数，还有其他相关的函数可用来提高网格的品质。

我们的应用程序能够提供用户自行透过一系列数据点所建构的多边形，或选择正规的矩形或圆形的方式，建立自己的区域边界； 边界可以是任意形状的，而且区域间可以包含其他内部边界或孔洞 (如图 5 所示)。 当建构好边界后，应用系统便会调用 d06ba 函数，将区域进行切割并生成边界的网格。

图 1 是我们应用程序的截图，用户以圆形的形态建构区域的边界 (我们可以观察到这是一个 11 边形)。

图 1：应用程序画面，图中显示 (左边绘图区) 区域的边界已经进行网格处理。控制边界的相关设定选单在上方， 右边的区域提供用户自行编辑边界条件，指定网格算法以及求解偏微分方程。右下角按钮提供说明文件。

已经有了网状的边界后，下一步便是生成区域的内部网格。我们前面提到，NAG 工具箱包含了三个不同的网格生成算法： Incremental (d06aa)、Delaunay-Voronoi (d06ab) 以及 Advancing Front (d06ac)。 每一个函数的参数都能够调整网格的粗糙程度 (如图 2)。 虽然在我们的应用程序中没有提供，但是在使用 Delaunay-Voronoi 与 Advancing Front 算法时，可以指定内部的许多点来调整精细网格程度； 这样的功能在某些问题上，是一个很好的解决方案，例如确定飞机机翼尾部的乱流。图 2-4 显示我们例子中使用三个算法的内部网格生成结果。
图 2：使用 Incremental 算法产生内部网格的三个例子。这算法中的参数可以控制内部三角形的大小；比较左图与中间的差异，中间的图中我们将此参数值减少。另一个参数可控制边界点间的距离；您可以比较右图与中间的差异，右图减少此参数的数值。

图 3：使用 Delaunay-Voronoi 算法产生内部网格

图 4：使用 Advancing front 算法产生内部网格

F04 函数 - 联立线性方程式

这些内部网格生成算法可以处理各种形状的区域，但是其边界必须明确的定义出来。特别是边界的边缘不能够交叉。 在我们例子中，我们将数据点输入后便调用 NAG 工具箱中的 f04aa 函数求解联立线性方程式。

以下是程序片段：

function [cross] = simeqn(coor, grad, i, j, x, y)

% (x,y) 是新边界线坐标
% coor 变量为所有边界线的坐标
% grad 变量为所有边界的斜率
% i 与 j 变量为求解之索引值

% 输入变量
% a 为方程式之 jacobian
% b 为包含边界线方程之矩阵
a = [-grad(i-1), 1; -grad(j), 1];
b = [y - grad(i-1)*x; coor(2,j) - grad(j)*coor(1,j)];

% 调用 NAG 函数求解方程
[a, c, ifail] = f04aa(a, b);

% 自变数 'c' 中，将两条线之交叉点坐标拿出
xx = c(1);
yy = c(2);

% 接下来的程序要确定以上的多边形边界方程解是否交叉

ab = norm(coor(:,j) - [xx;yy]) + norm(coor(:,j+1) - [xx;yy]);
distab = norm(coor(:,j) - coor(:,j+1));
AB = norm([x;y] - [xx;yy]) + norm(coor(:,i-1) - [xx;yy]);
distAB = norm([x;y] - coor(:,i-1));

if ab < distab + 1e-10 && AB < distAB + 1e-10
   % 多边形之线段交叉
    cross = 1;
else
    % 线未交叉
    cross = 0;
end

end

这是一个调用 F04 函数的简单例子， 透过功能强大的算法求解各种线性代数问题。具体的说，矩阵的系数可能是实数、复数、对称、Hermitian、正定、Toeplitz 正定或带状；矩阵也可以是矩形的，这种情况下可以求得最小二乘解。 F04 函数包含了两大类的函数：一个是黑盒子方式 (容易使用)；另一个则为通用型式 (提供更为弹性的功能，例如可同时提供右边矩阵求解)。

F04 函数也可以应用在 稀疏 (大部分的元素都为 0) 矩阵的计算中，然而在这样的情况下，我们会建议以 F11 函数求解，我们在以下部分讨论。

F11 函数 - 大尺度线性系统

我们在设定的区域中产生好了网格后，偏微分方程将会转化为联立的大型线性方程，系数的数量反应的是网格的大小 – 通常有数千阶。然而，系数的矩阵通常也是 带状的 (只有矩阵中对角线及少数非对角线的元素非 0)。此系统是属于稀方程组，可以使用 F11 函数求解，它只储存非 0 的元素。

偏微分方程求解最重要的先决条件就是边界条件，此会限制区域边界上的值 (或者是导数)。 在我们这个应用系统中，用户可以设定结果值 (一般称 Dirichlet 边界条件) 或导数 (Neumann 边界条件)。 边界条件的设定在前面程序 F11 函数的调用中可以设定。这里的例子是在区域中，透过拉普拉斯方程求解温度分布值，我们使用 Dirichlet 边界条件，在某些边界指定温度，并提供整个区域中的默认温度值 (图 5)。

F11 函数包含了许多求解大型稀疏线性系统的算法。如同 F04 的函数，他们也分成黑盒子 (易于使用) 以及通用 (较弹性) 两大类。 F11 同时也包含了许多公用的函数可以帮协助计算。 在我们的应用系统中，我们先调用 f11zb 公用函数，重新排定矩阵中非 0 的元素，这能提高计算效率。 再来调用 f11ja 与 f11jc 函数；第一个函数对矩阵做不完全 Cholesky 分解，同样的，这样的步骤可以提高计算的性能与速度。 f11jc 函数会调用其他 5 个 F11 函数，利用共轭梯度法 (如我们的例子) 或 Lanczosz 法，来求解实数对称稀疏线性系统。

经由计算，我们的应用系统可透果 MATLAB 的绘图工具，以动态方式拨放温度演变的过程。图 5 显示在内部有不规则孔洞，其周边默认值温度较高，其余部分温度较低的温度分布最终 (稳态) 计算结果。

图 5：拉普拉斯方程解 范围中包含 (不规则形状) 坑洞。内部边界温度为 100，边界外缘没有边界条件。

结论

本文中的应用是在 NAG MATLAB  工具箱最新版 22 版中的范例。在 Help 中 Toolboxes 内 (NAG demos: Applications 的 Mesh Generation 目录中。请注意：本系列中的文章可参考 这里； 同时也可在 NAG demos: MATLAB GUI 目录中找到。

NAG MATLAB 工具箱提供 32 位与 64 位的 Linux 与 Microsoft Windows 版本。其他细节请参考 http://cn.nag-gc.com/numeric/MB/start.asp。