Matlab與化學(xué)化工計算ppt課件

.2010,計算機在化學(xué)化工中的應(yīng)用七Matlab與化學(xué)化工計算,.2010,本節(jié)要點,本章背景Matlab基礎(chǔ)方程組求解數(shù)據(jù)插值作業(yè),.2010,問題的提出,MATLAB語言與其它語言的關(guān)系仿佛和C語言與匯編語言的關(guān)系一樣計算機語言的發(fā)展標(biāo)志著計算機語言向“智能化”方向發(fā)展，被稱為第四代編程語言,.2010,1Matlab基礎(chǔ)知識,.2010,1.1Matlab簡介,1967年由ClereMaler用FORTRAN語言設(shè)計和編寫1984年Mathworks公司用C語言完成了Matlab的商業(yè)化版本并推向市場經(jīng)過20余年的改進，Matlab已發(fā)展成為一個具有極高通用性的、帶有眾多實用工具的運算平臺，成為國際上廣泛認(rèn)可的優(yōu)秀科學(xué)計算軟件,.2010,Matlab的發(fā)展,1984年，MATLAB第1版(DOS版)1992年，MATLAB4.0版1994年，MATLAB4.2版1997年，MATLAB5.0版1999年，MATLAB5.3版2000年，MATLAB6.0版2001年，MATLAB6.1版2002年，MATLAB6.5版2004年，MATLAB7.0版,告別DOS版,1993年MathWorks公司從加拿大滑鐵盧大學(xué)購得Maple的使用權(quán)，推出了符號計算工具包,5.0的MATLAB擁有更豐富的數(shù)據(jù)類型和結(jié)構(gòu)、更友善的面向?qū)ο?、更加快速精良的圖形可視、更廣博的數(shù)學(xué)和數(shù)據(jù)分析資源、更多的應(yīng)用開發(fā)工具,.2010,Matlab的優(yōu)點,語法簡單易學(xué)，編程效率高高質(zhì)量、高可靠的數(shù)值計算能力強大的矩陣運算能力高級圖形和數(shù)據(jù)可視化處理能力提供600多個常用算法內(nèi)建函數(shù)，以及眾多面向應(yīng)用的工具箱,.2010,Matlab二維作圖,.2010,Matlab三維作圖,.2010,1.2Matlab的界面,.2010,1.3Matlab的幫助功能,聯(lián)機幫助系統(tǒng)命令窗口查詢helplookfor聯(lián)機演示系統(tǒng)Demos,“Help”下拉菜單中“FullProductFamilyHelp”命令打開聯(lián)機幫助系統(tǒng),若不知函數(shù)確切名，可“Lookfor關(guān)鍵詞”可查,.2010,help,Help全部主題,Help指定函數(shù),.2010,例7-1,查找包含“diff”關(guān)鍵詞的函數(shù)lookfordiffSETDIFFSetdifference.DIFFDifferenceandapproximatederivative.POLYDERDifferentiatepolynomial.DDE23Solvedelaydifferentialequations(DDEs)withconstantdelays.DDESDSolvedelaydifferentialequations(DDEs)withgeneraldelays.DEVALEvaluatethesolutionofadifferentialequationproblem.,用戶輸入的命令,查詢結(jié)果,.2010,2線性方程組求解,.2010,2.1線性方程組的一般形式,在應(yīng)用中，常常把線性方程組寫成AX=b的一般形式，其中,.2010,2.2線性方程組解的判斷,齊次線性方程組AX=0，其解的情況可以通過系數(shù)矩陣A的秩和未知數(shù)個數(shù)n的關(guān)系來判斷如果系數(shù)矩陣的秩為n，方程組只有零解，x=0如果系數(shù)矩陣的秩小于n，方程組有無窮多解如果系數(shù)矩陣的秩大于n，方程組無解,.2010,非其次線性方程組解的情況,非齊次線性方程組AX=b，根據(jù)系數(shù)矩陣A的秩、增廣矩陣B=Ab的秩和未知數(shù)個數(shù)n的關(guān)系來判斷其解的情況如果系數(shù)矩陣A的秩等于增廣矩陣B的秩且等于n，方程組有唯一解如果系數(shù)矩陣A的秩等于增廣矩陣B的秩且小于n，方程組有無窮多解如果系數(shù)矩陣A的秩小于增廣矩陣B的秩，方程組無解,.,例7-2判斷方程解的情況,解：在Matlab中輸入a=-1-24;211;11-1;rank(a)ans=2齊次線性方程組系數(shù)矩陣A的秩為2，小于未知數(shù)個數(shù)3，方程組有無窮多解,計算系數(shù)矩陣A的秩,;不能少,.2010,例7-2（2）,解：a=7028;0281;280196;b=1-39-7;%b為列向量，故輸入行向量后轉(zhuǎn)置rank(a)%計算系數(shù)矩陣A的秩ans=3rank(ab)%計算增廣矩陣Ab的秩ans=3非齊次線性方程組系數(shù)矩陣A的秩為3，增廣矩陣的秩為3，等于未知數(shù)個數(shù)3，方程組有唯一解。,“%”是Matlab的注釋符，%后的語句作為注釋處理,.2010,2.3線性方程組直接求解,例7-3求以下方程組的解步驟b1矩陣除法，驗證解a判斷解的情況b2逆矩陣法b3rref,.2010,例7-4,求下列方程組的解,視頻演示,.2010,3數(shù)據(jù)插值,.2010,3.1數(shù)據(jù)插值簡介,在工程領(lǐng)域，許多實驗數(shù)據(jù)常以列表函數(shù)或表格的形式存在，如水黏度隨溫度的列表函數(shù)在實際使用時，有時需要獲得介于表中兩個溫度結(jié)點之間（如15，25）的黏度值。而這些數(shù)據(jù)未在表中出現(xiàn)，需要我們根據(jù)已知的數(shù)據(jù)估算出表中未出現(xiàn)的溫度點的黏度數(shù)值，這一技術(shù)稱為插值技術(shù),.2010,插值的數(shù)學(xué)定義,已知由g(X)(可能未知或非常復(fù)雜)產(chǎn)生的n+1個離散數(shù)據(jù)(xi,yi),i=0,1,2,n，且這n+1個互異插值結(jié)點滿足a=x0x1x2xn=b，在插值區(qū)間a,b內(nèi)尋找一個相對簡單的函數(shù)f(x)，使其滿足插值條件f(xi)=yi,i=0,1,2,n。再利用已求得的f(x)計算任一非插值結(jié)點x*處的近似值y*=f(x*)。其中f(x)稱為插值函數(shù)，g(x)稱為被插值函數(shù)從計算的觀點看，插值就是用一個簡單函數(shù)在某種誤差范圍內(nèi)近似的代替原目標(biāo)函數(shù)關(guān)系式,.2010,3.2插值方法,線性插值二次插值其他插值方法最近（nearest）插值法樣條曲線（spline）法埃爾米特（Hermite）法,.2010,3.2.1線性插值,又稱兩點插值已知兩個數(shù)據(jù)點x0,y0,x1,y1(x0x1)，求對應(yīng)于x(x0xx1)的y值解法：由x0,y0,x1,y1構(gòu)造直線方程并求取在該點的函數(shù)值,線性插值的優(yōu)點是簡單，快捷，特別是對于插值結(jié)點間距較小的情況可以取得令人滿意的精度,.2010,3.2.2二次插值,又稱拉格朗日三點差值根據(jù)三個已知點x0,y0,x1,y1,x2,y2(x0x1x2)，構(gòu)造二次多項式插值函數(shù)y=a0+a1x+a2x2，并用該函數(shù)計算在x處的y值二次插值公式,.2010,3.3.1使用Matlab進行數(shù)據(jù)插值,一維插值只有一個自變量的插值Matlab提供的一維插值函數(shù)是interp1常用語法：YI=interp1(X,Y,XI,method)式中X,Y為已知數(shù)據(jù)點的x,y值；XI為待插值數(shù)據(jù)點的x值；YI為返回的插值結(jié)果；method用于指定所采用的插值方法,.2010,插值函數(shù)interp1提供的插值方法,.2010,不同方法插值的結(jié)果,在0,2區(qū)間內(nèi)生成11個等距的離散點，計算函數(shù)y=sin(x)的數(shù)值,分段三次Hermite插值,分段三次樣條插值,線性插值,最近插值,.2010,例7-5,用函數(shù)y=ex生成以下離散數(shù)據(jù)，使用Matlab的不同插值方法計算x=2.552.632.772.86處的函數(shù)值，并與真實值進行比較,視頻演示,.2010,例7-5計算結(jié)果的比較,最接近真實值,.2010,例7-6,已知水在20,21,22,23的飽和蒸汽壓分別為17.54,18.65,19.83,21.07mmHg,求20.5，21.5，22.5和24時水的飽和蒸汽壓各是多少？已知24時水的飽和蒸汽壓為22.38mmHg,視頻演示,.2010,3.3.2多維插值,具有多個自變量的插值二維插值三維插值高維插值,插值函數(shù),Method選項,.2010,例7-7,函數(shù)z=ex+sin(y)+y-1生成表7-6中的離散數(shù)據(jù)，應(yīng)用不同插值方法計算在x=0.36，y=1.9處的z值，并與真值作比較,.2010,例7-7插值結(jié)果,視頻演示,最接近真實值,.2010,4非線性方程的求解,.2010,4.1非線性方程數(shù)值求解,常見的非線性方程（組）有兩種形式或,直接迭代法、韋斯頓迭代法求解,牛頓迭代法求解,.2010,4.1.1直接迭代法,先設(shè)定x的初值x=x0，代入方程計算x1=f(x0)，再把x1作為新的初值代入方程計算x2=f(x1)直至求得的xn+1與xn足夠接近（稱為收斂），xn即為方程的根直接迭代法求解可寫為如下形式,直接迭代法的優(yōu)點是形式簡單，易于編程實現(xiàn)。缺點是計算量大、收斂速度慢。一般可通過改進初值、降低收斂要求等方法提高其收斂速度。也可采用其它方法進行求解,.2010,收斂判斷準(zhǔn)則,絕對偏差相對偏差半相對偏差,是用戶指定的一個很小的正數(shù)，確定適當(dāng)?shù)娜≈涤幸欢y度,優(yōu)點在于的選取不受方程根的數(shù)值大小的影響。一般取=0.001,是判斷收斂的一個較好的方法，的取值一般為0.00010.001,.2010,4.1.2韋格斯頓迭代法,韋斯頓法對直接迭代做了改進，使用前兩個計算點的信息進行求解。其迭代形式為韋格斯頓法迭代時需要前面兩個計算點的數(shù)據(jù)，可先執(zhí)行一次直接迭代法計算獲得,.2010,4.1.3牛頓迭代法,又稱切線法若將f(x)在其根附近進行泰勒級數(shù)展開，并取級數(shù)的線性部分作為f(x)的近似值，可得由上式可得牛頓迭代公式如函數(shù)的一階導(dǎo)數(shù)難以求得，可用差商作為近似導(dǎo)數(shù),.2010,4.2用Matlab求解非線性方程,Matlab求解非線性方程（組）的函數(shù)fzero：一元非線性方程的求解fsolve：用于非線性方程組的求解fzero函數(shù)用法x=fzero(fun,x0)fsolve函數(shù)用法x=fsolve(fun,x0),fun單變量實值函數(shù)，可以是Matlab內(nèi)部函數(shù)或用戶自定義函數(shù)x0若x0是一個單個的數(shù)值，系統(tǒng)會將其作為求解的初值，在其附近尋找解；若x0是一個二維向量，且fun(x0(1)和fun(x0(2)符號相反，Matlab將會在x0(1)和x0(2)區(qū)間內(nèi)尋找零點,fun用戶自定義函數(shù)，返回給定變量x時方程(組)的值y=fun(x)x0初值矩陣,對于fzero和fsolve函數(shù)，給定適當(dāng)?shù)某踔祵栴}的求解至關(guān)重要，若初值選擇不當(dāng)，將無法得到正確的解。一般可根據(jù)經(jīng)驗或簡化計算獲得合適的初值,.2010,例7-8,試用維里方程計算200，1.013MPa的異丙醇蒸汽的摩爾體積V與壓縮因子Z。已知異丙醇的維里系數(shù)實驗值B=-388cm3mol-1，C=-26000cm6mol-2,視頻演示,.2010,例7-9,600K下由CH3Cl和H2O反應(yīng)生成CH3OH，存在下列平衡CH3Cl(g)+H2O(g)=CH3OH(g)+HCl(g)(1)2CH3OH(g)=(CH3)2O(g)+H2O(g)(2)已知該溫度下Kp(1)=0.00154,Kp(2)=10.6。今以等摩爾的CH3Cl(g)和H2O(g)開始反應(yīng)，求CH3Cl的平衡轉(zhuǎn)化率,視頻演示,.2010,5常微分方程（組）求解,.2010,5.1化工中的常微分方程(組),微分方程中只有一個自變量的方程稱為常微分方程，自變量個數(shù)為兩個或兩個以上的微分方程稱為偏微分方程常微分方程初值問題：給定微分方程及初值條件邊值問題：給定微分方程及邊界條件常微分方程（組）的解法解析法和數(shù)值法（常用）,.2010,初值問題,記為或,.2010,5.2常微分方程（組）數(shù)值解法,歐拉公式梯形公式龍格-庫塔法常微分方程組的數(shù)值解法,.2010,5.2.1歐拉公式,若常微分方程初值問題的求解區(qū)間為，將其等分為m步，步長。記，相應(yīng)xn處的函數(shù)值為yn,，則yn可由下式計算向前歐拉公式向后歐拉公式中心歐拉公式,若yn+1同時出現(xiàn)在等號的兩側(cè)，稱為隱式歐拉公式，無法直接求解，一般需采用迭代法計算,.2010,5.2.2梯形公式,梯形公式也是隱式格式，需要迭代求解先用顯式公式算出初值，再用隱式公式進行一次或數(shù)次修正。這一過程稱為預(yù)估-校正過程公式為合并為,.2010,5.2.3龍格-庫塔法,工程應(yīng)用中求解常微分方程最常用的一種有效方法，其計算精度和運算速度較快，易于編程。常用的有二階、三階、四階龍格-庫塔公式二階龍格-庫塔公式，常見形式或,.2010,三階龍格-庫塔公式,常見形式,.2010,四階龍格-庫塔公式,常見形式,.2010,5.2.4常微分方程組的數(shù)值解法,將由m個一階方程組成的常微分初值問題其中,可由前邊所述的解常微分方程的各個方法求解,寫為向量形式,.2010,5.2.5高階常微分方程數(shù)值解法,可把高階常微分方程轉(zhuǎn)化為一階常微分方程組求解。例如三階常微分方程令將三階方程化為一階方程,.2010,5.3用Matlab求解常微分方程,.2010,例7-10,在間歇反應(yīng)器中進行液相反應(yīng)制備產(chǎn)物B，其反應(yīng)網(wǎng)絡(luò)如圖7-7所示。反應(yīng)溫度為224.6，反應(yīng)物X大量過剩。各反應(yīng)均為一級動力學(xué)關(guān)系：,各步反應(yīng)的k0i、Eai見表，試給出0-10000秒各產(chǎn)物的濃度變化規(guī)律。初始條件為：t=0，CA=1kmol/m3，CB=CC=CD=CE=0kmol/m3,反應(yīng)網(wǎng)絡(luò)圖,.2010,例7-10條件圖,反應(yīng)網(wǎng)絡(luò)圖,參數(shù)取值,.2010,例7-10分析步驟,分析步