微波仿真论坛附录COMSOLMultiphysics的MATLAB矢量计算基础

1、附录 COMSOL Multiphysics的MATLAB矢量计算基础W. B. J. ZIMMERMAN1，J. M. REES21Department of Chemical and Process Engineering, University of Sheffield,Newcastle Street, Sheffield S1 3JD United Kingdom2Department of Applied Mathematics, University of Sheffield, Hicks Building, Sheffield矢量计算支撑了偏微分方程和它们的数值近似求解。为了很

2、好的使用有限元方法，建模人员应该掌握矢量计算基础知识。本科毕业的工程师可能学过矢量计算的数学课程，但是由于没有碰到过矢量计算的实际应用，这时在工程建模中使用矢量计算就受到限制。本附录介绍了所有COMSOL MULTIPHYSICS WITH MATLAB中用到的矢量计算基础知识。所以也可以将该附录当作是COMSOL MULTIPHYSICS WITH MATLAB多变量微分计算的入门读本。当我们写该附录时曾经争论过是否将这部分内容直接加入到第一章（数值分析基础）中，因为导数的数值近似是偏微分方程求解的基础，而偏微分方程是COMSOL MULTIPHYSICS的基本运算单元。确实，在学习波谱法求

3、解偏微分方程时，基本理论就是“导数理论”如何使用波变换方法来近似导数。所以通过对比发现，有限元方法的基础就是数值微分。所以争论就不存在了，第一章主要是关于COMSOL MULTIPHYSICS直接计算的基本问题的。但是不管多有用，近似导数仍然只是建模的一个中间步骤，不是目标本身。我们这里只考虑用于矢量计算的MATLAB基础，本附录的重点在于特征值分析和逻辑表达式。这些在整本书中都有体现。应当注意到我们这里介绍的每个功能都可以在COMSOL Script中实现。本书中唯一不能在COMSOL Script中实现的Matlab命令就是fminsearch。1矢量回顾1.1 矢量表达FEMLAB可以处

4、理标量、矢量和矩阵数据，这里简单介绍一下矢量的表达（作为MATLAB矩阵数据类型的一个特例）。标量可以作为一个单独的数，但是矢量是具有大小和方向的。在如图1所示的右手坐标系系统中，向量a用以下形式表达： (1)这里，和是坐标方向的单位矢量，是向量在各轴方向上的分量。它们是对各单位矢量，和的投影。对于坐标系中的P点（x，y，z），矢量P对于初始坐标系统O的位置为： (2)MATLAB用分量的形式描述列矢量或行矢量：>> a = 1; 2; 3; % column vector>> a = 1 2 3; % row vector在行向量中，空白（任意连续空格）作为分界符。列

5、向量用分号或者回车符分界：>> a = 123;1.2 内积，矩阵乘法，单位矢量和矢积典型的内积（或点积）定义为： (3)这里是矢量和的夹角。为了在MATALB中达到相同的目的，我们使用*运算符：>> a = 1; 2; 3;>> b = -3 2 -1;>> b*aans =-2这是一个行向量（1×3矩阵）乘以列向量（3×1矩阵）的特殊情况。因为前者的列数和后者的行数相等，这两个矩阵是相容的，可以根据矩阵乘法通用法则计算。 (4)如果A是m×n矩阵，B是n×l矩阵，则AB是m×l矩阵。如果共用的

6、维数不同，那么矩阵不相容，不能定义乘法运算。MATLAB也可以将标量乘法作为特殊矩阵乘法来计算，但是必须考虑矩阵的相容性。例如>> a*bans =-3 2 -1-6 4 -2-9 6 -3出现了什么情况？很简单，a是3×1矩阵，乘以1×3矩阵b，得到的ab是3×3矩阵。 (5)对于向量，矩阵(ab)ik称为a和b的并积或并矢。这是矩阵外积的特殊情况，这里标量乘积也算内积。在MATLAB中通过转置运算符“”可以实现两个行向量或两个列向量的内积，它是一个一元运算符，容易被误解为英语中的缩写符号。ans =-2but>> a*bans =-3

7、2 -1-6 4 -2-9 6 -3仍然产生并矢。必须自己考虑矩阵的相容性。如果a和b是行向量，那么*a或a*将产生内积还是外积？出于这个原因，MATLAB提供了一个特殊的dot函数来取消这种相容性的差别。>> help dotDOT Vector dot product.C = DOT(A,B) returns the scalar product of the vectors A and B.A and B must be vectors of the same length. When A and B are bothcolumn vectors, DOT(A,B) is t

8、he same as A*B.DOT(A,B), for N-D arrays A and B, returns the scalar productalong the first non-singleton dimension of A and B. A and B musthave the same size.DOT(A,B,DIM) returns the scalar product of A and B in thedimension DIM.See also CROSS.Example.>> dot(a,b)ans =-2>> dot(1; 2; 3,-3

9、2 -1)ans =-2使用dot函数使得行/列向量的组合不受影响。向量值向量的值或模可以通过下式计算： (6)MATLAB通过下式计算向量的模。ans =3.7417or with the built-in command norm>> norm(a,2)ans =3.7417这里sqrt( )是预置的平方根函数。单位向量单位向量是模为1的向量。单位向量可以通过归一化处理得到，即： (7)例如，>> ahat=a/norm(a,2)ahat =0.26730.53450.8018以上除法是标量除法，是矢量的每个元素除以标量。COMSOL Multiphysics通常根

10、据预置的几何单位向量来描述边界。nx，ny，nz通常作为边界上指向外部的单位矢量。较少使用tx和ty作为边界曲线的切线矢量。通常二维表面有两个正交切线方向，所以有无穷多个切向矢量。叉积叉积定义如下： (8)这里是置换张量，当指数ijk为123的正置换时取1，当是123的负置换时取1，其它情况均为零。是包含a和b的平面内的单位标准向量。是第i个坐标方向的单位向量。MATLAB提供了一个特殊的函数来计算叉积。>> help crossCROSS Vector cross product.C = CROSS(A,B) returns the cross product of the ve

11、ctorsA and B. That is, C = A x B. A and B must be 3 elementvectors.C = CROSS(A,B) returns the cross product of A and B along thefirst dimension of length 3.C = CROSS(A,B,DIM), where A and B are N-D arrays, returns the crossproduct of vectors in the dimension DIM of A and B. A and B musthave the same

12、 size, and both SIZE(A,DIM) and SIZE(B,DIM) must be 3.See also DOT.For example,>> cross(a,b)ans =-8-88>> cross(b,a)ans =88-8我们看到叉积中的因子阶数决定了叉积的正负号，等效于改变了单位标准向量的方向。1.3 数组：简单数组，元包数组和结构体数组处理本质上就是从FEMLAB中提取数据。为方便起见，FEMLAB对于多物理场采用fem结构体组织模型，对于扩展多物理场采用xfem结构体组织模型。从结构体和元包数组中提取有用信息是访问FEMLAB模型（和结果

13、）的一个非常有用的方法。简单数组数组有维数（m×n×l）。矩阵是二维数组。每个维数都有长度。所以size( )和length( )是两个非常有用的命令。>> a = 1 2 3 4; 5 6 7 8;>> size(a)ans =2 4数组的尺寸本身就是一个行向量，它的长度等于数组的维数。>> length(a(1,:)ans =4第二个自变量的冒号（:）表示包括第二维的整个范围，在本例中表示元素1到4。>> length(a(:,3)ans =2>> a(1,2:4)ans =2 3 4实际上冒号相当于低维数的子

14、数组。a(1,:)是第一行；a(:,3)是a的第三列。a(1,2:4)给出了第1行2到4列元素组成的子数组。对于更高维数情况，提取出的子数组更加复杂，例如：>>b=ones(2,2,2)b(:,:,1) =1 11 1b(:,:,2) =1 11 1>> b(1,:)ans =1 1 1 1前两种情况下子数组是矩阵，第三种情况下最终两个维数聚集成一个简单的行向量。FORTRAN编程人员可能感觉单个元素的寻址比子数组更舒服，可能会使用循环结构来达到相同的目的。>> a(1,3)ans =3构造数组数组可以通过冒号元算符自动生成，即>> a=0: 0

15、.1: 1*pia =Columns 1 through 80 0.3142 0.6283 0.9425 1.2566 1.5708 1.8850 2.1991Columns 9 through11 2.5133 2.8274 3.1416这样产生了从0到间隔相等的11个值。a=linspace(0,pi,11) a =Columns 1 through 80 0.3142 0.6283 0.9425 1.2566 1.5708 1.88502.1991Columns 9 through11 2.5133 2.8274 3.1416linspace是一个通用的矩阵自动生成命令，它起着旧编程语言

16、中循环结构的作用。有的时候也用得到logspace函数。>> help linspaceLINSPACE Linearly spaced vector.LINSPACE(X1, X2) generates a row vector of 100 linearlyequally spaced points between X1 and X2.LINSPACE(X1, X2, N) generates N points between X1 and X2.For N < 2, LINSPACE returns X2.See also LOGSPACE, :.>> he

17、lp logspaceLOGSPACE Logarithmically spaced vector.LOGSPACE(X1, X2) generates a row vector of 50 logarithmicallyequally spaced points between decades 10X1 and 10X2. If X2is pi, then the points are between 10X1 and pi.LOGSPACE(X1, X2, N) generates N points.For N < 2, LOGSPACE returns 10X2.See also

18、LINSPACE, :.四种其它常见数组生成函数为zeros，ones，rand和eye。zeros生成全0数组；ones生成全1数组；rand生成均匀分布的随机数组；eye生成单位矩阵。>> help zerosZEROS Zeros array.ZEROS(N) is an N-by-N matrix of zeros.ZEROS(M,N) or ZEROS(M,N) is an M-by-N matrix of zeros.ZEROS(M,N,P,.) or ZEROS(M N P .) is an M-by-N-by-P-by-.array of zeros.ZEROS(S

19、IZE(A) is the same size as A and all zeros.>> help onesONES Ones array.ONES(N) is an N-by-N matrix of ones.ONES(M,N) or ONES(M,N) is an M-by-N matrix of ones.ONES(M,N,P,.) or ONES(M N P .) is an M-by-N-by-P-by-.array of ones.ONES(SIZE(A) is the same size as A and all ones.>> help randRAN

20、D Uniformly distributed random numbers.RAND(N) is an N-by-N matrix with random entries, chosen froma uniform distribution on the interval (0.0,1.0).RAND(M,N) and RAND(M,N) are M-by-N matrices with random entries.RAND(M,N,P,.) or RAND(M,N,P,.) generate random arrays.RAND with no arguments is a scalar

21、 whose value changes each time itis referenced. RAND(SIZE(A) is the same size as A.>> help eyeEYE Identity matrix.EYE(N) is the N-by-N identity matrix.EYE(M,N) or EYE(M,N) is an M-by-N matrix with 1s onthe diagonal and zeros elsewhere.EYE(SIZE(A) is the same size as A.到目前为止你可能已经注意到，如果知道命令的名字，M

22、atlab就可以给出相关语法。但是你怎么才能知道命令的名字呢？在COMSOL Script中，只要输入help就会产生一列命令类型。随后输入“help type”命令就可以显示出相关命令。例如，“help numerics”给出了一列数值命令。然后“help daspk”给出调用daspk命令的语法，这是单独为FEMLAB准备的刚性DAE求解器。标量数组运算对数组的标量四则运算等于对数组中元素的四则运算，例如：>> 3*aans =Columns 1 through 80 0.9425 1.8850 2.8274 3.76994.7124 5.6549 6.5973Columns

23、9 through11 7.5398 8.4823 9.4248>> 5+aans =Columns 1 through 85.0000 5.3142 5.6283 5.9425 6.25666.5708 6.8850 7.1991Columns 9 through 117.5133 7.8274 8.1416数组数组的基元运算数组对数组的四则运算是一个比较偏的内容。如果数组大小一致，那么可以对元素采用点运算符：>> b=linspace(1,11,11)b =1 2 3 4 5 6 7 8 9 10 11>> size(a)ans =1 11>>

24、; size(b)ans =1 11>> a.*bans =Columns 1 through 80 0.6283 1.8850 3.7699 6.28329.4248 13.1947 17.5929Columns 9 through1122.6195 28.2743 34.5575元包数组和结构体关于元包数组和结构体可以写整整一章的内容。对于这两种数据结构，应该注意到它们是不同种类混合数据类型的容器包括浮点数，复数，矩阵，字符串和其它数组和结构体。元包数组将数组中的每个元包作为指针。元包数组通过包含数组基元的括弧来直接定义。元包命令将会产生一个空的构架，你可以在其中分别填入元素或

25、数组。>> ca = every, good, boy, does, find, 3+3i, 0 1; -2 2 ca =Columns 1 through 6every good boy does find 3.0000+ 3.0000iColumn 72x2 double通过圆括号内的数组地址可以进行调用，返回元包的元素。>> ca(3)ans =boy通过大括号内的地址数可以返回元包中基元包含的内容。>> ca3ans =boy可能通过矩阵元包更容易理解这一点。>> ca(7)ans =2x2 double>> ca7ans =

26、0 1-2 2结构体更多的借鉴了域的概念，也常称作枚举法，在C语言中比较常见。使用结构体和元包数组的最大区别在于如果你改变结构体，增加或减少域并不会改变域的次序。但是减少或增加元包数组元素会改变元包编号，或在数组中产生“空穴”。COMSOL Multiphysics用户最常见到的结构体就是fem结构体，这是COMSOL Multiphysics对多物理场模型组织数据和计算结果的结构体。从COMSOL Multiphysics中导出fem结构体到MATLAB工作空间中产生以下电动流动模型。>> femfem =version: 1x1 structappl: 1x1 struct 1

27、x1 struct 1x1 structgeom: 1x1 geom2mesh: 1x1 femmeshshape: 1x7 cellgporder: 4 2cporder: 2 1simplify: offborder: 1outform: weakform: weakequ: 1x1 structbnd: 1x1 structpnt: 1x1 structexpr: zeta -zeta1*(Y+zetar*(1-Y) sig Y+sigr*(1-Y)elemcpl: 1x1 structdraw: 1x1 structconst: 1x34 cellglobalexpr: 1x12 ce

28、llxmesh: 1x1 com.femlab.xmesh.Xmeshsol: 1x1 femsol以上fem结构体的域列表给出了对各域中内容的描述。每个域都可以通过“点”号来寻址。>> fem.exprans =zeta 1x22 char sig Y+sigr*(1-Y)fem.expr是有四个元包的元包数组；该元包数组非常小，能够显示所有的内容。由于是个元包数组，可以通过大括号来调用元包中的内容。>> fem.expr2ans =-zeta1*(Y+zetar*(1-Y)正如我们看到的，fem结构体中含有元包数组，其它结构体，字符和数字。我们也可以用fem结构体来

29、构成元包数组，这就是COMSOL Multiphysics处理扩展多物理场时构建的xfem结构体，其中每个fem结构体对应一个逻辑几何结构。我们通常需要对fem和xfem结构体调用postinterp命令，来插值取得有限元离散中间点处的变量计算值：is,pe=postinterp(xfem,xx);u=postinterp(xfem,u1,is);通过对xfem结构体调用postinterp函数可以得到对模型和解的完整描述，对于不同的结构体会执行不同的命令分支。作为COMSOL Mulitphysics用户，有必要知道COMSOL Multiphysics模型和解中用到的MATLAB数据结构，

30、这样当我们希望进行特殊后处理或建模（而不是在COMSOL Multiphysics图形用户界面中实现）时，可以从中提取相应的数据。2标量和矢量场：MATLAB函数的表达物质属性通常取决于不同的位置或时间。在人们可见的长度尺度上，大多数物理量都看作连续函数在每个数值点都有相应值。这些量就被称作场。像温度和压力等描述单个值的量就是标量场。标量场是一个单一的数值，即 (9)三维矢量场需要三个分量： (10)F的每个分量都是位置的标量函数。例：图1 对于不同C值的等高线图1给出了常数C在不同值时的等高线图。MATLAB中没有数据类型来描述场量，这些量可以以函数形式来表达。在MATLAB中有三种主要的函

31、数表达方式：(1) inline函数，只在工作空间中定义>> myfun = 1+log(r);>> myfuni=inline(myfun,r)myfuni =Inline function:myfuni(r) = 1+log(r)>> a=feval(myfuni,1)a =1>> a=feval(myfuni,10)a =3.3026在COMSOL Multiphysics的Options菜单下有个新功能专门用来以公式形式定义inline函数。(2) m文件函数，它以文件形式存储在硬盘中，可以被工作空间或m文件代码调用。例如m文件函数tem

32、perat.m包含以下代码，保存在MATLAB当前文件夹中。function t=temperat(r)%TEMPERAT evaluates T = 1 + ln r% T = temperat(r)%t=1+log(r);>>a=temperat(10)a =3.3026 (3) 内插函数。有时候确定了函数在某些特定点的值，而在附近一些点的值就需要假设函数局部光滑来计算得到。内插需要一系列的MATLAB函数，但是最终产生一种函数形式。对于一维、二维和三维数据，MATLAB有预置的插值函数，分别称作interp1，interp2和interp3。COMSOL Multiphysi

33、cs的Options菜单下有个新功能以手动输入数据或给定格式的数据文件来定义内插函数。第二章122页举了个将水密度分布作为温度函数输入的例子。通常情况下，COMSOL Multiphysics中系数和边界条件等场变量都是通过预定义独立变量、因变量、派生变量串联方式输入的。例如，在通用模式下，可以输入单一因变量u和独立变量x的表达式：但是MATLAB m文件函数可以很容易的调用。这里要指出的很重要一点是，COMSOL Multiphysics通常以标量分量形式接受输入数据。如果需要一个向量或矩阵，就需要通过对标量分量形式输入，其中每个分量都可以是复杂函数。另外很重要的一点是，COMSOL Mul

34、tiphysics允许函数有矢量输入，即每个函数的自变量都可以是矢量。所以你在m文件函数中使用显式元素运算符，例如“.*”时，必须非常小心。COMSOL Multiphysics用fem结构体来记录计算结果，对于fem.mesh中定义的网格，自由度保存在fem.sol中。对于每个因变量和派生变量，COMSOL Multiphysics提供了一个特殊的内插函数来从fem.sol中取得插值。本书例题中大量使用postinterp内插函数。它甚至可以在m文件函数中自动调用mat文件中的fem结构体。3多变量计算中的微分3.1 标量场梯度如果，那么矢量 (11)被称作标量场的梯度，通常表示为grad。

35、梯度运算符是三维直角坐标系中的矢量运算符。 (12)COMSOL Multiphysics举例假设，则。但是MATLAB不直接处理这类符号计算（它的符号工具箱可以）。COMSOL Multiphysics粗略计算微分的数值近似解。所以标量场的梯度可以通过COMSOL Multiphysics的“基元”运算符得到。我们是如何轻易得到这一信息的呢？表1中给出了具体做法。应当注意到，由于实际上没有偏微分方程被求解，Neumann边界条件等价于中性或无条件边界。否则的话，只有边界数据满足0=phi-x2-y2时才是可能的解。表1 梯度模型模型导航栏二维，PDE模式，通用模式，独立变量：x，y 因变量：

36、phi选项设置Axes/Grid为-1，1×-1，1绘图矩形域-1，1×-1，1边界模式/边界设定设定四个边界均为Neumann边界条件子域模式/子域设定在域1中，设定0 0；da0；Fphi-x2-y2网格模式重绘默认网格（418节点，774基元）求解使用默认设定（非线性求解器）后处理选择arrow模式，如图A4现在将fem结构导入MATLAB中。我们通过MATLAB双线性回归插值得到近似数值解。>>x=0.5;y=0.5;xx,yy=meshgrid(-1:0.01:1,-1:0.01:1);xxx=xx(:); yy(:);phix=postinterp(

37、fem,phix,xxx);phiy=postinterp(fem,phiy,xxx);uu=reshape(phix,size(xx);vv=reshape(phiy,size(xx);u=interp2(xx,yy,uu,x,y);v=interp2(xx,yy,vv,x,y);u,vans = 1.0000 1.0000所有reshape命令都是确保数据满足interp2的正确格式。表2给出了其它点的梯度计算。表2 使用梯度模型对grad的数值计算Xyphixphiy0.50.51.00001.0000-0.250.75-0.50001.50000.75-0.51.5000-1.0000

38、0.25-0.750.5000-1.5000通过任意解法，FEM得到的一阶导数都十分准确。根据收敛准则，全局误差为O(10-16)。图2给出了梯度的数值近似解。图2 grad的矢量箭头图3.2 方向导数的方向导数就是沿给定方向的变化速率。如果是该方向的单位向量，则方向导数可以表示为： (13)坐标方向非常容易求解，即： (14)在第七章的ECT模型中我们使用方向导数计算电势的法向导数。显然，方向导数与通量的概念密切相关。对于“线性”特性系统（Fick定律，Fourier定律等），穿过材料表面的总通量正比于表面法向导数的积分。局部通量正比于某一点的法向导数。对于这一点，大多数矢量计算文章都已经证

39、明，的变化速率最大的方向就是grad的方向，且|grad|就是在该方向的变化速率。对于点(x,y)=(0.25,-0.75)，通过沿角度的步进法绘制，就能清楚的看到这一点。相关MATLAB代码如下：>> theta=linspace(0, pi, 100);dirder = zeros(size(theta);for k=1:length(theta)dirder(k)=cos(theta(k)*u+sin(theta(k)*v;endplot(theta, dirder)出于我的FORTRAN编程习惯，以上代码使用了循环结构。如果是对MATLAB编程比较熟悉，一般不使用循环法：&

40、gt;>dirder = u*cos(theta)+v*sin(theta);plot(theta,dirder)cos和sin函数会对向量thera中的每个元素进行计算，生成相同长度的向量。图3 方向角度随弧度变化图。注意存在最小的方向导数下降最陡的方向，对应于梯度方向。3.3 水平集合/水准面注意到方向导数穿过x轴，即某一个方向的方向导数为零在该方向没有任何变化速率。可以看出的方向垂直于梯度方向。所以在该方向局部为常数。绘制出每个为常数的曲线（二维）或曲面（三维），就会形成一簇曲线（曲面），称作的水平集合（见第八章）。在二维空间中，水平几何也被称作等高线。方向导数的概念非常类似于测量

41、图，就是陆地的海拔。等高线具有相同的海拔；方向导数是在方向的爬山速率，梯度就是最陡的爬山方向，爬山速率为|grad|。在流体动力学中，通常用流函数的等高线来描述，等高线就是正切于速度场的流线（稳定流动中的颗粒路径）。在第三章浮力对流的例子中介绍了如何计算流函数（见第三章方程(3)）。3.4 矢量场导数矢量微分运算符可以通过两种方式作用于矢量场：(1)内积形式，称作散度；(2)叉积形式，称作旋度。散度为： (15)旋度为： (16)是前面引入的置换张量。很容易会发现散度是标量，而旋度是矢量。运算符经常在传热或传质方程的对流项中见到。它显然不是散度，因为 (17)和方程(15)的标量相比较，它仍然

42、是个运算符。正如我们前面看到的，导数的数值近似是COMSOL Multiphysics的基本功能，所以我们可以用来近似计算散度和旋度。COMSOL Multiphysics举例假设，表3中给出具体计算步骤。表3 求解散度和旋度步骤模型导航栏三维，PDE模式，通用模式，独立变量：x，y，z因变量：u1，u2，u3选项设置Axes/Grid为0，1×0，1×0，1绘图绘制六面体BLK10，1×0，1×0，1边界模式/边界设定设定所有边界为Neumann边界条件子域模式/子域设定设定0 0 0; da1=0 0 0; F1=u1-x2设定0 0 0; da2=

43、0 0 0; F2=u2-3*x*y设定0 0 0; da3=0 0 0; F3=u3-x3网格模式设定网格缩放因子为3，重绘网格（201节点，719自由度）求解使用默认设定（非线性求解器）后处理(1)绘制散度u1x+u2y+u3z的彩色图(2)绘制旋度(u3y-u2z,u1z-u3x,u2x-u1y)的箭头图应当再次注意到，实际上没有求解任何偏微分方程，所以Neumann边界条件等价于中性和无条件边界。否则的话，只有边界数据满足时才是可能的解。通过符号计算很容易算得：数值解近似程度如何？尝试以下散度：>> xxx=0.42; 0.57; 0.33;postinterp(fem,u1x+u2y+u3z,xxx)ans = 2.1137>> 5*0.42ans = 2.1000显