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

为了很好的使用有限元方法，建模人员应该掌握矢量计算基础知识。

本科毕业的工程师可能学过矢量计算的数学课程，但是由于没有碰到过矢量计算的实际应用，这时在工程建模中使用矢量计算就受到限制。

本附录介绍了所有COMSOL MULTIPHYSICS WITH MATLAB 中用到的矢量计算基础知识。

所以也可以将该附录当作是COMSOL MULTIPHYSICS WITH MATLAB 多变量微分计算的入门读本。

当我们写该附录时曾经争论过是否将这部分内容直接加入到第一章（数值分析基础）中，因为导数的数值近似是偏微分方程求解的基础，而偏微分方程是COMSOL MULTIPHYSICS 的基本运算单元。

确实，在学习波谱法求解偏微分方程时，基本理论就是“导数理论”——如何使用波变换方法来近似导数。

所以通过对比发现，有限元方法的基础就是数值微分。

所以争论就不存在了，第一章主要是关于COMSOL MULTIPHYSICS 直接计算的基本问题的。

但是不管多有用，近似导数仍然只是建模的一个中间步骤，不是目标本身。

我们这里只考虑用于矢量计算的MATLAB 基础，本附录的重点在于特征值分析和逻辑表达式。

这些在整本书中都有体现。

应当注意到我们这里介绍的每个功能都可以在COMSOL Script 中实现。

本书中唯一不能在COMSOL Script 中实现的Matlab 命令就是fminsearch 。

1．矢量回顾1.1 矢量表达FEMLAB 可以处理标量、矢量和矩阵数据，这里简单介绍一下矢量的表达（作为MATLAB 矩阵数据类型的一个特例）。

标量可以作为一个单独的数，但是矢量是具有大小和方向的。

在如图1所示的右手坐标系系统中，向量a 用以下形式表达：123123(,,)a a a a a a =++=a i j k a (1)这里i ，j 和k 是坐标方向的单位矢量，1a ，2a ，3a 是向量a 在各轴方向上的分量。

它们是a 对各单位矢量i ，j 和k 的投影。

对于坐标系中的P 点（x ，y ，z ），矢量P 对于初始坐标系统O 的位置为：(,,)x y z x y z =++=r i j k(2)MATLAB 用分量的形式描述列矢量或行矢量：>> a = [1; 2; 3]; % column vector>> a = [1 2 3]; % row vector在行向量中，空白（任意连续空格）作为分界符。

列向量用分号或者回车符分界：>> a = [12 3];1.2 内积，矩阵乘法，单位矢量和矢积典型的内积（或点积）定义为：31122331cos i i i a b a b a b a b a b θ=⋅==++=∑a b(3)这里θ是矢量a 和b 的夹角。

为了在MATALB 中达到相同的目的，我们使用*运算符：>> a = [1; 2; 3]; >> b = [-3 2 -1]; >> b*a ans = -2这是一个行向量（1×3矩阵）乘以列向量（3×1矩阵）的特殊情况。

因为前者的列数和后者的行数相等，这两个矩阵是相容的，可以根据矩阵乘法通用法则计算。

()1nij jk ik j AB A B ==∑(4)如果A 是m ×n 矩阵，B 是n ×l 矩阵，则AB 是m ×l 矩阵。

如果共用的维数不同，那么矩阵不相容，不能定义乘法运算。

MATLAB 也可以将标量乘法作为特殊矩阵乘法来计算，但是必须考虑矩阵的相容性。

例如>> a*b ans =-3 2 -1 -6 4 -2 -9 6 -3出现了什么情况？很简单，a 是3×1矩阵，乘以1×3矩阵b ，得到的ab 是3×3矩阵。

()111nij jk i k ik j ab a b a b ===∑(5)对于向量，矩阵(ab )ik 称为a 和b 的并积或并矢。

这是矩阵外积的特殊情况，这里标量乘积也算内积。

在MATLAB 中通过转置运算符“’”可以实现两个行向量或两个列向量的内积，它是一个一元运算符，容易被误解为英语中的缩写符号。

ans =-2 but>> a*b’ans =-3 2 -1-6 4 -2-9 6 -3仍然产生并矢。

必须自己考虑矩阵的相容性。

如果a和b是行向量，那么'b*a 或a*'b将产生内积还是外积？出于这个原因，MATLAB提供了一个特殊的dot 函数来取消这种相容性的差别。

>> help dotDOT Vector dot product.C = DOT(A,B) returns the scalar product of the vectors A and B.A andB must be vectors of the same length. When A and B are bothcolumn vectors, DOT(A,B) is the 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 2 -1])ans =-2使用dot函数使得行/列向量的组合不受影响。

向量值向量的值或模可以通过下式计算：1/221njja=⎛⎫== ⎪⎝⎭∑a(6) MATLAB通过下式计算向量的模。

ans =3.7417or with the built-in command norm>> norm(a,2)ans =3.7417这里sqrt( )是预置的平方根函数。

单位向量单位向量是模为1的向量。

单位向量可以通过归一化处理得到，即：ˆ=aaa(7)例如，>> ahat=a/norm(a,2) ahat =0.2673 0.5345 0.8018以上除法是标量除法，是矢量的每个元素除以标量。

COMSOL Multiphysics 通常根据预置的几何单位向量来描述边界。

nx ，ny ，nz 通常作为边界上指向外部的单位矢量。

较少使用tx 和ty 作为边界曲线的切线矢量。

通常二维表面有两个正交切线方向，所以有无穷多个切向矢量。

叉积叉积定义如下：31ˆˆsin ijk j k i i na b e θε=⨯==∑a b a b (8)这里ijk ε是置换张量，当指数ijk 为123的正置换时取＋1，当是123的负置换时取－1，其它情况均为零。

ˆn是包含a 和b 的平面内的单位标准向量。

ˆi e 是第i 个坐标方向的单位向量。

MATLAB 提供了一个特殊的函数来计算叉积。

>> help crossCROSS Vector cross product.C = CROSS(A,B) returns the cross product of the vectors A and B. That is, C = A x B. A and B must be 3 element vectors.C = CROSS(A,B) returns the cross product of A and B along the first dimension of length 3.C = CROSS(A,B,DIM), where A and B are N-D arrays, returns the cross product of vectors in the dimension DIM of A and B. A and B must have the same size, and both SIZE(A,DIM) and SIZE(B,DIM) must be 3. See also DOT. For example, >> cross(a,b) ans =-8 -8 8 >> cross(b,a) ans =88-8我们看到叉积中的因子阶数决定了叉积的正负号，等效于改变了单位标准向量ˆn的方向。

1.3 数组：简单数组，元包数组和结构体数组处理本质上就是从FEMLAB中提取数据。

为方便起见，FEMLAB对于多物理场采用fem结构体组织模型，对于扩展多物理场采用xfem结构体组织模型。

从结构体和元包数组中提取有用信息是访问FEMLAB模型（和结果）的一个非常有用的方法。

简单数组数组有维数（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实际上冒号相当于低维数的子数组。

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前两种情况下子数组是矩阵，第三种情况下最终两个维数聚集成一个简单的行向量。