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顯