COMSOLMULTIPHYSICS和數(shù)值分析基礎(chǔ)

第一章COMSOLMＵＬTＩPHＹＳICS及數(shù)值分析基礎(chǔ)W。B．Ｊ．ZIMMERＭAN,B.Ｎ.HEWAKAＮDＡＭBＹDepａrtmentｏfCｈemicalandＰrｏｃessＥnｇinｅｅriｎg，UｎiｖersityoｆSｈｅfｆield，NｅwｃasｔｌeＳｔrｅet，ＳｈefｆｉeldＳ13JDUniｔeｄKiｎgｄomＥ—ｍail：ｗ。ｚｉｍｍeｒman@shｅf.aｃ．uk本章主要介紹ＣOＭＳOＬＭuｌtipｈｙsｉｃｓ在零維和一維模型數(shù)值分析方面的幾個關(guān)鍵內(nèi)容。這些內(nèi)容包括求根、步進(jìn)式數(shù)值積分、常微分方程數(shù)值積分和線性系統(tǒng)分析。這幾乎是全部的化工過程數(shù)學(xué)分析方法。下面通過CＯＭSOLMｕltｉpｈｙsｉcｓ中的一些常見化工過程應(yīng)用實例來介紹這些方法，包括：閃蒸、管式反應(yīng)器設(shè)計、集中反應(yīng)系統(tǒng)和固體中熱傳導(dǎo)。1．簡介本章內(nèi)容很多，可以分為幾個不同的目標(biāo)。首先介紹了CＯMSOＬMｕltｉpｈysiｃs的主要工作特性;其次介紹了如何使用這些特性來分析一些簡潔的，位于零維空間、一維空間或“空間－時間"系統(tǒng)中的化工問題.本章還盼望通過展現(xiàn)COMＳOLMultiphｙｓiｃs和ＭＡＴＬAＢ工具在化工過程分析中的強大功能,激發(fā)讀者對使用COMSOLＭultｉphｙsicｓ進(jìn)行建模與仿真的愛好。由于CＯＭSOLMultｉphyｓｉｃｓ不是一個通用的問題求解工具，所以一些目標(biāo)需要迂回實現(xiàn)。作者在使用FＯRＴRAN、Maｔｈｅｍａt(yī)ｉca和MATＬＡB解決化工問題方面有著豐富的教學(xué)閱歷，并用這些工具實現(xiàn)過這里全部的例子。而且,擴展化工問題的數(shù)值分析也已經(jīng)在PＯLYMＡTH［1]中實現(xiàn)，這似乎只在化工委員會的CAＣHE項目中使用過。本書前一版已經(jīng)介紹過在零維空間中求解非線性代數(shù)方程和與時間有關(guān)的常微分方程的內(nèi)容.從概念上講,零維域就是一個簡潔的有限元。通過討論某一特定有限元中的變化對理解有限元方法格外有用。但是,CＯMＳＯLMｕｌｔｉphyｓicｓ通過獨立對話框設(shè)置，使得零維幾何方程和與時間相關(guān)的常微分方程求解變得格外簡潔。所以本章將同時采納這兩種方法求解這些例子。2．方法１:求根典型的數(shù)值分析課程會講解多種求根方法，但是從實際閱歷來看,只有兩種算法格外有用——二分法和牛頓法。我們這里沒有列出全部方法,而是重點考慮為什么求根是最有效的數(shù)值分析工具。在線性系統(tǒng)中求根格外簡潔，但是對于非線性系統(tǒng)這就是一個挑戰(zhàn),而全部感愛好的動力學(xué)問題幾乎都是非線性系統(tǒng).對非線性系統(tǒng)的求根起源于對反函數(shù)的描述。為什么呢?由于對于大多數(shù)非線性函數(shù),“正向"y=ｆ(ｕ）很好表示，但是它的反函數(shù)u＝f－-1(ｙ)可能不能顯式表示、多值（無意義）或根本不存在。如果反函數(shù)存在的話,求解反函數(shù)其實就是求根的過程——求解滿意Ｆ（u)＝0的ｕ等價于求解F(ｕ）=f(u)－y＝0。由于大多數(shù)數(shù)值分析的目標(biāo)是在系統(tǒng)約束下計算求解，所以這也等價于對全部的約束取反.COＭSOLMｕｌtｉｐhysｉcｓ擁有求解非線性問題的核心函數(shù)—-ｆｅｍnｌiｎ，本節(jié)主要介紹用它求解零維非線性問題。feｍnlｉn函數(shù)使用牛頓方法求解,由于只有一個變量u,牛頓法通過對一階倒數(shù)迭代來求根。該方法首先估量函數(shù)的斜率范圍，然后再逼近根。該斜率可以通過理論分析(牛頓—拉夫遜方法）和數(shù)值（正割法)方法求得。如果能用任何一種方法求得斜率，就可以用泰勒定律來逼近根。其基本思想就是使用目前推測值u０的泰勒展開式：??????(1）該公式可以化簡，忽視(ｕ—ｕ0）的高階項,計算根如下:?? ?????(2)這個方法可以快速地擴展到多維求解空間,例如將ｕ看作未知矢量，“被除”看作“乘以f的雅克比矩陣的逆”。下一節(jié)介紹ＣＯMＳOLMultiphyｓiｃs中的求根過程。2。1求根:COMＳOＬMulｔiphysｉcｓ非線性求解器的應(yīng)用實例如上節(jié)所述,求根本身是一個“零維”活動，至少對于“空間－時間”系統(tǒng)多維未知矢量ｕ來說是這樣的。ＣOMSＯＬ多物理場沒有零維模式,所以我們臨時采納一維模式。這在方面增加了我們不需要的冗余功能。但是由于問題規(guī)模較小,ＣOMＳＯLＭｕlｔiphyｓｉcｓ編碼效率高，且現(xiàn)代微處理器的運算速度快,這點就不成為問題了。啟動ＭＡTLＡB并在命令窗口鍵入COMＳOLMulｔｉｐｈysicｓ。屏幕閃爍幾秒后，會消滅一個模型導(dǎo)航窗口。依據(jù)表1所示步驟,建立一個零維應(yīng)用模式來解決非線性多項式方程：??????? (3）通過在“Phｙsiｃs”菜單、“Subｄomａinｓｅttingｓ"選項中的設(shè)定，使得表１中的每個子域都滿意該方程.注意左上角,這是以矢量符號給出的方程。在一維模式下,這個方程可以簡化為:? ?? (4）顯然，α、γ和β在一維簡化里是多余的。既然我們是求零維的根,全部公式４左側(cè)的系數(shù)都可以設(shè)為0。通過重新組合多項式,我們發(fā)現(xiàn)ａ=4和ｆ＝u＋ｕ＋２。注意我們是通過設(shè)定最大基元尺寸為1、將域離散化為只有一個基元的域,從而得到零維空間的！表1系數(shù)模式下的求根實例。文件名稱:ｒｏotｆiｎdｅr。mｐh．MoｄeｌNaviｇator選擇1D空間維數(shù)，COMＳOLMulｔiｐhｙsｉcs：PDEmｏdes：ＰＤE，cｏeffｉcienｔ形式Ｄraw菜單Sｐeｃiｆyoｂjeｃｔｓ:Ｌｉnｅ。Ｃoordｉnａt(yī)ｅs彈出菜單。x:01名稱:inｔeｒｖａl,完成Pｈysics菜單：Bouｎdａrｙsｅttｉｎgs選擇域：1和2（按住Ctrl鍵同時選中)選擇Nｅumanｎ邊界條件采納默認(rèn)設(shè)置q=０ｇ＝0,完成Physiｃs菜單:Ｓｕbｄoｍａinsｅttings選擇域：1設(shè)定:c＝0;a=4;ｆ＝ｕ^3+u^2+2；＝0應(yīng)用，選擇Init選項卡：設(shè)定ｕ（ｔ０)=－2,完成Meｓh菜單：Meｓhpａraｍeters設(shè)定最大基元尺寸1點擊remｅsｈ，完成Sｏｌve菜單：Soｌｖeｒpaｒａmeｔeｒs求解器選擇Stationaryｎoｎliｎeａr,求解,完成Pｏｓt－prｏcｅssｉｎgPointEｖalutｉon選擇邊界1，表達(dá)式:u，完成設(shè)定初始推測值為ｕ（ｔ）＝－2，下面我們來尋找最接近這個值的根.你可能對為什么設(shè)置ａ=４而不是把全部有關(guān)參數(shù)放進(jìn)ｆ中感到奇怪，那是由于對公式(4)右側(cè)有限元離散后得不到一個奇異剛度矩陣。后處理階段在輸出窗口中顯示出結(jié)果:Ｖalue：-2.7320５1，Expreｓｓioｎ:u，Bounｄａry:1。解析解得出的根為-1—，數(shù)值解能夠很好的與其吻合。依據(jù)代數(shù)中二次公式的解知道，另外一個根是—1＋，第三個根是1.回到子域設(shè)置中來，如果最初設(shè)置的推測值為u(t0）＝-0．5,則ＣOMSOLMultiｐhyｓics解出u=0．7６2051,又一個很好的近似數(shù)值解。如果u(t0)=１.2作為最初猜想時，得到ｕ＝１。該練習(xí)體現(xiàn)了非線性求解器和問題的兩個特性——（１）非線性問題可以有多個解；(２）最初推測值對于求解很關(guān)鍵。牛頓法通常收斂到最接近的解,但是在高階非線性問題中,這點可能不成立。這些特性在多維求解空間和“空間-時間”相關(guān)系統(tǒng)中同樣適用。COMＳOＬＭultiｐhysiｃｓ模型ｍｐｈ文件（ｒｏotfindｅｒ.ｍph)中包含了全部的MATLAB／COMSOＬSｃrｉpt源代碼和FEMLAB的擴展功能,以便再次恢復(fù)到ＦEMLAB當(dāng)前圖形用戶界面。該文件在網(wǎng)頁HYPERLIＮＫ”htｔp：/／ｅｙｒie。shef.ac.uｋ／femｌaｂ＂ｈｔtp：//eyriｅ。sheｆ.ａc.uｋ/ｆeｍlａb中可以得到.只需打開“Menu”菜單,選擇“Opｅｎmｏｄeｌｍ-fiｌe”選項，在彈出的“Opｅnfile”對話框中選中它,你就可以在“Ｓubdｏｍainsettings"中快速設(shè)定非線性函數(shù),給定初始推測值，然后用非線性求解器得到收斂解.但是如果函數(shù)中沒有線性項可以放在公式（4）左側(cè)怎么辦？例如，F(xiàn)ＥMLAB將taｎ（u）－u2+５=0放在公式(4）左側(cè)時，會生成奇異剛度矩陣。建議在子域中設(shè)定一個u的二次派生系數(shù)，c＝１。通過與Neｕmann邊界條件相結(jié)合，這個人為集中因素不能轉(zhuǎn)變基元中解為常數(shù)的事實，但是它可以避開剛度矩陣奇異。在一般模式下求根tan(u）—ｕ2+５=０中奇異剛度矩陣的問題可以通過一般模式來避開，它可以求解下式：????? ?(5)其中Г(u，ｕｘ）和（4）式中的系數(shù)項功能類似,但是被求解器處理方式不同。在系數(shù)模式中，認(rèn)為系數(shù)與ｕ無關(guān),除非使用數(shù)值雅可比矩陣，這會產(chǎn)生一些非線性依靠關(guān)系-—迭代次數(shù)增加。一般模式下構(gòu)建剛度矩陣時，精確雅可比矩陣會將Г和F對u進(jìn)行微分。通常一般模式比使用數(shù)值雅可比矩陣的系數(shù)模式需要的迭代次數(shù)少。下面使用的雅可比矩陣，不需要任何特殊處理來避開剛度矩陣奇異.總的來說,一般模式在處理非線性問題時比系數(shù)模式更加有用。從我的個人觀點看來,系數(shù)模式是ＣＯMSＯＬＭuｌtiｐhyｓiｃｓ的一個遺傳特性--MAＴLAＢ的偏微分方程工具箱，它在很多方面是ＦＥＭＬAB和ＣOMSＯLＭｕltipｈyｓiｃs的先驅(qū)，它廣泛使用了系數(shù)形式.此外，數(shù)值雅克比矩陣系數(shù)公式是一個長期存在的標(biāo)準(zhǔn)有限元方法，所以相對于其它代碼，它是一種格外有用的公式.表2給出了應(yīng)用一般模式的步驟——對我們剛才的步驟做了小小改動.表２一般模式下的求根實例。文件名：rｔfingen.ｍpｈMｏdeｌNaｖｉgaｔor1D，CＯMSOLMultiphyｓics:PＤEmｏdｅs,ｇｅｎerａl模式Oｐｔioｎs設(shè)定Ａxes/Ｇrid為［０，１]Dｒａw名稱：Inｔeｒvａｌ；Ｓtaｒｔ＝0;Stop=1Ｐhyｓicｓ菜單:BouｎdａｒySettｉnｇs設(shè)定兩個端點（域）為Neumａnn邊界條件Phｙsics菜單:SubｄｏｍａｉｎSettinｇs設(shè)定Г=0;=0；F=ｕ３＋u2—4*u+２初值u（t０）＝－0。5Meｓh模式設(shè)定最大基元尺寸，generaｌ＝1;點擊ＲeｍeｓhSolve使用默認(rèn)設(shè)置(ｎonｌineaｒsolver，ｅｘacｔJacobiａｎ）Pｏsｔ－Pｒocessinｇ經(jīng)過5次迭代，結(jié)果收斂。點擊圖片讀出u=0.7３20５1。轉(zhuǎn)變初始條件找到另外兩個根雖然這里給出了單變量簡潔函數(shù)求根的模板（ｒtｆｉｎdgeｎ．ｍｐh)，但是實際上ＭAＴＬAＢ有更簡潔的計算子程序,通過內(nèi)置函數(shù)ｆzeｒo和函數(shù)聲明來求根，現(xiàn)在ＣＯMＳOLMｕｌtiphysics圖形用戶界面可以使用幫助方程中的幫助變量來求解幾何約束,該變量稱為狀態(tài)變量。作為對一般模式求根模型的補充，使用狀態(tài)變量v依據(jù)表3的步驟來求解一個非線性方程。表３在ＯDE設(shè)置中求根Physics菜單：ODＥｓetｔｉｎｇs名稱：ｖ。方程：ｔａnh(v）—v^2＋5，完成Soｌve菜單:Solverｐaｒａmｅtｅrs選擇Staｔionａrynｏnlｉneａｒ求解器，求解，完成Post－ｐroｃｅssiｎg選擇Poiｎteｖａluａｔｉｏn,邊界２，表達(dá)式：ｖReporｔ窗口值:—2．00８８1９，表達(dá)式:v，邊界:2下一節(jié)我們將非線性求根方法應(yīng)用到一個常見的化工實例中——閃蒸,它用到了COＭＳOＬMｕltiｐhｙsicｓ圖形用戶界面的更多新特性。２。2求根：閃蒸實例化學(xué)熱力學(xué)中廣泛存在求根法的應(yīng)用實例，由于化學(xué)平衡和質(zhì)量守衡的約束條件很充分，再加上狀態(tài)方程,就產(chǎn)生了和問題中未知變量個數(shù)相同的約束條件。在本節(jié)中,我們以閃蒸為例來介紹簡潔的求根方法，構(gòu)建只有一個自由度的系統(tǒng)，這里以相分率作為未知變量。表4閃蒸單元的初始組分以及平衡時的安排系數(shù)K組分XI6５℃，3.4bar下的Kｉ乙烷0。00７91６。2丙烷0．１２８15.２i－丁烷０．０８４92。６n－丁烷0．26901.98i－正戊烷0．05890。91n-正戊烷0．１36１0．72己烷０．３15１0。２８液相碳?xì)浠衔锘旌衔锿ㄟ^閃蒸到65℃、3．4baｒ狀態(tài)。液相組成和每種組分閃蒸“Ｋ"值見表4.我們盼望知道閃蒸過程中氣相和液相產(chǎn)物的成分和液相的比例。表4給出了原始成分。組分i的質(zhì)量平衡滿意以下關(guān)系? ?? （6)Ｘi是閃蒸前液體的莫爾分?jǐn)?shù),xi是閃蒸后液體莫爾分?jǐn)?shù),yi是蒸發(fā)產(chǎn)物莫爾分?jǐn)?shù),是液相產(chǎn)物與總供應(yīng)摩爾流率的比值。平衡系數(shù)的定義為：Ｋｉ＝ｙi／xi。將該方程消去質(zhì)量平衡中的xi,得到一個關(guān)于ｙi和Xｉ的方程：???? ?? （７）由于各yi相加為１，所以我們有的非線性方程：??????(8）其中n是組分?jǐn)?shù)量。該函數(shù)可用ｒｏot（s)求解，將計算結(jié)果回代就可以得到全部產(chǎn)物的莫爾分?jǐn)?shù)。牛頓-拉夫遜方法需要知道當(dāng)前導(dǎo)數(shù)來決定計算方向,在COMＳOLＭｕltiphｙｓiｃs中將解析計算該值。通過下式很容易得到牛頓-拉夫遜迭代方法:??? ??(９)現(xiàn)在再回到ＣOＭＳOＬMｕlｔiphysics求根計算。作為一個練習(xí)，我們用一般ＰＤＥ模式來計算求解。我們可以只載入ｒootfindｅr。mph或ｒtfindgen.mｐｈ并加以修改來計算這個問題，但是熟識ＣＯMSOＬMulｔiphｙsics功能也是格外重要的目的.表５閃蒸實例ModelNaviｇaｔｏr選擇１D,COMＳOLMultｉｐｈyｓics:ＰDＥmｏdes:PＤＥ，geneｒal形式Draｗ菜單Specifyobｊｅcts:Linｅ，Ｃoｏrdiｎａｔes彈出菜單,x：01名稱：ｉnｔerｖal，完成Ｏptiｏｎｓ菜單：Consｔaｎts輸入表4中數(shù)據(jù)X１,０．0079Ｘ2，0．1281X3,0．０849X４,0.２6９0X5,0．0589X6,0.1３6１X７，0.3151Oｐｔiｏnｓ菜單:ＳcalarExpressｉons為（8)式中右側(cè)各項定義表達(dá)式t１-Ｘ1／（１－u＊（1—１/K１））t2－X2/(1－ｕ*(１－1／Ｋ２）)ｔ3-X３／（1－u＊（１—1/K3))ｔ4-X４／（1—u＊（１－１/K４))ｔ５—Ｘ5／（1-ｕ＊（1-1/K5）)t6-Ｘ６/(１—ｕ*（1-1/K6）)t7－X7／(1－u＊(１－1／Ｋ7)）Physiｃｓ菜單：Boｕnｄarｙsettｉｎｇs選擇域：１和２(按住Cｔｒｌ鍵)選擇Ｎeumanｎ邊界條件采納默認(rèn)設(shè)置q=０,ｇ=0完成Ｐhｙｓｉcｓ菜單:Ｓｕbdomaｉnseｔtinｇs選擇域:１設(shè)定F=1＋t１＋ｔ2+…+t7；=0選擇Iｎiｔ選項卡，設(shè)定u（t0）=０.5Mesh菜單：Mｅsｈparａmｅｔeｒs設(shè)定最大基元尺寸為1點擊Rｅmesｈ，完成Solｖe菜單：Solverｐaraｍeteｒｓ選擇Stationaｒynｏｎｌineａr求解器求解,完成Post—ｐｒoｃeｓsｉｎg：PointEvaｌｕation選擇邊界：1,表達(dá)式:u,完成Reｐoｒt窗口,值：０。45８５0９，表達(dá)式：ｕ啟動COＭＳOLＭｕlｔｉｐhyｓｉcs,打開模型導(dǎo)航欄.如果你已經(jīng)打開COMＳOLMultｉｐhyｓics對話框,那么將你的工作空間保存為mph文件或者將命令保存為m文件，然后打開“Fｉle”菜單選擇“Ｎeｗ”。依據(jù)表5的步驟建立閃蒸模型。注意本例中有兩個新增內(nèi)容--“Optｉons:Constａｎｔs”和“Optｉons：Expresｓｉons”.在這里可以定義常數(shù)，然后在ＣOMＳOLＭｕltiｐhyｓｉｃｓ中任何可以使用純數(shù)字的輸入框中使用定義的常數(shù)。表達(dá)式和常數(shù)類似，可以在任何ＣOMＳOLＭulｔiｐhｙｓics數(shù)字輸入框中使用,但是不同之處在于表達(dá)式依靠于因變量。定義的常數(shù)和表達(dá)式同樣可以在后處理過程中使用。后處理過程顯示t1和t2為：Ｖalue:—0．01３865,Ｅxｐrｅｓｓiｏｎ:ｔ1，Bouｎdaｒy：1Vａlue:—０．２0344１,Eｘpｒeｓsｉoｎ:t2，Ｂouｎｄary：１另外求解器日志中常常包含有用的信息。打開“solｖeｒ"菜單，在最底部選擇“ＶｉｅｗＬog”，會彈出求解器日志對話框。它顯示了求解器何時開頭計算，執(zhí)行了什么命令，如何從初始狀態(tài)得到計算結(jié)果。這里經(jīng)過三次迭代，肯定誤差達(dá)到１0—9。如果你的解不是不收斂或收斂很慢的話,很容易得到這樣的信息。練習(xí)計算方程的根.由于該函數(shù)是三次多項式，所以存在解析的無理數(shù)解，。ueｘｐ（u）—1計算方程的根。這是個超越方程，意味著不存在解析的無理數(shù)解.如果你使用系數(shù)模式,設(shè)ｃ=1來幫助收斂。3．方法２:步進(jìn)法數(shù)值積分?jǐn)?shù)值積分是數(shù)值分析的支柱.在有計算機之前,科學(xué)計算的首要工作是制作特殊方程手冊,其中幾乎包括了全部特殊常微分方程的解。那么用什么方法計算得到的？就是用一維數(shù)值積分。一維數(shù)值積分分為兩種：初值問題(IVP)和邊界值問題（ＢVＰ）。我們下一節(jié)再介紹邊界值問題.最容易積分的是初值問題，由于全部的初始條件都在一點定義,可以依據(jù)一階導(dǎo)數(shù)直接步進(jìn)積分。顯然，如果常微分方程是一階的，例如： ? ??（10)那么當(dāng)△ｔ→0時,（1０)中其次項嚴(yán)格成立。它滿意歐拉法的條件,是積分一階常微分方程最簡潔的方法.對于一維情況，可以簡潔地依據(jù)ｆ在（xn,yｎ)點的導(dǎo)數(shù)向前步進(jìn)積分，這里n是指積分過程中的第n次離散化步驟。所以，??? ??(１１）這里假設(shè)導(dǎo)數(shù)變化不超過步進(jìn)量ｈ，這只對線性函數(shù)才嚴(yán)格成立。對于任何有曲率的函數(shù),該假設(shè)都不成立。下面考慮在圖1中，如果步進(jìn)量h取得過大會造成什么錯誤.顯然，歐拉法需要改進(jìn)的一點就是增大步進(jìn)量，由于它需要小步進(jìn)量以保證精準(zhǔn)性。歐拉法也被稱作“一階"精準(zhǔn)性，由于誤差只能降低到一階h的程度.圖1數(shù)值積分中的歐拉步進(jìn)龍格庫塔方法所以如果我們盼望增大步進(jìn)量，那么就需要一個“更高階的方法"，隨著步進(jìn)量的減小,誤差快速降低.k階方法誤差大約在hk量級。假設(shè)我們忽視了曲率,可以通過計算斜率f（x）在ｘｎ和xn+１之間幾個中間點的值來計算ｙ（ｘ）的曲率?？梢允紫纫罁?jù)初始導(dǎo)數(shù)計算間隔中點的值，然后再用整個間隔中點的微分值來獲得二階精準(zhǔn)性。? ?? ?（１2)這樣我們通過對兩個函數(shù)的計算又得到一階精準(zhǔn)性.所以,如果使用一階方法，Ｎ次計算得到Ｏ（1/N）誤差,而二階方法2Ｎ次計算得到O(1/４Ｎ２)誤差，如果使用一階方法的話,這需要Ｎ2次計算。更高階龍格庫塔方法我們可以做的更好嗎?顯然，如果用三個中間點可以獲得三階精準(zhǔn)性，用四個中間點可以獲得四階精準(zhǔn)性等等。不過更高階的方法需要更多的編程工作，所以也要考慮我們的時間.但是從本質(zhì)上講，我們計算函數(shù)的k階導(dǎo)數(shù)不肯定光滑?？赡茉谠黾印敖凭珳?zhǔn)性”的同時,高階導(dǎo)數(shù)的整體誤差會變得很大,如果是這樣，每次連續(xù)步進(jìn)都可能導(dǎo)致誤差的急劇增大。這表明高階方法沒有低階方法穩(wěn)定。常微分方程通常選用四階龍格庫塔方法積分，這樣程序比較緊湊,精準(zhǔn)性高，也有很好的穩(wěn)定性。其它方法在數(shù)值積分編程的時候，還需要注意兩個問題:數(shù)值穩(wěn)定性：即使使用高階精準(zhǔn)性的方法，你的積分結(jié)果與檢驗值偏差仍然很大，這可能就是由于你的數(shù)值方法不穩(wěn)定。你可以通過降低步進(jìn)量來獲得數(shù)值穩(wěn)定性，但是這也意味著更長的計算時間。如果你的計算量很大，計算時間太長,可以選用半隱式方法,例如猜測校正法等。剛性系統(tǒng)：物理機理起作用時，剛性系統(tǒng)通常會有兩個完全不同的長度或時間尺度。剛性系統(tǒng)可能會消滅上面提到的“系統(tǒng)非穩(wěn)定性”現(xiàn)象，或者非物理現(xiàn)象的振蕩。可以參考Geａr[2］的書來解決這個問題。３.１數(shù)值積分簡潔實例高階常微分方程通常通過降階和步進(jìn)法求解，考慮如下方程:???? ? （13)除非q(ｘ）和r(x)是常數(shù),那么你很幸運能夠從大多數(shù)分析方法教材中找到答案。也有一些特殊的q（x)和r(x)可以求得解析解,但是最好能計算出各種情況下的數(shù)值解。為什么呢？由于為了搞清楚y（x）,你需要畫出曲線,得到解析解時你也需要花費肯定的精力在作圖上.那么應(yīng)該怎么做?首先將上面的二階系統(tǒng)降低為兩個一階系統(tǒng)：上式的常微分方程能夠仿照（11)和（１2）進(jìn)行時間步進(jìn)法數(shù)值積分。舉個簡潔例子：?? ????（１4）降階產(chǎn)生了兩個一階常微分方程：??? ? ?(15）假設(shè)初始條件為u１=1，u2=０，可以建立零維空間系統(tǒng)來積分這一對常微分方程。打開ＣOMＳOＬMｕltiphysｉcs，等待模型導(dǎo)航欄窗口。如果已經(jīng)打開一個CＯMＳOＬMulｔiphｙｓｉｃｓ窗口，把工作空間保存為mｐh文件或者把命令保存為ｍ文件，然后在“Fｉｌe"菜單中選擇“Neｗ”選項。依據(jù)表6的步驟進(jìn)行設(shè)置。表6時間步進(jìn)積分的簡潔實例ＭｏｄｅlＮaｖigａt(yī)or選擇１D,COＭSOＬMｕltiphｙｓicｓ：PＤEｍoｄeｓ:ＰDE，genｅｒal模式，插入因變量：u1u２Drａw菜單Speｃｉｆyｏｂｊｅcts:Lｉｎe。Coordinaｔes彈出菜單,x:01名稱:intｅrvａl,完成Physｉcｓ菜單：Bｏuｎｄaｒysetｔinｇs選擇域：１和２（按住Ｃｔrl鍵)選擇Ｎｅｕmann邊界條件保持默認(rèn)選項ｑ＝０，ｇ=0完成Phｙｓiｃs菜單：Subdomaｉnsｅttiｎｇs選擇域：１設(shè)定Ｆ１=—ｕ２,F2=u1選擇Iｎiｔ選項卡，設(shè)定u1（ｔw)=1Mｅsｈ菜單:Meｓｈpａrameterｓ設(shè)置最大基元尺寸為1點擊Ｒemesh，完成Ｓoｌｖe菜單：Sｏｌｖｅrpａｒameters選擇ｔｉｍｅdｅpeｎdent求解器,在Geneｒａl選項卡中輸入Ｔimｅｓ：ｌinspaｃe(0,２*ｐｉ，5０)求解，完成Post－pｒoｃessing:Crosｓ-secｔioｎｐｌoｔparａｍeｔersＰｏinｔ選項卡，接受ｕ１的默認(rèn)設(shè)置Geｎerａl選項卡，接受全部時間的默認(rèn)設(shè)置完成這個例子給了我們兩個獨立變量u１、u２和空間坐標(biāo)x。注意“Sｕbdoｍaiｎsｅｔtｉnｇs”窗口左上角方程中的矢量標(biāo)志.由于是矢量變量,全部輸入選項卡都是矢量(F）或者矩陣（dａ)輸入。MＡTＬＡB命令linspace（0，2*pｉ，５0）生成50個基元的矢量，給出從0到２π的數(shù)據(jù),u1(t=2π）=1.0０441４.假設(shè)解析解是ｕ１（t=2π)＝1,結(jié)果不夠精準(zhǔn)（0．４％誤差)。之前FＥMＬAB允許用戶在MＡＴＬAＢ和FEＭＬAＢ中選擇針對時間積分的預(yù)置的求解器。也允許用戶在求解器變量對話框的通用選項卡中調(diào)整容錯值(相對和肯定)。將相對誤差調(diào)整為0。0０１，肯定誤差調(diào)整為0.00０1,得到一個相對較好的計算結(jié)果u１（t＝2π）＝0.998０27.注意累積的全局誤差與求解方法精準(zhǔn)性０。0０1的量級全都.圖２一個周期的u１(ｔ)圖3一個周期的ｕ２(t)圖2和圖3表明如果設(shè)定足夠小的步進(jìn)時間，F(xiàn)EＭＬAB能夠完成高精度的正弦和余弦函數(shù)數(shù)值積分。雖然我們一般認(rèn)為正弦和余弦是“解析函數(shù)”，但是通過比較，還是能清楚地看到解析函數(shù)和需要數(shù)值積分的函數(shù)是似是而非的——沒有比Ｂｅｓsel函數(shù)和橢圓方程更有解析性的函數(shù)。練習(xí)使用物理場菜單中的常微分方程設(shè)置和相同的初時條件,試著積分方程（１5）,假設(shè)狀態(tài)變量為v1和v２?；跁r間的常微分方程和代數(shù)方程的唯一區(qū)分就是必須使用符號代替ｄv１/dt和v1ｔ等。４．管式反應(yīng)器數(shù)值積分本節(jié)介紹在化工管式反應(yīng)器設(shè)計時，如何同時求解一階非線性常微分方程組。通常管式反應(yīng)器設(shè)計的關(guān)鍵要素是反應(yīng)器的長度。一個氣態(tài)酒精脫水管式反應(yīng)器，工作在2bar和１５0℃下。反應(yīng)方程式為:已有試驗表明反應(yīng)速率可以表示如下，其中CA是酒精濃度(mｏl/L)，R是酒精的消耗速率（ｍoｌ/s／m3)：反應(yīng)器直徑為0.0５ｍ，入口酒精流速為１0ｇ／s。問：如何確定反應(yīng)器的長度，進(jìn)而獲得各種酒精轉(zhuǎn)化率？我們盼望確定酒精出口摩爾分?jǐn)?shù)分別為0．５,０．4，0.３,0．２和0．1時的反應(yīng)器長度?；ぴO(shè)計理論假設(shè)反應(yīng)熱很小，反應(yīng)器內(nèi)為柱狀流，抱負(fù)氣體,可以確定反應(yīng)流體由四個常微分方程掌握，其中獨立變量為CＡ，CＷ（水濃度),Ｖ（速度）和ｘ（反應(yīng)器長度）: ?? ?? (16)最后一個方程說明表面速度使得反應(yīng)器長度和反應(yīng)物在反應(yīng)器內(nèi)的存在時間之間存在一個平衡.初始入口條件為：???? ? (１7)方法從初時條件和化學(xué)計量系數(shù)可以看出,CW＝CE(酒精濃度），由于溫度和壓力設(shè)為常數(shù)，所以濃度C也為常數(shù),可以看作是抱負(fù)氣體，有：? ?????（１8）初始進(jìn)口速度Ｖ0可以由給定流速、入口密度(酒精分子量為46kg/kmol）和管的橫截面積計算出。該方程需要對時間t數(shù)值積分到需要的酒精摩爾分?jǐn)?shù)。使用較為簡潔的歐拉方法或者龍格庫塔方法積分。讀者也許注意到,積分ＣA可能不需要其它變量，但是CＷ，V和x與時間ｔ嚴(yán)格相關(guān)。但是由于我們需要求出對應(yīng)摩爾分?jǐn)?shù)的x值，所以最好同時求解四個變量。部分結(jié)果計算得到的數(shù)值解為：? ?? ???(19）其中CA/Ｃ如圖4所示。圖4酒精濃度隨時間t變化曲線ＣOMSOLＭultｉphysics計算我們盼望再次建立虛擬的零維模擬空間，這次使用四個獨立變量。打開CＯMSOＬMultiphyｓｉｃs，等待模型導(dǎo)航欄窗口。依據(jù)表7的步驟建立管式反應(yīng)器模型。該模型給出四個獨立變量u1、u２、ｕ3、u４和一個空間坐標(biāo)x。由于變量比較多,建議用變量名稱進(jìn)行常數(shù)設(shè)定.進(jìn)一步來說,由于使用了速度定律，用矢量表達(dá)式更為便利.表７COMＳＯLMｕlｔiphyｓicｓ中管式反應(yīng)器設(shè)計建模步驟ＭoｄelNａｖｉｇatoｒ選擇１D空間維數(shù),COＭSＯLMｕltｉphysｉcｓ：PDEmodes：PDＥ，ｇeneral形式設(shè)定獨立變量:u1u2u3ｕ4選擇Elｅmｅnt:Lagrangｅ—Lineａr，完成Ｄｒａw菜單Ｓpeciｆyoｂjects：Lｉne，Cｏordinaｔes彈出菜單，x：01名稱：intｅrval，完成Opｔiｏｎｓ菜單:Ｃonstanｔs如下填寫變量表格：名稱EｘprｅssioｎP20０00０T423R8３14ＭＭ46Fｌowrate0。01Dia0。0５CP/（RT)arｅapi＊Dｉａ^2／4rhoMＭ＊CvｅｌFlowraｔe/rho/ａera完成Options菜單：ScａlarExｐｒessions定義ratｅ=5２.7＊ｕ1^2／(１+０.０１３/u１）Physiｃs菜單：Ｂouｎdaｒysettinｇs選擇域:1和2（按住Cｔrl鍵）選擇Neｕmａnｎ邊界條件默認(rèn)設(shè)置ｑ=0g＝0,完成Phyｓｉcs菜單：Sｕｂdomainsｅtｔinｇs選擇域１F選項卡；設(shè)定F1=-ｒate*（1+u１/C);F2=rａｔe＊(1-u２/C)；Ｆ3=ｒaｔe*u3/Ｃ；F4=u3/ＣＩnit選項卡；設(shè)定ｕ1（ｔ0）＝Ｃ;ｕ３(t0)=veｌ，完成Ｍｅsh菜單:Ｍｅsｈpａramｅters設(shè)定最大基元尺寸為1點擊Ｒemｅsh，完成Ｓolve菜單:Solvｅrｐarameterｓ選擇timedependenｔ求解器，在Gｅneｒaｌ選項卡中輸入Tiｍｅｓ:linspａｃｅ（0，10，１00）求解，完成Post—ｐroｃessingCｒosｓ－ｓecｔionｐｌotpaｒａｍｅtｅrsＰｏint選項卡,接受u１的默認(rèn)設(shè)置Ｇeneral選項卡，接受全部時間的默認(rèn)設(shè)置完成試著在整個時間范圍內(nèi)畫出u１,ｕ2,u3和u4的曲線，與圖4是否保持全都？與之前計算結(jié)果是否相符？練習(xí)計算以下方程組中（ｘ=1)的值，并畫出隨x變化曲線，ｘ取0到３區(qū)間。連續(xù)攪拌反應(yīng)器中的一階可逆反應(yīng)系統(tǒng)可以用線性常微分方程組表示。例如,考慮如下異構(gòu)反應(yīng):正向反應(yīng)速率為k1和k3,相應(yīng)逆向反應(yīng)速率為k２和ｋ4。一階動力學(xué)常微分方程組為:?????（２0）你也許會感到奇怪，由于該方程組為線性,它有通用的解析解。雖然是通用解析解,但是對動力學(xué)系統(tǒng)的反映并不多。假設(shè)剛開頭時，CA＝１,ｋ1＝1Hz,k2=0Ｈz，ｋ3＝2Hz,ｋ４=３Ｈｚ。針對該初值問題，畫出濃度隨時間變化的曲線。5。方法3：常微分方程數(shù)值積分前一節(jié)介紹了步進(jìn)法數(shù)值積分，通常也稱作“時間步進(jìn)",但是在反應(yīng)器設(shè)計中，多是空間積分問題。步進(jìn)法是挨次求解未知項，而其它常見積分方法是同時求解某一網(wǎng)格點的全部未知獨立變量。步進(jìn)法只能求解初值問題（ＩＶP),初值個數(shù)必須嚴(yán)格等于常微分方程組的階數(shù)。但是對于二階或者更高階的系統(tǒng),可能會用到其次類邊界條件-—邊界值問題，在一維情況下，這些邊界條件只消滅在域的起點和終點，這就是兩點邊界值問題。步進(jìn)法也能夠牽強求解邊界值問題——人為設(shè)定初值，通過嘗試和誤差修正找到滿意邊界值問題的初始條件。對于更高維數(shù)的偏微分方程，邊界值問題發(fā)生在整個域的每一個邊界上。有限元方法的一個最主要優(yōu)點就是能夠很容易地求解兩點邊界值問題。例如，一維反應(yīng)和集中方程：????? ??(2１）這里u是組分濃度，D是集中系數(shù),L是域的長度,R(u）是反應(yīng)消耗速率,ｘ是無量綱空間坐標(biāo)。如果未知函數(shù)u(ｘ）在x=xj=j△x網(wǎng)格點能夠用不連續(xù)值uj=u（xj）近似,那么依據(jù)中心差分，方程變?yōu)椋? ???（22）其中Ｒj=R（uj）,Mij是對角元素為－２，上、下對角元素為１的三對角矩陣：? ?????（23）。通過矩陣轉(zhuǎn)置和對ｕiｎ進(jìn)行積分可以求解該方程,這里n指第n次推測:?????? (24）Rj=Ｒ（uiｎ－1)。無論是初值問題還是邊界值問題，都可以將(23)式中矩陣Ｍ的行向量變得符合邊界條件.(2３）式假設(shè)在ｘ=0和ｘ=１處都有ｕ＝0。這是狄利克雷邊界條件，也是有限微分法的自然邊界條件-—說它自然是由于如果在設(shè)定邊界條件時沒有轉(zhuǎn)變（23)式中的行向量,那么就會消滅這樣的邊界條件。下面介紹如何用CＯMSOLMｕltipｈysics在一維域上求解方程（21),對于一階反應(yīng)R（u)=ku,典型的無量綱變量Dａｍkｏhler數(shù)為:?????(２5)邊界條件為ｘ=０時u＝1,在u=1處沒有流量。下面結(jié)合MATLＡＢ來做這個練習(xí),進(jìn)而了解ＣＯMＳOLMulｔｉphyｓiｃs如何表示計算數(shù)據(jù)和模型輸出的結(jié)構(gòu).在wｉndows中,結(jié)合MＡTLAB的COMSOＬＭｕltiphyｓics是一個可選的桌面圖標(biāo)，在UＮIX/ｌinux中,該功能可以通過在命令欄中執(zhí)行以下語句實現(xiàn)：Coｍsｏlmaｔlａbpath－ｍlnodesktop－ｍlnosplaｓｈ其中“mａt(yī)lａb”告知femlab打開一個ｍaｔlaｂ命令窗口,“ｐａt(yī)ｈ"依據(jù)COMＳOＬ命令庫建立maｔlab命令窗口。首先運行CＯＭSOLMｕｌｔｉｐhｙｓｉｃs和模型導(dǎo)航欄窗口，然后依據(jù)表8中的步驟建立模型.該實例給了一個獨立變量ｕ和一個一維空間坐標(biāo)x。ｈ和ｒ是系數(shù)模式中狄利克雷邊界條件的兩個句柄。如果你盼望邊界速度u為常數(shù)Ｕ0,可以通過設(shè)定ｈ＝1,ｒ=U０來實現(xiàn)。點擊三角形按鈕生成默認(rèn)網(wǎng)格（１５個），然后點擊兩個嵌套三角形按鈕加密網(wǎng)格（30個).表８反應(yīng)－集中系統(tǒng)中兩點邊界值問題的常微分方程實例MｏdeｌNaｖigaｔor選擇１D空間維數(shù)，CＯＭSOLMuｌtｉphysics：PDEｍｏdｅs:PDＥ,coｅｆｆicieｎt形式設(shè)置因變量:u選擇Ｅlement：Ｌaｇraｎge-Linｅａr，完成Draw菜單Speｃiｆyobjects:Ｌinｅ，Cooｒｄinatｅs彈出菜單,ｘ：０1名稱:ｉnteｒvaｌ，完成Phyｓiｃs菜單：Ｂｏｕndaryｓｅttｉngs選擇域1選擇Diｒicｈｌet邊界條件，設(shè)定ｈ=1;ｒ=1選擇域2選擇Ｎeuｍann邊界條件,完成Phｙｓｉcs菜單：Suｂdomａiｎsettiｎｇｓ選擇域1設(shè)定Ｃ＝—１;f=0。83３*ｕ;=0選擇Ｉｎｉt選項卡;設(shè)定u（ｔ０）＝1-x完成Ｍeshinｇ點擊三角形按鈕進(jìn)行網(wǎng)格劃分Solve菜單:ＳｏlveｒｐａrametｅrsStationarylineａｒ求解器默認(rèn)設(shè)置,完成Ｇenｅｒal選項卡，設(shè)定解的形式為“Ｃｏｅfｆicienｔ"求解(＝）Pｏｓt－processｉng：Pｏｉntｅｖalution選擇邊界２,輸入表達(dá)式ｕ。Rｅｐorｔs窗口顯示:值：０.86116７,表達(dá)式：ｕ，邊界:2依據(jù)計算結(jié)果可以畫出圖５，顯然邊界條件需要滿意：x＝0處u=１，在ｘ=１處斜率減為零。但是COＭSＯLMuｌtipｈｙsics有沒有依據(jù)我們設(shè)想的求解呢？圖5穩(wěn)態(tài)下的濃度曲線現(xiàn)在用穩(wěn)態(tài)非線性求解器再計算一次。首先查看求解器菜單中的日志，迭代13次收斂.你注意到最終值從0。8６降到0.６9了嗎?為什么會這樣呢？線性（系數(shù)型）求解器只在初始時刻u（t0)=1—x時計算一次R（u）,所以方程（2４）只需要一次迭代。而非線性求解器在每次迭代時計算一次Ｒ(u）,計算時使用上次迭代的u值.所以非線性求解器在向收斂靠近時，經(jīng)過幾次迭代可能完全“忘記”了初時推測。檢驗一下這種解釋對不對。轉(zhuǎn)變初值將會轉(zhuǎn)變穩(wěn)定的線性解.返回子域設(shè)置對話框，設(shè)定初值為ｕ（t0)＝1,最終得到的u（x=1）是多少？再試一下穩(wěn)態(tài)非線性求解器，能得到和其它初值一樣的結(jié)果嗎？該例說明為方程選擇正確求解器的重要性.如果函數(shù)f依靠于任何非獨立變量，就應(yīng)該選擇非線性求解器。線性求解器更快,但是它假設(shè)偏微分方程的系數(shù)不依靠于非獨立變量ｕ(否則為非線性問題）。如果不確定，還是應(yīng)該用非線性求解器。別忘了,(21）式滿意R（ｕ）=ｋu，是線性問題，但是COMSＯLＭultiphysicｓ使用非線性求解器才得到正確解!收斂比較慢也是模型形式導(dǎo)致的結(jié)果——選中雅克比求解器選項，非線性求解器只需要兩次迭代就可以得出正確結(jié)果.我們一開頭認(rèn)為（2２)式是該問題有限微分矩陣方程，但是實際上最后用(２４）式描述CＯMSＯLMuｌtiphｙsics的有限元問題。由于我們這里使用拉格朗日線性微元,在這種特殊情況下有限元和有限微分矩陣全都。依據(jù)這一點,我們看一下MＡTＬＡＢ/ＣＯＭSOLSｃｒｉpｔ如何計算該ＣOＭＳOLＭｕｌtｉphysics問題.打開“Fｉle”菜單,選擇“EｘporｔFEMStｒuctuｒｅasｆem”。這就把現(xiàn)在的計算結(jié)果轉(zhuǎn)化為MATＬAB/ＣOＭSOLScｒiｐｔ工作空間的數(shù)據(jù)結(jié)構(gòu),然后用MAＴLAＢ/ＣOＭSOLScript預(yù)置的函數(shù)和命令來計算它。在MAＴLＡB/CＯMＳOLScｒｉpｔ工作空間輸入以下命令：〉〉x=feｍ.mesｈ。p；〉＞u＝fｅm．sｏl.u;＞〉plｏt(x，ｕ)將會彈出一個包含網(wǎng)格節(jié)點u值的MAＴLＡB／CＯMＳOLSｃript圖形窗口。不要懷疑你的圖形失真。這是由于CＯＭSＯLMultiphyｓics存儲網(wǎng)格點和對應(yīng)結(jié)果的方式是為了使矩陣變得稀疏和緊湊?？梢杂肕ATLAB排序命令對網(wǎng)格點進(jìn)行排序，進(jìn)而使圖形有意義。在ＭAＴLＡＢ工作空間輸入以下命令：〉>［ｘx，idx］＝sort（ｘ）；>＞ploｔ（xｘ,ｕ(iｄx)）最后考慮MAＴLＡB對有限元矩陣的訪問。Fem文件結(jié)構(gòu)沒有包含有限元剛度矩陣,但是它包含重建這個矩陣的ＦEMLAＢ函數(shù)的必要信息。這對有限元方法格外重要,這個FＥMLAＢ函數(shù)是ａsseｍble,輸入以下命令：>＞［Ｋ，L,Ｍ，N］=assembｌｅ（fｅm）；>＞K’／15這張圖包含上次COMSＯＬMｕltipｈysics的計算結(jié)果，類似圖5。實際上,我們之所以能夠快速搞清楚ｆem求解結(jié)果的結(jié)構(gòu)，是由于這是只有一個獨立變量的一維問題。多變量和多維產(chǎn)生的網(wǎng)格和結(jié)果的數(shù)據(jù)結(jié)構(gòu)只能用COＭSOＬMuｌtiphｙｓiｃｓ的工具／函數(shù)來讀取。圖６給出了由MAＴLAＢ命令spｙ(K’)產(chǎn)生的松散結(jié)構(gòu)。很有必要將它與（２3)式相互比較,由于能夠清楚地看到具有統(tǒng)一網(wǎng)格劃分的一維拉格朗日線性基元和生成矩陣方程的有限微分法是類似的。圖6由ＭATLAB命令sｐy（)產(chǎn)生的松散結(jié)構(gòu)下面看一下MＡＴＬAB中矩陣的松散結(jié)構(gòu),全部的元素只有1，-2和1，位于三條對角線上。這是有限元法的剛度矩陣，對于未知類型，它等于(23)式.如果返回“Ｓuｂｄｏmａｉnsｅtｔings”對話框的“ｅlｅment”選項卡,選擇拉格朗日二次元素，再次求解并導(dǎo)出FEＭ，生成上面的K，你會發(fā)現(xiàn)雖然矩陣也很疏松,但是和拉格朗日線性元素完全不同。練習(xí)偏微分方程的系數(shù)形式有一項αｕ，令α=0．8３3，f=0，再次計算以上“反應(yīng)—集中"的例子，比較穩(wěn)態(tài)線性和非線性求解器得出的結(jié)果。你能解釋為什么這樣的表達(dá)式會導(dǎo)致這樣的結(jié)果嗎？這樣的表達(dá)式對剛度矩陣Ｋ又有什么影響?如果把Da選為某特定值,例如Da△x２=2，使得對角元素幾乎為零，你認(rèn)為求解時會消滅什么困難？6。方法４:線性系統(tǒng)分析MＡTLAＢ和COＭSＯＬMｕｌｔiphysics的核心是線性系統(tǒng)分析.本節(jié)將簡要回顧線性算法理論—-本科生工程數(shù)學(xué)中的“矩陣方程”就是這類典型問題。好消息是不需要你自己動手編程處理矩陣，這就是MATLAB的用途：供應(yīng)一個工程矩陣計算子程序庫的用戶界面。應(yīng)當(dāng)注意到，對于本章的例子,ＣOMSＯLSｃｒiｐt也能很好地實現(xiàn)這個功能.以前科學(xué)計算主要集中在矩陣計算的高效和松散方法上。Golub和ＶａnLoａn[3］的書是一本很好的專家級矩陣計算指南,但是對于ＭAＴLAＢ入門水平，Hａnsｅｌmaｎ和Littlｅｆieｌd[4]的書是個不錯的選擇。簡潔來講,標(biāo)準(zhǔn)的矩陣方程如下：???? （26)這里有N個與Ｍ方程有關(guān)的未知變量ｘj。系數(shù)ａiｊ和右側(cè)常數(shù)項ｂi都是未知數(shù)。在工程應(yīng)用中，通常將模型轉(zhuǎn)化為這樣的線性方程系統(tǒng).例如，質(zhì)量和能量守恒通常都會生成這樣的線性方程系統(tǒng)。求解可能性當(dāng)N=M時,約束條件和未知數(shù)相同，所以可能求出唯一解ｘj。但是如果一個或者多個方程彼此線性相關(guān)(行退化）,或者全部的方程都只含由某些特定變量組成的未知數(shù)（列退化）,就不能得出唯一解。對于稀疏矩陣，行退化和列退化相同，退化方程會導(dǎo)致奇異。但是從數(shù)值求解講，至少還有兩點會導(dǎo)致錯誤：如果方程彼此不是嚴(yán)格線性相關(guān),一些方程也可能在計算過程中由于舍入誤差而導(dǎo)致變得格外接近線性相關(guān)。求解過程中累積的舍入誤差可能會“沉沒”真實解，這在N比較大的時候容易發(fā)生。這種情況下，計算過程沒有出錯，但是計算結(jié)果會不滿意原始方程。線性系統(tǒng)指南沒有“典型”的線性方程系統(tǒng),但是一個粗略的想法就是舍入誤差變得不行忽視：當(dāng)Ｎ在20—５0之間時,可以用單精度常規(guī)方法求解,不需要對以上兩種錯誤進(jìn)行修正.當(dāng)N在幾百數(shù)量級時可以用雙精度求解。當(dāng)N在幾千數(shù)量級時，如果系數(shù)矩陣稀疏（幾乎都是零)，由于系數(shù)特性可以求解.對于稀疏矩陣,ＭAＴLAＢ有特殊的數(shù)據(jù)形式和一套針對性的函數(shù).但是在工程實際和物理過程中,有很多問題本身就是奇異或者格外接近奇異的。你可能會發(fā)現(xiàn)N=１0也很難求解.用奇異值分解法有時可以解決奇異值問題，將其分解為非奇異值。數(shù)值線性代數(shù)的常規(guī)任務(wù)方程（2６)可以寫為緊湊矩陣方程形式:???? ? ?（27)求解未知矢量x，這里Ａ是方陣，ｂ是已知矢量。A為常數(shù),求解多個ｂ矢量時的解。計算方陣A的逆矩陣A—１.計算方陣A的行列式。如果Ｍ〈N，或者Ｍ＝Ｎ但是方程退化，這時有效方程數(shù)比未知量少-—欠定系統(tǒng)。這種情況下可能沒有解或者有多個解。解空間由特殊解xp乘以任何線性相關(guān)的Ｎ－M維矢量構(gòu)成，稱為矩陣Ａ的化零空間.我們的任務(wù)是通過奇異值分解找出解空間。如果M〉N，通常沒有符合(26）的矢量解x。常常會消滅這種過定系統(tǒng),我們的任務(wù)主要是尋找最接近于滿意方程的折衷解。通常接近程度以（２６）方程左右兩側(cè)差值的“最小二乘法"來表示。MＡＴLＡB矩陣計算通過ｉｎv(矩陣）命令可以很容易獲得逆矩陣。矩陣方程的解以矩陣除法運算符“\”表示:>＞A=［3—10；—16—2；0－２10］;〉〉B＝[1；５；２6］；〉〉Ｘ=Ａ\BＸ=１.０0002．00０03。００００行列式行列式多用在穩(wěn)定理論和評估矩陣奇異性上。為什么需要知道行列式？大多數(shù)時候你是想知道行列式是否為零。但是,當(dāng)行列式為零、或者數(shù)值上格外接近于零時，由于前面提到的“舍入誤差”緣由很難進(jìn)行數(shù)值計算。這是奇異值分解的另一個實例.MATLＡB用det（A）命令計算奇異值.在MAＴＬAB命令行中手工輸入以下矩陣，或者從maｔｒix２。daｔ文件中復(fù)制以下矩陣：〉＞A＝［0。4５，—０.24４1１1,－0．０1９３３7３，0.32397２,-0。１１８829;－0．2４４11１，０．６８403６,－０．103427,０。2055６9,0。0029２382;-０.0１9３３73,—0。103427,0．8２９5，０.0１８9674，—0。０１1169；0.3２３972，０。20５56９，０.01８9674,0.６５9479,0。1９７３88；-０。１18829，０．０02923８2，-０.011１69，0.19738８,０．7769８5]用det命令可以計算出行列式:＞>dｅt(Ａ）ａns=—1。9682e—０08主軸理論：特征值和特征向量MATLＡB有預(yù)置的計算矩陣特征值和特征向量的函數(shù):＞〉eiｇ(Ａ)aｎs＝-0.００00０。70０00。8000０.9０001.０00０當(dāng)使用以下形式調(diào)用eig（)函數(shù)時，它會以矩陣V的列向量形式返回特征向量:>＞［V，Ｄ]=eiｇ（A)Ｖ＝—0．6836０。０00０-０。54６9－０。４7８5—0。０６８４—０.418１0。61６20.18310。4５3０—0.4５４７－０．0８３70.400３０。６189-0。6２320。2４７90.５４０90．2582－0。2４15－0.4042-０.647４—０.２416—0。62720．47５5—0.１1９０—０。5５50D＝－0。00０0０000０0.７0０0００0０００。８0000000０0．9000０0０0０1。0０00ｅｉｇs()函數(shù)是eig（）函數(shù)的變形，它能夠計算稀疏矩陣的特征值和特征向量。在下一節(jié)中將結(jié)合COMSOLMultｉphyｓｉcｓ介紹它的應(yīng)用。如果矩陣A的行列式格外接近于零，那么一個特征值也格外接近于零。特征向量是矩陣A的化零空間——給出映射到零的方向：＞〉Ａ＊Ｖ（:,1）ans＝1．0ｅ-0０7＊0。2669０．16３30.０32７－0．２1120.０９43全部特征向量可以通過乘以特征值等于矩陣Ａ的特性來證明這一點，例如:＞〉Ａ*Ｖ（：，2）.／V（:，2）aｎｓ=0．７００00。700０0．7０000。７0０00.7000在ＭＡＴLAＢ中,“./"符號表示元素對元素的除法，上面的冒號表示矩陣全部的列。由于系統(tǒng)接近奇異,全部包含該矩陣的解都很差，對這一點我們不應(yīng)該感到奇怪。例如:>>B=［0；１；0；1;０];>〉A(chǔ)\Banｓ=1．０e+０0６＊２．14８7１.３14２０.２6３1-1。70０１0.７593由于A的全部元素都是個位數(shù)，被除數(shù)B矢量元素也是個位數(shù)，方程（２7)的解也應(yīng)該由個位數(shù)組成,而不是百萬量級。對于化學(xué)工程師，這意味著對于一個質(zhì)量守恒系統(tǒng)，入口流速為1ｋg/ｈr，由于質(zhì)量守恒限制出口也為1ｋg／hr，但是反應(yīng)器中流速為100０00０kg／hｒ。這是不行能的。但這就是由近奇異矩陣給出的計算結(jié)果。奇異值分解（SVD)奇異值分解在很多情況下可以給出一個較好的解。類似于特征值和特征向量的主軸理論,全部矩陣都只有唯一的分解：? ??? ?(28）這里的U和V都是正交方陣,diaｇ是包含奇異值的對角矩陣.在解出Ｕ，V和ｄｉag的基礎(chǔ)上，(27）式可以輕易解得? ? ???(29)U和V為正交矩陣意味著其轉(zhuǎn)置矩陣也是其逆矩陣.對角陣的逆矩陣為對角元素取倒數(shù).所以求解的唯一問題就是何時一個或多個奇異值(ｄiａgj)接近于零.這時（1／diagj)是一個格外大的數(shù),導(dǎo)致求解結(jié)果錯誤.一個比較好的近似方法是將這些奇異值的(1／ｄｉagj）設(shè)為零!???? ? （３０）幾乎為零元素的替代向量應(yīng)該在數(shù)量級上最小,近似滿意方程。在我們矩陣A的例子中,如果只有一個輸出,ＭATLAB命令ｓvｄ()給格外異值，如果有三個輸出,命令sｖd(）給出:＞>［Ｕ,D,Ｖ]=svｄ（A）Ｕ＝－０。06８４-0。4785０。５46９０．０000-0。6836-0。4５４70．４５3０－０。1831－0．６162—０．４1８10.2４79－0.6232-０.6189—0．４0０３-0.0837－0。6４7４-0．4０４２0。2415—0．25820．540９－0.5550—0．１１90—０.４75５0。６2７2－0。241６D=1.00０0000000．９00000０0０0。８００000０00０.700０000０－０．00001。００0０V=－0。０６84-0．478５０．546９0.00０0－0。68３6—0．454７0．45３0—0.1831－0.6１6２-0。4１８10.２4７９-0．6２32—0.61８9－０.４003-0。0837—０.6474—0.4０４２０。2４15—0.25820.５409－０．5550—0。119０－0.47550。6２７２-０.2416由于Ａ為對稱矩陣,必定有Ｕ=V。最小數(shù)量級的奇異值分解命令如下：＞>ss=[1．1．／０。９1./0.８１./0．７０］;＞〉diｎｖ=diag（ss)；〉>V＊ｄinv*U’＊Ｂanｓ=０。０８９３１．2８２00。1４79１。０３1７-0.213０這是一個物理上更容易接受的答案,例如，前面商量的質(zhì)量守恒系統(tǒng)中的內(nèi)部氣體流速。由于有限元方法基于矩陣,所以線性系統(tǒng)理論對于COMSOＬＭultｉphｙsｉcs建模格外重要。當(dāng)產(chǎn)生的剛度矩陣變得接近奇異時,COMSOLＭｕｌtiｐhｙsｉcｓ可能不能給出滿意要求的結(jié)果.MATＬＡB中的矩陣計算和稀疏化方法可輕易地用于檢查COMSOLMｕltiphｙsics結(jié)果的合理性，也可以用于對系統(tǒng)自然動力學(xué)的特征分析。這些想法將在下一節(jié)中作為ＣOMSOＬＭｕｌtiphｙsｉcs的一個模型例子加以介紹。6.1非均勻介質(zhì)中的熱傳導(dǎo)常常遇到的穩(wěn)態(tài)熱傳導(dǎo)方程，是可以分析求解的最簡潔的偏微分方程:泊松方程。但是對于簡潔幾何情況,會比較難處理。作者認(rèn)為某些通解很難推導(dǎo)和收斂［5]。所以，數(shù)值解比通解更好,進(jìn)一步講，傳熱過程的任何轉(zhuǎn)變都可能會破壞分析結(jié)果.本節(jié)我們考慮非均勻厚板中的一維傳熱問題,在板中有分布熱源：? ?????（3１)打開耦合MATLAB/COＭＳOＬＳｃｒipt的COＭSOＬMuｌｔipｈysiｃs模型導(dǎo)航欄：依據(jù)表9中的步驟建立有分布熱源的介質(zhì)傳熱問題。依據(jù)計算結(jié)果可以畫出一條幾乎線性的曲線,斜率為－１。由于這是個線性問題,求解器日志顯示兩步迭代收斂。計算結(jié)果得T｜x＝0。5=０.4７4097,而解析解為Ｔ｜x=0.5～0。475：?????? （３2）下面嘗試ｋ=１－ｘ／２。這種情況下也有解析解,但是由于數(shù)值簡潔性,需要采納對數(shù)和分支切割。解析解為Ｔ|x＝０。5~0。550，你的計算結(jié)果呢?表9非均勻介質(zhì)中的熱傳導(dǎo)——分布熱源情況ModelNaviｇatｏr選擇1D空間維數(shù)，CＯMSOLMuｌｔiphｙsics:ＨｅaｔTransfｅｒ:Conｄｕctioｎ:Steadyｓtａt(yī)eanａlysiｓ,設(shè)定因變量：ｕ選擇Elｅｍeｎt:Ｌagｒanｇｅ—Linｅar，完成Ｄraw菜單Ｓpｅｃiｆｙobjｅcts：Linｅ,Ｃｏoｒdinaｔes彈出菜單,x:0１名稱：inｔeｒｖaｌ，完成Pｈysiｃs菜單：Boundarysettingｓ選擇邊界1設(shè)定邊界條件為:溫度；T0＝1選擇邊界2設(shè)定邊界：T1＝0,完成Phｙｓicｓ菜單：Sｕｂdomａｉnｓｅｔtings選擇域1設(shè)定ｋ=1;Ｑ=－ｘ＊(１—x）選擇Iｎit選項卡；設(shè)定Ｔ(ｔ0）=１-x,完成Ｍesｈinｇ點擊三角形按鈕,劃分網(wǎng)格Ｓolver點擊求解按鈕（＝)Post－ｐrocｅsｓｉngDａt(yī)ａdisｐlay設(shè)定x=０。５值:0．４74０97［K],表達(dá)式:Ｔ，位置：(０。5）下面來看線性系統(tǒng)理論。打開“Ｆile"菜單，選擇“ＥxｐoｒtFEＭstｒuctureaｓfem”。這是我們其次次導(dǎo)出為FEＭ結(jié)構(gòu)，有必要簡潔介紹一下。對于任何MATLAB/ＣOMSOLＳcript變量,在命令行中輸入變量名字,ＭATLＡB/CＯMSOＬScrｉpt將會給出變量值或其數(shù)據(jù)結(jié)構(gòu)。試一下：>〉ｆｅmfeｍ=version：［１x1ｓｔruｃt]ａpｐl:｛［１x１ｓｔrucｔ］｝geom：［1x1gｅoｍ1]meｓh:［１x１fｅｍmesｈ]bｏｒdｅr：１oｕｔｆorm：’general’form：’generaｌ’unitｓ:’SI’equ：[1x1struct]bnd:［1ｘ1strucｔ]dｒａｗ：[１x1struct]xmeｓｈ:[1x１cｏｍ.fｅmｌaｂ。xmesｈ.Xmｅsh］soｌ:[1x1ｆｅmsoｌ]當(dāng)然,你也可以進(jìn)一步查看FEM結(jié)構(gòu),了解整個樹的干、枝、葉.我們已經(jīng)剪除了fｅm.soｌ。u，fem。ｍesh。p1不同的分支包含了完整的COＭSOLMｕlｔiphｙｓicｓ模型-—為fem.gｅom中保存的幾何結(jié)構(gòu)定義了方程(fem.appl{1｝．equ）和邊界條件（fｅm．a(chǎn)ppl｛1}．eｑｕ）.與其它變量相同,這些結(jié)構(gòu)可以被ＭAＴLAＢ/CO