1、附錄 COMSOL Multiphysics的MATLAB矢量計(jì)算基礎(chǔ)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矢量計(jì)算支撐了偏微分方程和它們的數(shù)值近似求解。為了很

2、好的使用有限元方法，建模人員應(yīng)該掌握矢量計(jì)算基礎(chǔ)知識(shí)。本科畢業(yè)的工程師可能學(xué)過(guò)矢量計(jì)算的數(shù)學(xué)課程，但是由于沒(méi)有碰到過(guò)矢量計(jì)算的實(shí)際應(yīng)用，這時(shí)在工程建模中使用矢量計(jì)算就受到限制。本附錄介紹了所有COMSOL MULTIPHYSICS WITH MATLAB中用到的矢量計(jì)算基礎(chǔ)知識(shí)。所以也可以將該附錄當(dāng)作是COMSOL MULTIPHYSICS WITH MATLAB多變量微分計(jì)算的入門(mén)讀本。當(dāng)我們寫(xiě)該附錄時(shí)曾經(jīng)爭(zhēng)論過(guò)是否將這部分內(nèi)容直接加入到第一章（數(shù)值分析基礎(chǔ)）中，因?yàn)閷?dǎo)數(shù)的數(shù)值近似是偏微分方程求解的基礎(chǔ)，而偏微分方程是COMSOL MULTIPHYSICS的基本運(yùn)算單元。確實(shí)，在學(xué)習(xí)波譜法求

3、解偏微分方程時(shí)，基本理論就是“導(dǎo)數(shù)理論”如何使用波變換方法來(lái)近似導(dǎo)數(shù)。所以通過(guò)對(duì)比發(fā)現(xiàn)，有限元方法的基礎(chǔ)就是數(shù)值微分。所以爭(zhēng)論就不存在了，第一章主要是關(guān)于COMSOL MULTIPHYSICS直接計(jì)算的基本問(wèn)題的。但是不管多有用，近似導(dǎo)數(shù)仍然只是建模的一個(gè)中間步驟，不是目標(biāo)本身。我們這里只考慮用于矢量計(jì)算的MATLAB基礎(chǔ)，本附錄的重點(diǎn)在于特征值分析和邏輯表達(dá)式。這些在整本書(shū)中都有體現(xiàn)。應(yīng)當(dāng)注意到我們這里介紹的每個(gè)功能都可以在COMSOL Script中實(shí)現(xiàn)。本書(shū)中唯一不能在COMSOL Script中實(shí)現(xiàn)的Matlab命令就是fminsearch。1矢量回顧1.1 矢量表達(dá)FEMLAB可以處

4、理標(biāo)量、矢量和矩陣數(shù)據(jù)，這里簡(jiǎn)單介紹一下矢量的表達(dá)（作為MATLAB矩陣數(shù)據(jù)類(lèi)型的一個(gè)特例）。標(biāo)量可以作為一個(gè)單獨(dú)的數(shù)，但是矢量是具有大小和方向的。在如圖1所示的右手坐標(biāo)系系統(tǒng)中，向量a用以下形式表達(dá)： (1)這里，和是坐標(biāo)方向的單位矢量，是向量在各軸方向上的分量。它們是對(duì)各單位矢量，和的投影。對(duì)于坐標(biāo)系中的P點(diǎn)（x，y，z），矢量P對(duì)于初始坐標(biāo)系統(tǒng)O的位置為： (2)MATLAB用分量的形式描述列矢量或行矢量：>> a = 1; 2; 3; % column vector>> a = 1 2 3; % row vector在行向量中，空白（任意連續(xù)空格）作為分界符。列

5、向量用分號(hào)或者回車(chē)符分界：>> a = 123;1.2 內(nèi)積，矩陣乘法，單位矢量和矢積典型的內(nèi)積（或點(diǎn)積）定義為： (3)這里是矢量和的夾角。為了在MATALB中達(dá)到相同的目的，我們使用*運(yùn)算符：>> a = 1; 2; 3;>> b = -3 2 -1;>> b*aans =-2這是一個(gè)行向量（1×3矩陣）乘以列向量（3×1矩陣）的特殊情況。因?yàn)榍罢叩牧袛?shù)和后者的行數(shù)相等，這兩個(gè)矩陣是相容的，可以根據(jù)矩陣乘法通用法則計(jì)算。 (4)如果A是m×n矩陣，B是n×l矩陣，則AB是m×l矩陣。如果共用的

6、維數(shù)不同，那么矩陣不相容，不能定義乘法運(yùn)算。MATLAB也可以將標(biāo)量乘法作為特殊矩陣乘法來(lái)計(jì)算，但是必須考慮矩陣的相容性。例如>> a*bans =-3 2 -1-6 4 -2-9 6 -3出現(xiàn)了什么情況？很簡(jiǎn)單，a是3×1矩陣，乘以1×3矩陣b，得到的ab是3×3矩陣。 (5)對(duì)于向量，矩陣(ab)ik稱(chēng)為a和b的并積或并矢。這是矩陣外積的特殊情況，這里標(biāo)量乘積也算內(nèi)積。在MATLAB中通過(guò)轉(zhuǎn)置運(yùn)算符“”可以實(shí)現(xiàn)兩個(gè)行向量或兩個(gè)列向量的內(nèi)積，它是一個(gè)一元運(yùn)算符，容易被誤解為英語(yǔ)中的縮寫(xiě)符號(hào)。ans =-2but>> a*bans =-3

7、2 -1-6 4 -2-9 6 -3仍然產(chǎn)生并矢。必須自己考慮矩陣的相容性。如果a和b是行向量，那么*a或a*將產(chǎn)生內(nèi)積還是外積？出于這個(gè)原因，MATLAB提供了一個(gè)特殊的dot函數(shù)來(lái)取消這種相容性的差別。>> 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函數(shù)使得行/列向量的組合不受影響。向量值向量的值或模可以通過(guò)下式計(jì)算： (6)MATLAB通過(guò)下式計(jì)算向量的模。ans =3.7417or with the built-in command norm>> norm(a,2)ans =3.7417這里sqrt( )是預(yù)置的平方根函數(shù)。單位向量單位向量是模為1的向量。單位向量可以通過(guò)歸一化處理得到，即： (7)例如，>> ahat=a/norm(a,2)ahat =0.26730.53450.8018以上除法是標(biāo)量除法，是矢量的每個(gè)元素除以標(biāo)量。COMSOL Multiphysics通常根

10、據(jù)預(yù)置的幾何單位向量來(lái)描述邊界。nx，ny，nz通常作為邊界上指向外部的單位矢量。較少使用tx和ty作為邊界曲線的切線矢量。通常二維表面有兩個(gè)正交切線方向，所以有無(wú)窮多個(gè)切向矢量。叉積叉積定義如下： (8)這里是置換張量，當(dāng)指數(shù)ijk為123的正置換時(shí)取1，當(dāng)是123的負(fù)置換時(shí)取1，其它情況均為零。是包含a和b的平面內(nèi)的單位標(biāo)準(zhǔn)向量。是第i個(gè)坐標(biāo)方向的單位向量。MATLAB提供了一個(gè)特殊的函數(shù)來(lái)計(jì)算叉積。>> 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我們看到叉積中的因子階數(shù)決定了叉積的正負(fù)號(hào)，等效于改變了單位標(biāo)準(zhǔn)向量的方向。1.3 數(shù)組：簡(jiǎn)單數(shù)組，元包數(shù)組和結(jié)構(gòu)體數(shù)組處理本質(zhì)上就是從FEMLAB中提取數(shù)據(jù)。為方便起見(jiàn)，F(xiàn)EMLAB對(duì)于多物理場(chǎng)采用fem結(jié)構(gòu)體組織模型，對(duì)于擴(kuò)展多物理場(chǎng)采用xfem結(jié)構(gòu)體組織模型。從結(jié)構(gòu)體和元包數(shù)組中提取有用信息是訪問(wèn)FEMLAB模型（和結(jié)果

13、）的一個(gè)非常有用的方法。簡(jiǎn)單數(shù)組數(shù)組有維數(shù)（m×n×l）。矩陣是二維數(shù)組。每個(gè)維數(shù)都有長(zhǎng)度。所以size( )和length( )是兩個(gè)非常有用的命令。>> a = 1 2 3 4; 5 6 7 8;>> size(a)ans =2 4數(shù)組的尺寸本身就是一個(gè)行向量，它的長(zhǎng)度等于數(shù)組的維數(shù)。>> length(a(1,:)ans =4第二個(gè)自變量的冒號(hào)（:）表示包括第二維的整個(gè)范圍，在本例中表示元素1到4。>> length(a(:,3)ans =2>> a(1,2:4)ans =2 3 4實(shí)際上冒號(hào)相當(dāng)于低維數(shù)的子

14、數(shù)組。a(1,:)是第一行；a(:,3)是a的第三列。a(1,2:4)給出了第1行2到4列元素組成的子數(shù)組。對(duì)于更高維數(shù)情況，提取出的子數(shù)組更加復(fù)雜，例如：>>b=ones(2,2,2)b(:,:,1) =1 11 1b(:,:,2) =1 11 1>> b(1,:)ans =1 1 1 1前兩種情況下子數(shù)組是矩陣，第三種情況下最終兩個(gè)維數(shù)聚集成一個(gè)簡(jiǎn)單的行向量。FORTRAN編程人員可能感覺(jué)單個(gè)元素的尋址比子數(shù)組更舒服，可能會(huì)使用循環(huán)結(jié)構(gòu)來(lái)達(dá)到相同的目的。>> a(1,3)ans =3構(gòu)造數(shù)組數(shù)組可以通過(guò)冒號(hào)元算符自動(dòng)生成，即>> 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這樣產(chǎn)生了從0到間隔相等的11個(gè)值。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是一個(gè)通用的矩陣自動(dòng)生成命令，它起著舊編程語(yǔ)言

16、中循環(huán)結(jié)構(gòu)的作用。有的時(shí)候也用得到logspace函數(shù)。>> 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, :.四種其它常見(jiàn)數(shù)組生成函數(shù)為zeros，ones，rand和eye。zeros生成全0數(shù)組；ones生成全1數(shù)組；rand生成均勻分布的隨機(jī)數(shù)組；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.到目前為止你可能已經(jīng)注意到，如果知道命令的名字，M

22、atlab就可以給出相關(guān)語(yǔ)法。但是你怎么才能知道命令的名字呢？在COMSOL Script中，只要輸入help就會(huì)產(chǎn)生一列命令類(lèi)型。隨后輸入“help type”命令就可以顯示出相關(guān)命令。例如，“help numerics”給出了一列數(shù)值命令。然后“help daspk”給出調(diào)用daspk命令的語(yǔ)法，這是單獨(dú)為FEMLAB準(zhǔn)備的剛性DAE求解器。標(biāo)量數(shù)組運(yùn)算對(duì)數(shù)組的標(biāo)量四則運(yùn)算等于對(duì)數(shù)組中元素的四則運(yùn)算，例如：>> 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數(shù)組數(shù)組的基元運(yùn)算數(shù)組對(duì)數(shù)組的四則運(yùn)算是一個(gè)比較偏的內(nèi)容。如果數(shù)組大小一致，那么可以對(duì)元素采用點(diǎn)運(yùn)算符：>> 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元包數(shù)組和結(jié)構(gòu)體關(guān)于元包數(shù)組和結(jié)構(gòu)體可以寫(xiě)整整一章的內(nèi)容。對(duì)于這兩種數(shù)據(jù)結(jié)構(gòu)，應(yīng)該注意到它們是不同種類(lèi)混合數(shù)據(jù)類(lèi)型的容器包括浮點(diǎn)數(shù)，復(fù)數(shù)，矩陣，字符串和其它數(shù)組和結(jié)構(gòu)體。元包數(shù)組將數(shù)組中的每個(gè)元包作為指針。元包數(shù)組通過(guò)包含數(shù)組基元的括弧來(lái)直接定義。元包命令將會(huì)產(chǎn)生一個(gè)空的構(gòu)架，你可以在其中分別填入元素或

25、數(shù)組。>> 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通過(guò)圓括號(hào)內(nèi)的數(shù)組地址可以進(jìn)行調(diào)用，返回元包的元素。>> ca(3)ans =boy通過(guò)大括號(hào)內(nèi)的地址數(shù)可以返回元包中基元包含的內(nèi)容。>> ca3ans =boy可能通過(guò)矩陣元包更容易理解這一點(diǎn)。>> ca(7)ans =2x2 double>> ca7ans =

26、0 1-2 2結(jié)構(gòu)體更多的借鑒了域的概念，也常稱(chēng)作枚舉法，在C語(yǔ)言中比較常見(jiàn)。使用結(jié)構(gòu)體和元包數(shù)組的最大區(qū)別在于如果你改變結(jié)構(gòu)體，增加或減少域并不會(huì)改變域的次序。但是減少或增加元包數(shù)組元素會(huì)改變?cè)幪?hào)，或在數(shù)組中產(chǎn)生“空穴”。COMSOL Multiphysics用戶(hù)最常見(jiàn)到的結(jié)構(gòu)體就是fem結(jié)構(gòu)體，這是COMSOL Multiphysics對(duì)多物理場(chǎng)模型組織數(shù)據(jù)和計(jì)算結(jié)果的結(jié)構(gòu)體。從COMSOL Multiphysics中導(dǎo)出fem結(jié)構(gòu)體到MATLAB工作空間中產(chǎn)生以下電動(dòng)流動(dòng)模型。>> 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結(jié)構(gòu)體的域列表給出了對(duì)各域中內(nèi)容的描述。每個(gè)域都可以通過(guò)“點(diǎn)”號(hào)來(lái)尋址。>> fem.exprans =zeta 1x22 char sig Y+sigr*(1-Y)fem.expr是有四個(gè)元包的元包數(shù)組；該元包數(shù)組非常小，能夠顯示所有的內(nèi)容。由于是個(gè)元包數(shù)組，可以通過(guò)大括號(hào)來(lái)調(diào)用元包中的內(nèi)容。>> fem.expr2ans =-zeta1*(Y+zetar*(1-Y)正如我們看到的，fem結(jié)構(gòu)體中含有元包數(shù)組，其它結(jié)構(gòu)體，字符和數(shù)字。我們也可以用fem結(jié)構(gòu)體來(lái)

29、構(gòu)成元包數(shù)組，這就是COMSOL Multiphysics處理擴(kuò)展多物理場(chǎng)時(shí)構(gòu)建的xfem結(jié)構(gòu)體，其中每個(gè)fem結(jié)構(gòu)體對(duì)應(yīng)一個(gè)邏輯幾何結(jié)構(gòu)。我們通常需要對(duì)fem和xfem結(jié)構(gòu)體調(diào)用postinterp命令，來(lái)插值取得有限元離散中間點(diǎn)處的變量計(jì)算值：is,pe=postinterp(xfem,xx);u=postinterp(xfem,u1,is);通過(guò)對(duì)xfem結(jié)構(gòu)體調(diào)用postinterp函數(shù)可以得到對(duì)模型和解的完整描述，對(duì)于不同的結(jié)構(gòu)體會(huì)執(zhí)行不同的命令分支。作為COMSOL Mulitphysics用戶(hù)，有必要知道COMSOL Multiphysics模型和解中用到的MATLAB數(shù)據(jù)結(jié)構(gòu)，

30、這樣當(dāng)我們希望進(jìn)行特殊后處理或建模（而不是在COMSOL Multiphysics圖形用戶(hù)界面中實(shí)現(xiàn)）時(shí)，可以從中提取相應(yīng)的數(shù)據(jù)。2標(biāo)量和矢量場(chǎng)：MATLAB函數(shù)的表達(dá)物質(zhì)屬性通常取決于不同的位置或時(shí)間。在人們可見(jiàn)的長(zhǎng)度尺度上，大多數(shù)物理量都看作連續(xù)函數(shù)在每個(gè)數(shù)值點(diǎn)都有相應(yīng)值。這些量就被稱(chēng)作場(chǎng)。像溫度和壓力等描述單個(gè)值的量就是標(biāo)量場(chǎng)。標(biāo)量場(chǎng)是一個(gè)單一的數(shù)值，即 (9)三維矢量場(chǎng)需要三個(gè)分量： (10)F的每個(gè)分量都是位置的標(biāo)量函數(shù)。例：圖1 對(duì)于不同C值的等高線圖1給出了常數(shù)C在不同值時(shí)的等高線圖。MATLAB中沒(méi)有數(shù)據(jù)類(lèi)型來(lái)描述場(chǎng)量，這些量可以以函數(shù)形式來(lái)表達(dá)。在MATLAB中有三種主要的函

31、數(shù)表達(dá)方式：(1) inline函數(shù)，只在工作空間中定義>> 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菜單下有個(gè)新功能專(zhuān)門(mén)用來(lái)以公式形式定義inline函數(shù)。(2) m文件函數(shù)，它以文件形式存儲(chǔ)在硬盤(pán)中，可以被工作空間或m文件代碼調(diào)用。例如m文件函數(shù)tem

32、perat.m包含以下代碼，保存在MATLAB當(dāng)前文件夾中。function t=temperat(r)%TEMPERAT evaluates T = 1 + ln r% T = temperat(r)%t=1+log(r);>>a=temperat(10)a =3.3026 (3) 內(nèi)插函數(shù)。有時(shí)候確定了函數(shù)在某些特定點(diǎn)的值，而在附近一些點(diǎn)的值就需要假設(shè)函數(shù)局部光滑來(lái)計(jì)算得到。內(nèi)插需要一系列的MATLAB函數(shù)，但是最終產(chǎn)生一種函數(shù)形式。對(duì)于一維、二維和三維數(shù)據(jù)，MATLAB有預(yù)置的插值函數(shù)，分別稱(chēng)作interp1，interp2和interp3。COMSOL Multiphysi

33、cs的Options菜單下有個(gè)新功能以手動(dòng)輸入數(shù)據(jù)或給定格式的數(shù)據(jù)文件來(lái)定義內(nèi)插函數(shù)。第二章122頁(yè)舉了個(gè)將水密度分布作為溫度函數(shù)輸入的例子。通常情況下，COMSOL Multiphysics中系數(shù)和邊界條件等場(chǎng)變量都是通過(guò)預(yù)定義獨(dú)立變量、因變量、派生變量串聯(lián)方式輸入的。例如，在通用模式下，可以輸入單一因變量u和獨(dú)立變量x的表達(dá)式：但是MATLAB m文件函數(shù)可以很容易的調(diào)用。這里要指出的很重要一點(diǎn)是，COMSOL Multiphysics通常以標(biāo)量分量形式接受輸入數(shù)據(jù)。如果需要一個(gè)向量或矩陣，就需要通過(guò)對(duì)標(biāo)量分量形式輸入，其中每個(gè)分量都可以是復(fù)雜函數(shù)。另外很重要的一點(diǎn)是，COMSOL Mul

34、tiphysics允許函數(shù)有矢量輸入，即每個(gè)函數(shù)的自變量都可以是矢量。所以你在m文件函數(shù)中使用顯式元素運(yùn)算符，例如“.*”時(shí)，必須非常小心。COMSOL Multiphysics用fem結(jié)構(gòu)體來(lái)記錄計(jì)算結(jié)果，對(duì)于fem.mesh中定義的網(wǎng)格，自由度保存在fem.sol中。對(duì)于每個(gè)因變量和派生變量，COMSOL Multiphysics提供了一個(gè)特殊的內(nèi)插函數(shù)來(lái)從fem.sol中取得插值。本書(shū)例題中大量使用postinterp內(nèi)插函數(shù)。它甚至可以在m文件函數(shù)中自動(dòng)調(diào)用mat文件中的fem結(jié)構(gòu)體。3多變量計(jì)算中的微分3.1 標(biāo)量場(chǎng)梯度如果，那么矢量 (11)被稱(chēng)作標(biāo)量場(chǎng)的梯度，通常表示為grad。

35、梯度運(yùn)算符是三維直角坐標(biāo)系中的矢量運(yùn)算符。 (12)COMSOL Multiphysics舉例假設(shè)，則。但是MATLAB不直接處理這類(lèi)符號(hào)計(jì)算（它的符號(hào)工具箱可以）。COMSOL Multiphysics粗略計(jì)算微分的數(shù)值近似解。所以標(biāo)量場(chǎng)的梯度可以通過(guò)COMSOL Multiphysics的“基元”運(yùn)算符得到。我們是如何輕易得到這一信息的呢？表1中給出了具體做法。應(yīng)當(dāng)注意到，由于實(shí)際上沒(méi)有偏微分方程被求解，Neumann邊界條件等價(jià)于中性或無(wú)條件邊界。否則的話(huà)，只有邊界數(shù)據(jù)滿(mǎn)足0=phi-x2-y2時(shí)才是可能的解。表1 梯度模型模型導(dǎo)航欄二維，PDE模式，通用模式，獨(dú)立變量：x，y 因變量：

36、phi選項(xiàng)設(shè)置Axes/Grid為-1，1×-1，1繪圖矩形域-1，1×-1，1邊界模式/邊界設(shè)定設(shè)定四個(gè)邊界均為Neumann邊界條件子域模式/子域設(shè)定在域1中，設(shè)定0 0；da0；Fphi-x2-y2網(wǎng)格模式重繪默認(rèn)網(wǎng)格（418節(jié)點(diǎn)，774基元）求解使用默認(rèn)設(shè)定（非線性求解器）后處理選擇arrow模式，如圖A4現(xiàn)在將fem結(jié)構(gòu)導(dǎo)入MATLAB中。我們通過(guò)MATLAB雙線性回歸插值得到近似數(shù)值解。>>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命令都是確保數(shù)據(jù)滿(mǎn)足interp2的正確格式。表2給出了其它點(diǎn)的梯度計(jì)算。表2 使用梯度模型對(duì)grad的數(shù)值計(jì)算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通過(guò)任意解法，F(xiàn)EM得到的一階導(dǎo)數(shù)都十分準(zhǔn)確。根據(jù)收斂準(zhǔn)則，全局誤差為O(10-16)。圖2給出了梯度的數(shù)值近似解。圖2 grad的矢量箭頭圖3.2 方向?qū)?shù)的方向?qū)?shù)就是沿給定方向的變化速率。如果是該方向的單位向量，則方向?qū)?shù)可以表示為： (13)坐標(biāo)方向非常容易求解，即： (14)在第七章的ECT模型中我們使用方向?qū)?shù)計(jì)算電勢(shì)的法向?qū)?shù)。顯然，方向?qū)?shù)與通量的概念密切相關(guān)。對(duì)于“線性”特性系統(tǒng)（Fick定律，F(xiàn)ourier定律等），穿過(guò)材料表面的總通量正比于表面法向?qū)?shù)的積分。局部通量正比于某一點(diǎn)的法向?qū)?shù)。對(duì)于這一點(diǎn)，大多數(shù)矢量計(jì)算文章都已經(jīng)證

39、明，的變化速率最大的方向就是grad的方向，且|grad|就是在該方向的變化速率。對(duì)于點(diǎn)(x,y)=(0.25,-0.75)，通過(guò)沿角度的步進(jìn)法繪制，就能清楚的看到這一點(diǎn)。相關(guān)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編程習(xí)慣，以上代碼使用了循環(huán)結(jié)構(gòu)。如果是對(duì)MATLAB編程比較熟悉，一般不使用循環(huán)法：&

40、gt;>dirder = u*cos(theta)+v*sin(theta);plot(theta,dirder)cos和sin函數(shù)會(huì)對(duì)向量thera中的每個(gè)元素進(jìn)行計(jì)算，生成相同長(zhǎng)度的向量。圖3 方向角度隨弧度變化圖。注意存在最小的方向?qū)?shù)下降最陡的方向，對(duì)應(yīng)于梯度方向。3.3 水平集合/水準(zhǔn)面注意到方向?qū)?shù)穿過(guò)x軸，即某一個(gè)方向的方向?qū)?shù)為零在該方向沒(méi)有任何變化速率。可以看出的方向垂直于梯度方向。所以在該方向局部為常數(shù)。繪制出每個(gè)為常數(shù)的曲線（二維）或曲面（三維），就會(huì)形成一簇曲線（曲面），稱(chēng)作的水平集合（見(jiàn)第八章）。在二維空間中，水平幾何也被稱(chēng)作等高線。方向?qū)?shù)的概念非常類(lèi)似于測(cè)量

41、圖，就是陸地的海拔。等高線具有相同的海拔；方向?qū)?shù)是在方向的爬山速率，梯度就是最陡的爬山方向，爬山速率為|grad|。在流體動(dòng)力學(xué)中，通常用流函數(shù)的等高線來(lái)描述，等高線就是正切于速度場(chǎng)的流線（穩(wěn)定流動(dòng)中的顆粒路徑）。在第三章浮力對(duì)流的例子中介紹了如何計(jì)算流函數(shù)（見(jiàn)第三章方程(3)）。3.4 矢量場(chǎng)導(dǎo)數(shù)矢量微分運(yùn)算符可以通過(guò)兩種方式作用于矢量場(chǎng)：(1)內(nèi)積形式，稱(chēng)作散度；(2)叉積形式，稱(chēng)作旋度。散度為： (15)旋度為： (16)是前面引入的置換張量。很容易會(huì)發(fā)現(xiàn)散度是標(biāo)量，而旋度是矢量。運(yùn)算符經(jīng)常在傳熱或傳質(zhì)方程的對(duì)流項(xiàng)中見(jiàn)到。它顯然不是散度，因?yàn)?(17)和方程(15)的標(biāo)量相比較，它仍然

42、是個(gè)運(yùn)算符。正如我們前面看到的，導(dǎo)數(shù)的數(shù)值近似是COMSOL Multiphysics的基本功能，所以我們可以用來(lái)近似計(jì)算散度和旋度。COMSOL Multiphysics舉例假設(shè)，表3中給出具體計(jì)算步驟。表3 求解散度和旋度步驟模型導(dǎo)航欄三維，PDE模式，通用模式，獨(dú)立變量：x，y，z因變量：u1，u2，u3選項(xiàng)設(shè)置Axes/Grid為0，1×0，1×0，1繪圖繪制六面體BLK10，1×0，1×0，1邊界模式/邊界設(shè)定設(shè)定所有邊界為Neumann邊界條件子域模式/子域設(shè)定設(shè)定0 0 0; da1=0 0 0; F1=u1-x2設(shè)定0 0 0; da2=

43、0 0 0; F2=u2-3*x*y設(shè)定0 0 0; da3=0 0 0; F3=u3-x3網(wǎng)格模式設(shè)定網(wǎng)格縮放因子為3，重繪網(wǎng)格（201節(jié)點(diǎn)，719自由度）求解使用默認(rèn)設(shè)定（非線性求解器）后處理(1)繪制散度u1x+u2y+u3z的彩色圖(2)繪制旋度(u3y-u2z,u1z-u3x,u2x-u1y)的箭頭圖應(yīng)當(dāng)再次注意到，實(shí)際上沒(méi)有求解任何偏微分方程，所以Neumann邊界條件等價(jià)于中性和無(wú)條件邊界。否則的話(huà)，只有邊界數(shù)據(jù)滿(mǎn)足時(shí)才是可能的解。通過(guò)符號(hào)計(jì)算很容易算得：數(shù)值解近似程度如何？嘗試以下散度：>> xxx=0.42; 0.57; 0.33;postinterp(fem,u1x+u2y+u3z,xxx)ans = 2.1137>> 5*0.42ans = 2.1000顯