1、第8章 Matlab數(shù)值微分和積分化工學(xué)院 伍水生2019-11第8章 Matlab數(shù)值微分和積分化工學(xué)院 伍水生本章知識(shí)要點(diǎn)極限 (limit)數(shù)值微分(polyder, fnder)數(shù)值積分(quad, quadl, fnint)本章知識(shí)要點(diǎn)數(shù)值微分、數(shù)值積分在化工計(jì)算中的作用狀態(tài)方程計(jì)算流體的焓和熵微分法反應(yīng)動(dòng)力學(xué)方程擬合等溫活塞流反應(yīng)器的設(shè)計(jì)計(jì)算微觀離析反應(yīng)器的計(jì)算數(shù)值微分、數(shù)值積分在化工計(jì)算中的作用狀態(tài)方程計(jì)算流體的焓和單變量函數(shù)的極限格式1： L= limit( fun, x, x0)格式2： L= limit( fun, x, x0, left 或 right)8.1 數(shù)值微分8

2、.1.1 極限問題的解析解單變量函數(shù)的極限8.1 數(shù)值微分8.1.1 極限問題的解析例: 試求解極限問題 syms x a b; L=limit(x*(1+a/x)x*sin(b/x),x,inf) L = exp(a)*b例:求解單邊極限問題 syms x; limit(exp(x3)-1)/(1-cos(sqrt(x-sin(x),x,0,right) ans =12例: 試求解極限問題多變量函數(shù)的極限：格式： L1=limit(limit(f,x,x0),y,y0) 或 L1=limit(limit(f,y,y0), x,x0) 如果x0 或y0不是確定的值，而是另一個(gè)變量的函數(shù)，如x-

3、g(y)，則上述的極限求取順序不能交換。多變量函數(shù)的極限：例：求出二元函數(shù)極限值 syms x y a; f=exp(-1/(y2+x2) *sin(x)2/x2*(1+1/y2)(x+a2*y2); L=limit(limit(f,x,1/sqrt(y),y,inf)L =exp(a2)例：求出二元函數(shù)極限值8.2 函數(shù)導(dǎo)數(shù)的解析解函數(shù)的導(dǎo)數(shù)和高階導(dǎo)數(shù)格式： y=diff(fun,x) %求導(dǎo)數(shù)(默認(rèn)為1階) y= diff(fun,x,n) %求n階導(dǎo)數(shù)例： 求一階導(dǎo)數(shù)： syms x; f1=diff(sin(x)/(x2+4*x+3),x)8.2 函數(shù)導(dǎo)數(shù)的解析解函數(shù)的導(dǎo)數(shù)和高階導(dǎo)數(shù)s

4、yms x;f1=diff(sin(x)/x)得結(jié)果: f1=cos(x)/x-sin(x)/x2解: 輸入命令: syms x; f1=diff(log(sin(x)得結(jié)果: f1=cos(x)/sin(x).在matlab中，函數(shù)lnx用log(x)表示，而log10(x)表示 lgx 。syms x;解: 輸入命令: 在matlab中，函數(shù)l多元函數(shù)的偏導(dǎo)：格式： f=diff(diff(f,x,m),y,n) 或 f=diff(diff(f,y,n),x,m)例： 求偏導(dǎo)數(shù) syms x y; z=(x2-2*x)*exp(-x2-y2-x*y); zx=diff(z,x)zx =(2

5、*x-2)*exp(-x2-y2-x*y)+(x2-2*x)*(-2*x-y)*exp(-x2-y2-x*y)多元函數(shù)的偏導(dǎo)： zy=diff(z,y)zy =(x2-2*x)*(-2*y-x)*exp(-x2-y2-x*y)直接繪制三維曲面 x,y=meshgrid(-3:.2:3,-2:.2:2); z=(x.2-2*x).*exp(-x.2-y.2-x.*y); surf(x,y,z), axis(-3 3 -2 2 -0.7 1.5) 第8章_matlab數(shù)值微分與積分課件隱函數(shù)的偏導(dǎo)數(shù)：格式：F=diff(f,xj)/diff(f,xi)隱函數(shù)的偏導(dǎo)數(shù)：例： z=f(x,y)=-x2

6、-y2-x*y; syms x y; f=x2-y2-x*y; diff(f,x)/diff(f,y)或者： syms x y; diff(x2-y2-x*y,x)/diff(x2-y2-x*y,y)ans =(-2*x-y)/(-2*y-x) 例： z=f(x,y)=-x2-y2-x*y;對(duì)于列表型函數(shù)往往需要用數(shù)值方法計(jì)算函數(shù)的微分?jǐn)?shù)值微分的基本方法差分利用插值（擬合）多項(xiàng)式求微分利用三次樣條插值（擬合）函數(shù)求微分?jǐn)?shù)值微分可以放大誤差，應(yīng)謹(jǐn)慎使用對(duì)于列表型函數(shù)往往需要用數(shù)值方法計(jì)算函數(shù)的微分函數(shù)diff對(duì)于向量X，diff(X)表示了X(2)-X(1) X(3)-X(2) . X(n)-X

7、(n-1).對(duì)于矩陣X，diff(X)表示了X(2:n,:) - X(1:n-1,:) diff(x,n,dim)得到矩陣x在dim維上的n階差值 diff(1:10).2)1 3 5 7 9 11 13 15 17 19x=1 3 8; 2 4 6diff(x,1,1)diff(x,1,2)diff(x,2,2)diff(x,3,2)1 1 -22 5;2 2Empty matrix: 2-by-030函數(shù)diff對(duì)于向量X，diff(X)表示了X(2)-X(Matlab數(shù)值微分實(shí)現(xiàn)方法有限差分法：用差分函數(shù)diff()近似計(jì)算導(dǎo)數(shù)多項(xiàng)式擬合方法：離散數(shù)據(jù) 多項(xiàng)式擬合函數(shù) 導(dǎo)函數(shù)pp 在xi

8、的導(dǎo)數(shù)值 三次樣條插值方法樣條擬合方法離散數(shù)據(jù) 樣條插值函數(shù)cs導(dǎo)函數(shù)pp 在xi的導(dǎo)數(shù)值 離散數(shù)據(jù) 樣條擬合函數(shù)sp導(dǎo)函數(shù)pp 在xi的導(dǎo)數(shù)值 Matlab數(shù)值微分實(shí)現(xiàn)方法有限差分法：用差分函數(shù)diff(例題某一液相反應(yīng)濃度隨時(shí)間變化的實(shí)驗(yàn)數(shù)據(jù)如下表：t/min00.20.61.02.05.010.0C(g/L)5.193.772.301.570.80.250.094試t=0.1,0.4min時(shí)的反應(yīng)速度 x=0 0.2 0.6 1 2 5 10; C=5.19 3.77 2.30 1.57 0.8 0.25 0.094; f=spline(x,C); dC=fnder(f); dC1=fn

9、val(dC,0.1); dC2=fnval(dC,0.4); dC1, dC2三次樣條插值方法:例題某一液相反應(yīng)濃度隨時(shí)間變化的實(shí)驗(yàn)數(shù)據(jù)如下表：t/min0 x=0 0.2 0.6 1 2 5 10; C=5.19 3.77 2.30 1.57 0.8 0.25 0.094; f=polyfit(x,C,3); dC=polyder(f); dC1=polyval(dC,0.1); dC2=polyval(dC,0.4);dC1, dC2多項(xiàng)式擬合方法： x=0 0.2 0.6 1 2 5 10;多項(xiàng)式擬8.2 數(shù)值積分?jǐn)?shù)值積分在數(shù)值計(jì)算中有著重要作用,許多數(shù)值計(jì)算問題可以轉(zhuǎn)化為數(shù)值積分問題

10、,如常微分方程初值問題等通?？捎帽平囗?xiàng)式Pn(x)來代替被積函數(shù)f(x)，計(jì)算積分構(gòu)造數(shù)值積分的方法很多，主要有Newton-Cotes系列數(shù)值積分法、Gauss積分法和Romberg積分法等8.2 數(shù)值積分?jǐn)?shù)值積分在數(shù)值計(jì)算中有著重要作用,許多數(shù)值計(jì)數(shù)值積分MATLAB函數(shù)公式quad自適應(yīng)Simpson求積公式（低階）quadl自適應(yīng)Lobatto求積公式；精度高，最常用trapz梯形求積公式；速度快，精度差cumtrapz梯形法求一個(gè)區(qū)間上的積分曲線cumsum等寬距法求一個(gè)區(qū)間上的積分曲線，精度很差fnint利用樣條函數(shù)求不定積分；與spline，ppval配合使用，主要應(yīng)用于表格“

11、函數(shù)”積分?jǐn)?shù)值積分MATLAB函數(shù)公式quad自適應(yīng)Simpson求積梯形法數(shù)值積分：trapz()調(diào)用格式：z=trapz(x,y) 用梯形求積方法計(jì)算y的積分近似值。對(duì)于向量y，trapz(y)返回y的積分；對(duì)于矩陣y，trapz(y)返回一行向量，向量中的各元素為矩陣y的對(duì)應(yīng)列 向量的積分值；梯形法數(shù)值積分：trapz()調(diào)用格式：z=trapz(x格式： S=trapz(x,y)例： x1=0:pi/30:pi; y=sin(x1) cos(x1) sin(x1/2); S1=trapz(x1,y) S1 = 1.9982 0.0000 1.9995第8章_matlab數(shù)值微分與積分課

12、件等同于： x=0:pi/30:pi; y1=sin(x); y2=cos(x) ; y3=sin(x1/2); S1=trapz(x,y1) ; S2=trapz(x,y2) ; S3=trapz(x,y3) ;S1, S2, S3等同于：調(diào)用格式：q=quad(fun,a,b)q=quad(fun,a,b,tol)輸入?yún)?shù)：fun被積函數(shù)。在定義fun時(shí)，被積函數(shù)表達(dá)式必須是向量形式，即表達(dá)式必須使用 點(diǎn)運(yùn)算符(.*、./和.)以支持向量a,b即積分限a,btol 絕對(duì)誤差限，默認(rèn)值為1.e-6輸出參數(shù)：q 積分結(jié)果自適應(yīng)Lobatto法數(shù)值積分：quadl() 調(diào)用格式同quad自適應(yīng)S

13、impson法數(shù)值積分：quad()調(diào)用格式：q=quad(fun,a,b)自適應(yīng)Simp例：第一種，一般定義函數(shù)方法 y=quad(c3ffun,0,1.5)y = 0.9661第一種，一般定義函數(shù)方法 y=quad(c3ffun用 inline 函數(shù)定義被積函數(shù)： f=inline(2/sqrt(pi)*exp(-x.2),x); y=quad(f,0,1.5)y = 0.9661y=quad(2/sqrt(pi)*exp(-x.2), 0,1.5)y = 0.9661第8章_matlab數(shù)值微分與積分課件3.4.5 雙重積分問題的數(shù)值解矩形區(qū)域上的二重積分的數(shù)值計(jì)算格式： 矩形區(qū)域的雙重

14、積分： y=dblquad(Fun,xm,xM,ym,yM) 限定精度的雙重積分： y=dblquad(Fun,xm,xM,ym,yM， )3.4.5 雙重積分問題的數(shù)值解矩形區(qū)域上的二重積分的數(shù)值計(jì)例：求解 f=inline(exp(-x.2/2).*sin(x.2+y),x,y); y=dblquad(f,-2,2,-1,1)y = 1.57449318974494或者：function z=fun(x,y)z=exp(-x.2/2).*sin(x.2+y);y=dblquad(fun,-2,2,-1,1)或者： y=dblquad(exp(-x.2/2).*sin(x.2+y),-2,2

15、,-1,1)例：求解3.4.6 三重定積分的數(shù)值求解格式:I=triplequad(Fun,xm,xM,ym,yM,zm,zM,tol ,quadl) 其中quadl為具體求解一元積分的數(shù)值函數(shù),也可選用quad或自編積分函數(shù),但調(diào)用格式要與quadl一致。3.4.6 三重定積分的數(shù)值求解例： triplequad(inline(4*x.*z.*exp(-x.*x.*y-z.*z), x,y,z), 0, 1, 0, pi, 0, pi,1e-7,quadl)ans = 1.7328或者： triplequad(4*x.*z.*exp(-x.*x.*y-z.*z), 0, 1, 0, pi,

16、0, pi,1e-7,quadl)ans = 1.7328例：基于樣條插值的數(shù)值積分運(yùn)算格式： f=fnint(S)其中S為樣條函數(shù)。例：考慮 中較稀疏的樣本點(diǎn),用樣條積分的方式求出定積分及積分函數(shù)。 x=0,0.4,1 2,pi; y=sin(x); sp1=csapi(x,y); a=fnint(sp1,1); 建立三次樣條函數(shù) xx=fnval(a,0,pi); xx(2)-xx(1)ans = 2.0191基于樣條插值的數(shù)值積分運(yùn)算反應(yīng)器停留時(shí)間分布的混合特性在t0的時(shí)刻，在一容器入口處突然向流進(jìn)容器的流體脈沖注入一定量的示蹤劑，同時(shí)在容器出口處測(cè)量流出物料中示蹤劑濃度隨時(shí)間的變化，實(shí)

17、驗(yàn)數(shù)據(jù)如下表：t/s020406080100120140160180200C(kmol/m3)10300000.45.516.211.11.70.10試計(jì)算流體在容器中的平均停留時(shí)間以及擴(kuò)散準(zhǔn)數(shù)。數(shù)學(xué)模型：平均停留時(shí)間方差： 擴(kuò)散特征數(shù)為： 反應(yīng)器停留時(shí)間分布的混合特性在t0的時(shí)刻，在一容器入口處突t = 0:20:200;C = 0, 0, 0, 0, 0.4, 5.5, 16.2, 11.1, 1.7, 0.1, 0; sp1= spline(t,C); %樣條函數(shù)擬合 sp= spaps(t,C,1); %導(dǎo)函數(shù) t0 = 0;tf = t(end); IC = quadl(Func,t

18、0,tf,sp); ICt = quadl(Func1,t0,tf,sp); ICt2 = quadl(Func2,t0,tf,sp); tm1 = ICt/IC; ss1 = ICt2/IC - tm12; DL2uL1 = ss1/tm12/2; tm1 , ss1 , DL2uL1function y = Func(x,sp) % f= Cy = fnval(sp,x); function y = Func1(x,sp) % f= Cty = fnval(sp,x).*x;function y = Func2(x,sp) % f= Ct2y = fnval(sp,x).*x.2;t = 0:20:200;C = 0, 0, 0,微分法進(jìn)行動(dòng)力學(xué)數(shù)據(jù)分析反應(yīng)物A在一等溫間歇反應(yīng)器中發(fā)生反應(yīng)為： ，測(cè)量得到反應(yīng)器中不同時(shí)間下反應(yīng)物A的濃度CA如下表所示。試根據(jù)表中數(shù)據(jù)確定其反應(yīng)速率方程。數(shù)學(xué)模型：t/s0204060120180300CA(mol/L)10865321微分法進(jìn)行動(dòng)力學(xué)數(shù)據(jù)分析反應(yīng)物A在一等溫間歇反應(yīng)器中發(fā)