數(shù)值微分與數(shù)值積分

第5章

數(shù)值微分與數(shù)值積分1/13/20231微積分是高等數(shù)學(xué)中的重要內(nèi)容，在化學(xué)工程上有許多非常重要的應(yīng)用微積分的數(shù)值方法，不同于高等數(shù)學(xué)中的解析方法，尤其適合求解沒(méi)有或很難求出微分或積分表達(dá)式的實(shí)際化工問(wèn)題的計(jì)算，例如：列表函數(shù)求微分或積分引言1---數(shù)值微積分方法不同于解析方法1/13/20232數(shù)值微分和數(shù)值積分與插值和擬合往往是密不可分的如在進(jìn)行數(shù)值微分時(shí)，常針對(duì)離散的數(shù)據(jù)點(diǎn)，利用插值和擬合可以減少誤差；而數(shù)值積分的基本思路也來(lái)自于插值法。例如如果所積函數(shù)的形式比較復(fù)雜或以表格形式給出，則可通過(guò)構(gòu)造一個(gè)插值多項(xiàng)式來(lái)代替原函數(shù)，從而使問(wèn)題簡(jiǎn)化引言2---數(shù)值微分和數(shù)值積分與插值和擬合關(guān)系密切1/13/20233§5.1數(shù)值微分

化工領(lǐng)域的實(shí)際問(wèn)題中時(shí)常需要求列表函數(shù)在節(jié)點(diǎn)和非節(jié)點(diǎn)處的導(dǎo)數(shù)值，這是數(shù)值微分所要解決的問(wèn)題。數(shù)值微分方法可近似求出某點(diǎn)的導(dǎo)數(shù)值例如在反應(yīng)動(dòng)力學(xué)的研究中，根據(jù)實(shí)驗(yàn)數(shù)據(jù)確定反應(yīng)的動(dòng)力學(xué)方程：這里實(shí)驗(yàn)測(cè)得一批離散點(diǎn)，要計(jì)算只能借助數(shù)值微分求導(dǎo)解決表5-1反應(yīng)動(dòng)力學(xué)實(shí)驗(yàn)數(shù)據(jù)

0導(dǎo)入1/13/202340.1建立數(shù)值微分公式的三種思路

常用三種思路建立數(shù)值微分公式：從微分定義出發(fā)，通過(guò)近似處理，得到數(shù)值微分的近似公式從插值近似公式出發(fā)，對(duì)插值公式的近似求導(dǎo)可得到數(shù)值微分的近似公式先用最小二乘擬合方法根據(jù)已知數(shù)據(jù)得近似函數(shù)，再對(duì)此近似函數(shù)求微分可得到數(shù)值微分的近似公式。然后對(duì)各方法數(shù)值微分后得到的多項(xiàng)式求值，即可求出任意點(diǎn)處的任意階微分1/13/202351方法概述在微積分中，一階微分的計(jì)算可以在二相鄰點(diǎn)x+h和x間函數(shù)取下列極限求得：

取其達(dá)到極限前的形式，就得到以下微分的差分近似式：

注：高階微分項(xiàng)可以利用低階微分項(xiàng)來(lái)計(jì)算，如二階微分式可以表示為：對(duì)應(yīng)的差分式有：§5.1.1差分近似微分上式中三種不同表示形式依次是一階前向差分、一階后向差分和一階中心差分來(lái)近似表示微分。其中一階中心差分的精度較高。1/13/202362差分的Matlab實(shí)現(xiàn)在Matlab中，可用diff函數(shù)進(jìn)行離散數(shù)據(jù)的近似求導(dǎo)調(diào)用形式：Y=diff(X,n)其中：X表示求導(dǎo)變量，可以是向量或矩陣。如是矩陣形式則按各列作差分；

n表示n階差分，即差分n次

Y是X的差分結(jié)果注：用diff函數(shù)進(jìn)行離散數(shù)據(jù)的近似求導(dǎo)與前向差分近似，但誤差較大。最好將數(shù)據(jù)利用插值或擬合得到多項(xiàng)式，然后對(duì)近似多項(xiàng)式進(jìn)行微分1/13/20237例5.1：丁二烯的氣相二聚反應(yīng)如下：實(shí)驗(yàn)在一定容器的反應(yīng)器中進(jìn)行，3260C時(shí)，測(cè)得物系中丁二烯的分壓(mmHg)與時(shí)間的關(guān)系如表5-2所示。。表5-2丁二烯二聚反應(yīng)實(shí)驗(yàn)數(shù)據(jù)用數(shù)值微分法計(jì)算所列時(shí)刻每一瞬間的反應(yīng)速率

1/13/20238解：程序如下：t=[0:5:90];pA=[632.0590.0552.0515.0485.0458.0435.0414.0396.0378.0362.0348.0336.0325.0314.0304.0294.0284.0274.0];dt=diff(t);％求時(shí)間t的差分dpA=diff(pA);％求壓力的差分q=dpA./dt％q為數(shù)值微分結(jié)果執(zhí)行結(jié)果：q=Columns1through8-8.4000-7.6000-7.4000-6.0000-5.4000-4.6000-4.2000-3.6000Columns9through16-3.6000-3.2000-2.8000-2.4000-2.2000-2.2000-2.0000-2.0000Columns17through18-2.0000-2.00001/13/20239§5.1.2三次樣條插值函數(shù)求微分若三次樣條插值函數(shù)收斂于，那么導(dǎo)數(shù)收斂于，因此用樣條插值函數(shù)作為函數(shù)，不但彼此的函數(shù)值非常接近，而且導(dǎo)數(shù)值也很接近。

用三次樣條插值函數(shù)求數(shù)值導(dǎo)數(shù)是可靠的，這是化工計(jì)算中求數(shù)值微分的有效方法。的近似1/13/2023101方法概述用三次樣條插值函數(shù)建立的數(shù)值微分公式為：

求導(dǎo)得：

（5-2）;式（5-1）和（5-2）不但適用于求節(jié)點(diǎn)處的導(dǎo)數(shù)，而且可求非節(jié)點(diǎn)處的導(dǎo)數(shù)。（5-1）其中，1/13/2023112三次樣條插值函數(shù)求微分的Matlab函數(shù)Matlab求離散數(shù)據(jù)的三次樣條插值函數(shù)微分方法分三個(gè)步驟：Step1:對(duì)離散數(shù)據(jù)用csapi函數(shù)（或spline函數(shù)）得到其三次樣條插值函數(shù)調(diào)用形式pp=csapi(x,y)其中：x,y分別為離散數(shù)據(jù)對(duì)的自變量和因變量；

pp為得到的三次樣條插值函數(shù)。Step2:用fnder函數(shù)求三次樣條插值函數(shù)的導(dǎo)數(shù)調(diào)用形式fprime=fnder(f,dorder)其中：f為三次樣條插值函數(shù)；

dorder為三次樣條插值函數(shù)的求導(dǎo)階數(shù)；

fprime為得到的三次樣條插值函數(shù)的導(dǎo)函數(shù)。Step3:可用fnval函數(shù)求導(dǎo)函數(shù)在未知點(diǎn)處的導(dǎo)數(shù)值調(diào)用形式

v=fnval(fprime,x)其中：fprime為三次樣條插值函數(shù)導(dǎo)函數(shù)；

x為未知點(diǎn)處自變量值；

v為未知點(diǎn)處的導(dǎo)數(shù)值。1/13/202312例5.2：某液體冷卻時(shí)，溫度隨時(shí)間的變化數(shù)據(jù)如表5-3所示：

表5-3冷卻溫度隨時(shí)間的變化數(shù)據(jù)試分別計(jì)算t＝2，3，4min及t＝1.5，2.5，4.5min時(shí)的降溫速率。解：三次樣條插值函數(shù)求數(shù)值微分的程序如下：t=[0:5];T=[92,85.3,79.5,74.5,70.2,67];cs=csapi(t,T);%生成三次樣條插值函數(shù)pp=fnder(cs);%生成三次樣條插值函數(shù)的導(dǎo)函數(shù)t1=[2,3,4,1.5,2.5,4.5];dT=fnval(pp,t1);%計(jì)算導(dǎo)函數(shù)在t1處的導(dǎo)數(shù)值disp('相應(yīng)時(shí)間時(shí)的降溫速率:')disp([t1;dT])執(zhí)行結(jié)果：相應(yīng)時(shí)間時(shí)的降溫速率:2.00003.00004.00001.50002.50004.5000-5.3722-4.6722-3.8389-5.7972-4.9889-3.2222注：前者是計(jì)算節(jié)點(diǎn)處的一階導(dǎo)數(shù)，后者是計(jì)算非節(jié)點(diǎn)處的一階導(dǎo)數(shù)1/13/202313§5.1.3最小二乘法樣條擬合函數(shù)求微分

在實(shí)際化工應(yīng)用中，當(dāng)來(lái)自實(shí)驗(yàn)觀測(cè)的離散數(shù)據(jù)不可避免地含有較大隨機(jī)誤差時(shí)，此時(shí)用插值公式求數(shù)值微分雖然樣本點(diǎn)處誤差較小，但可能會(huì)使非樣本點(diǎn)處產(chǎn)生較大誤差為此，可采用最小二乘法樣條擬合實(shí)驗(yàn)數(shù)據(jù)，獲得一個(gè)函數(shù)模型，然后再對(duì)其求導(dǎo)數(shù)。與樣條插值不同，樣條擬合不要求曲線經(jīng)過(guò)全部的數(shù)據(jù)點(diǎn)，這樣處理求導(dǎo)結(jié)果會(huì)有很大改善1/13/2023141方法概述

可用最小二乘法擬合成平滑B樣條曲線，即對(duì)于離散數(shù)據(jù)（所求的K次樣條擬合函數(shù)滿足：其中：再對(duì)擬合函數(shù)作平滑處理后求導(dǎo)，即可求出任意點(diǎn)處微分。）為權(quán)重系數(shù)，默認(rèn)為11/13/202315Matlab求離散數(shù)據(jù)的最小二乘法平滑B樣條擬合函數(shù)求微分共三個(gè)步驟：Step1:對(duì)離散數(shù)據(jù)用spaps函數(shù)得到最小二乘平滑B樣條擬合函數(shù)。調(diào)用格式：sp=spaps(x,y,tol)其中：

x,y――要處理的離散數(shù)據(jù)（）

tol――光滑時(shí)的允許精度，通常?。?0-2－10-4）Step2:可用fnder函數(shù)求最小二乘平滑B樣條擬合函數(shù)的導(dǎo)數(shù);Step3:可用fnval函數(shù)求導(dǎo)函數(shù)在未知點(diǎn)處的導(dǎo)數(shù)值。

fnder()和fnval（）調(diào)用形式前文中已經(jīng)介紹過(guò)。2最小二乘法平滑B樣條擬合函數(shù)求微分的Matlab實(shí)現(xiàn)1/13/202316由離散數(shù)據(jù)求數(shù)值微分的四種方法及有關(guān)MATLAB函數(shù)：（1）差分法用差分函數(shù)diff()近似計(jì)算導(dǎo)數(shù)，即dy=diff(y)./diff(x)。對(duì)于向量X，diff(X)表示了[X(2)-X(1)X(3)-X(2)...X(n)-X(n-1)].對(duì)于矩陣X，diff(X)表示了[X(2:n,:)-X(1:n-1,:)]小結(jié)注：此法用一階差分，精度較差，若改用二階差分，可大大提高精度，但編程麻煩些。（2）多項(xiàng)式擬合方法向量p表示的多項(xiàng)式擬合函數(shù)導(dǎo)函數(shù)pppp在xi的導(dǎo)數(shù)值。其中：函數(shù)polyfit()和polyval()在前文中已介紹導(dǎo)函數(shù)polyder()的調(diào)用格式為：pp=polyder(p)

離散數(shù)據(jù)該函數(shù)對(duì)向量p表示的多項(xiàng)式函數(shù)進(jìn)行求導(dǎo)，返回導(dǎo)函數(shù)為pp。1/13/202317（3）三次樣條插值方法

（4）樣條擬合方法（最小二乘法）

其中：函數(shù)csaps()、spap2()、spaps()和fnval()已在前文中介紹，

求導(dǎo)函數(shù)fnder()的調(diào)用格式為：pp=fnder(f,dorder)

該函數(shù)計(jì)算函數(shù)f的dorder階導(dǎo)函數(shù)，默認(rèn)階數(shù)dorder=1。離散數(shù)據(jù)三次樣條插值函數(shù)cscs的導(dǎo)數(shù)pppp在xi的導(dǎo)數(shù)值離散數(shù)據(jù)樣條擬合函數(shù)spsp的導(dǎo)數(shù)pppp在xi的導(dǎo)數(shù)值。

1/13/202318例5.3：反應(yīng)物A在一等溫間歇反應(yīng)器中發(fā)生的反應(yīng)為：A測(cè)量得到的反應(yīng)器中不同時(shí)間下反應(yīng)物A的濃度表5-4間歇反應(yīng)器動(dòng)力學(xué)數(shù)據(jù)系統(tǒng)的動(dòng)力學(xué)模型為：，試根據(jù)表中數(shù)據(jù)確定其反應(yīng)速率方程。產(chǎn)物如表5-4所示。1/13/202319解：（1）系統(tǒng)的動(dòng)力學(xué)模型為非線性形式，可將其線性化。對(duì)方程兩邊取對(duì)數(shù)：令則原模型變?yōu)椋海?）計(jì)算t=[0204060120180300];CA=[10865321];sp=spaps(t,CA,0.006);%生成光滑B樣條擬合函數(shù)pp=fnder(sp);%生成光滑B樣條擬合函數(shù)的導(dǎo)函數(shù)dCAdt=fnval(pp,t);%計(jì)算導(dǎo)函數(shù)在t處的數(shù)值微分%繪制圖形ti=linspace(t(1),t(end),200);CAi=fnval(sp,ti);plot(t,CA,'bo',ti,CAi,'r-'),xlabel('t'),ylabel('CA')及得到擬合曲線的圖形程序如下：1/13/202320（3）線性擬合程序如下：y=log(-dCAdt);x=log(CA);p=polyfit(x,y,1);k=exp(p(2)),m=p(1)執(zhí)行結(jié)果：k=0.0059m=1.2904所以本例的反應(yīng)速率方程為：1/13/2023211)被積函數(shù)以一組數(shù)據(jù)形式表示；2)被積函數(shù)過(guò)于特殊或原函數(shù)不能用初等函數(shù)表示，積分表中無(wú)法找到可沿用的現(xiàn)成公式；3)有的原函數(shù)十分復(fù)雜難以計(jì)算。對(duì)于積分：但是在工程技術(shù)和科學(xué)研究中,常會(huì)見到以下現(xiàn)象:0導(dǎo)入5.2數(shù)值積分1/13/202322

數(shù)值積分就是一種常用的近似計(jì)算方法。數(shù)值積分方法不受以上幾個(gè)問(wèn)題的限制，在化學(xué)化工領(lǐng)域應(yīng)用甚廣，如反應(yīng)熱效應(yīng)計(jì)算、熱容計(jì)算、熵的計(jì)算、反應(yīng)活化能的計(jì)算等。以上這些現(xiàn)象,Newton-Leibniz很難發(fā)揮作用，只能建立積分的近似計(jì)算方法1/13/202323如：某反應(yīng)器中進(jìn)行的反應(yīng)，計(jì)算出口物料的焓值：1/13/2023240.1數(shù)值積分的基本思路和方法

常用的數(shù)值積分的基本思路來(lái)自于插值法，它通過(guò)構(gòu)造一個(gè)插值多項(xiàng)式作為的近似表達(dá)式，用的積分值作為的近似積分值。數(shù)值積分的方法很豐富，常用的插值型求積公式有兩類：一類是等距節(jié)點(diǎn)的牛頓－柯特斯求積公式；另一類是不等距節(jié)點(diǎn)的高斯型求積公式。1/13/202325

Newton-Cotes公式是指等距節(jié)點(diǎn)下使用Lagrange插值多項(xiàng)式建立的數(shù)值求積公式各節(jié)點(diǎn)為則：§5.2.1牛頓－柯特斯（Newton-Cotes）求積公式1方法概述1/13/202326這里是既不依賴于被積函數(shù)，也不依賴于積分區(qū)間的常數(shù)，稱為柯特斯系數(shù)。式（5-3）稱為牛頓－柯特斯求積公式。其中：

一般地：令，得：

（5-3）1/13/202327在Newton-Cotes公式中,n=1,2,4時(shí)的公式是最常用也是最重要三個(gè)公式,稱為低階公式（當(dāng)n=0時(shí)的公式為矩形公式）1）梯形(trapezia)公式Cotes系數(shù)為求積公式為1/13/202328上式稱為梯形求積公式,（5-4）1/13/2023292）Simpson公式Cotes系數(shù)為1/13/202330上式稱為Simpson求積公式，也稱三點(diǎn)公式或拋物線公式求積公式為式（5-4）和式（5-5）是化工領(lǐng)域常用的兩個(gè)求積公式。與梯形法求積公式相比，Simpson法求積公式是一個(gè)較高精度的求積公式。用式（5-3）還可得到更高階數(shù)（5-5）1/13/202331

Newton-Cotes公式當(dāng)n大于7時(shí)，公式的穩(wěn)定性將無(wú)法保證，因此,在實(shí)際應(yīng)用中一般不使用高階Newton-Cotes公式而是采用低階復(fù)合求積法在實(shí)際計(jì)算中為了保證計(jì)算的精度，往往首先用分點(diǎn)xk=a+kh,

(k=0,1,…,n)將區(qū)間［a,b］分成n個(gè)相等的子區(qū)間，而后對(duì)每個(gè)子區(qū)間再應(yīng)用梯形公式或Simpson公式，分別得到：2復(fù)化法求積公式復(fù)化Simpson公式：1/13/202332盡管復(fù)化法求積公式具有很高的精度，但是它必須采用等步長(zhǎng)方法，從而限制了它的效率。這里我們介紹一種更加靈活選取步長(zhǎng)的方法，即自適應(yīng)步長(zhǎng)法。（5-8）。（5-6）

（5-7），令，當(dāng)

上Simpson積分達(dá)到精度時(shí)，可認(rèn)為區(qū)間取計(jì)算需滿足的精度為以Simpson積分法為例，某區(qū)間，記，考慮該區(qū)間上的Simpson積分和二等分以后的兩個(gè)Simpson積分和：3自適應(yīng)求積公式1/13/202333這樣重復(fù)下去，直至每個(gè)分段部分達(dá)到相應(yīng)精度（步長(zhǎng)為時(shí)精度為）。這樣，不同段的步長(zhǎng)可能是不一樣的，積分結(jié)果為每一小段的積分總和。自適應(yīng)步長(zhǎng)Simpson法從開始按式（5-6）～（5-8）的方法檢驗(yàn)。若滿足精度，則以為計(jì)算結(jié)果；若不滿足精度，則分成兩個(gè)小區(qū)間各自重復(fù)逐步上述過(guò)程，每個(gè)小區(qū)間精度為1/13/2023344Newton-Cotes求積公式的Matlab實(shí)現(xiàn)

常用的三種Newton-Cotes系列數(shù)值積分法的相應(yīng)Matlab函數(shù)如下:1)復(fù)合梯形法數(shù)值積分：trapz()調(diào)用形式：Z=trapz(X,Y)其中:X,Y－分別代表數(shù)目相同的向量或數(shù)組，而Y與X的關(guān)系可以是一函數(shù)型態(tài)（如y=sin(x)）或是不以函數(shù)描述的離散型態(tài)；

Z－代表返回的積分值；注:離散型態(tài)數(shù)據(jù)用trapz函數(shù)時(shí)，還需設(shè)定X在區(qū)間[a,b]之間離散點(diǎn)的間隔；缺省參數(shù)X時(shí)，表示X被等分，每份寬為1。1/13/2023352）自適應(yīng)Simpson法數(shù)值積分：quad（）基本調(diào)用格式：q=quad(fun,a,b)

或q=quad(fun,a,b,tol,trace,p1,p2,…)其中:fun－被積函數(shù)。可以是內(nèi)置函數(shù)、m文件或函數(shù)句柄,函數(shù)表達(dá)式中的必須使用點(diǎn)運(yùn)算符號(hào)。

a,b－分別是積分的下限和上限；

q－積分結(jié)果。

tol－默認(rèn)誤差限，默認(rèn)值為1.e-6.trace-取0表示不用圖形顯示積分過(guò)程，非0表示用圖形顯示積分過(guò)程

p1,p2,…直接傳遞給函數(shù)fun的參數(shù)。3）自適應(yīng)Lobatto法數(shù)值積分：quadl（）

quadl是高階的自適應(yīng)Newton-Cotes數(shù)值積分法函數(shù)，它比quad函數(shù)更有效，精度更高。其使用方法和quad（）完全相同。需要了解更多的內(nèi)容可查閱Matlab功能函數(shù)庫(kù)（funfun）。1/13/202336例5.4：真實(shí)氣體的逸度可用下式計(jì)算：

現(xiàn)測(cè)得00C下氫氣的有關(guān)數(shù)值如表5-5所示，試求1000atm下的逸度。表5-500C下氫氣的相關(guān)數(shù)值-真實(shí)氣體的實(shí)測(cè)體積和按理想氣體定律計(jì)算得到的體積之間的差值。R-氣體常數(shù)[]其中：-逸度；-壓力，atm；T-絕對(duì)溫度，K。

1/13/202337解：本題是離散型數(shù)據(jù)，可用trapz函數(shù)求解數(shù)值積分。（1）計(jì)算數(shù)值積分，程序如下：%計(jì)算數(shù)值積分P=0:100:1000;a1=[15.4615.4615.4615.6115.8515.9316.0916.1316.1616.1616.14];a1=-a1;A=trapz(P,a1);%梯形法計(jì)算數(shù)值積分A=-A執(zhí)行結(jié)果：A=15865（2）計(jì)算逸度值，程序如下：%計(jì)算逸度lf=log10(1000)+A./(2.303.*82.06.*273.2);%lf代表f=10.^lf執(zhí)行結(jié)果：f=2.0290e+003所以f=2029.00atm1/13/202338例5.5：氯仿-苯雙組分精餾系統(tǒng)的氣液平衡數(shù)據(jù)如表5-6所示。規(guī)定進(jìn)料和塔頂?shù)慕M成分別是

，精餾段的回流比為精餾段理論板數(shù)的模型為表5-6精餾段氣液平衡數(shù)據(jù)，試用Matlab計(jì)算所需的精餾段理論板數(shù)。1/13/202339解：因模型中的的函數(shù)關(guān)系是以表格形式給出，我們?nèi)粢眯疗丈ǖ扔?jì)算較精確的精餾段理論板數(shù)的時(shí)候，需先采用擬合法將離散數(shù)據(jù)（，）擬合成多項(xiàng)式，再將多項(xiàng)式代入被積函數(shù)求積。計(jì)算程序如下：functionli55clearall;xi=[0.1780.2750.3720.4560.6500.844];yi=[0.2430.3820.5180.6160.7950.931];sp=spline(xi,yi);%用spline（）擬合多項(xiàng)式%畫出擬合曲線，直觀比較擬合效果xplot=linspace(xi(1),xi(end),100);yplot=fnval(sp,xplot);plot(xi,yi,'o',xplot,yplot,'-');N=quad(@func1,0.4,0.9,[],[],sp);%計(jì)算精餾段理論板數(shù)N=round(N+0.5)%數(shù)據(jù)取整functionf=func1(x,sp)y=fnval(sp,x);f=1./(y-x-(0.9-y)./5);執(zhí)行結(jié)果：N=

5所以，精餾段理論板數(shù)為5塊。1/13/202340§5.2.2高斯－勒讓德公式

高斯法是一種精度較高的求積分法，它的一般公式是 式中ω(x)是一個(gè)權(quán)重函數(shù)，Aj為系數(shù)，xj為橫坐標(biāo)上的節(jié)點(diǎn) 高斯－勒讓德多項(xiàng)式的權(quán)重函數(shù)ω(x)=1，因而，一個(gè)n點(diǎn)高斯－勒讓德求積公式具有如下形式：

1方法概述

1/13/202341右邊的f(xj)是函數(shù)f(x)在節(jié)點(diǎn)xj處的值，節(jié)點(diǎn)xj是勒讓德多項(xiàng)式Pn(x)的根。

Aj為系數(shù)，其值為式中P’n(x)是勒讓德多項(xiàng)式Pn(x)的一階導(dǎo)數(shù)1/13/202342pn(x)的前幾項(xiàng)表達(dá)式為

由高斯－勒讓德多項(xiàng)式得出的2～6點(diǎn)