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

第5章

數(shù)值微分與數(shù)值積分6/22/20251微積分是高等數(shù)學(xué)中的重要內(nèi)容，在化學(xué)工程上有許多非常重要的應(yīng)用微積分的數(shù)值辦法，不同于高等數(shù)學(xué)中的解析辦法，特別適合求解沒有或很難求出微分或積分體現(xiàn)式的實際化工問題的計算，例如：列表函數(shù)求微分或積分引言1---數(shù)值微積分辦法不同于解析辦法6/22/20252數(shù)值微分和數(shù)值積分與插值和擬合往往是密不可分的如在進行數(shù)值微分時，常針對離散的數(shù)據(jù)點，運用插值和擬合能夠減少誤差；而數(shù)值積分的基本思路也來自于插值法。例如如果所積函數(shù)的形式比較復(fù)雜或以表格形式給出，則可通過構(gòu)造一種插值多項式來替代原函數(shù)，從而使問題簡化引言2---數(shù)值微分和數(shù)值積分與插值和擬合關(guān)系親密6/22/20253§5.1數(shù)值微分

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

t1t2…tnpA1pA2…PAn0導(dǎo)入6/22/202540.1建立數(shù)值微分公式的三種思路慣用三種思路建立數(shù)值微分公式：從微分定義出發(fā)，通過近似解決，得到數(shù)值微分的近似公式從插值近似公式出發(fā)，對插值公式的近似求導(dǎo)可得到數(shù)值微分的近似公式先用最小二乘擬合辦法根據(jù)已知數(shù)據(jù)得近似函數(shù)，再對此近似函數(shù)求微分可得到數(shù)值微分的近似公式。然后對各辦法數(shù)值微分后得到的多項式求值，即可求出任意點處的任意階微分6/22/202551辦法概述在微積分中，一階微分的計算能夠在二相鄰點x+h和x間函數(shù)取下列極限求得：取其達成極限前的形式，就得到下列微分的差分近似式：注：高階微分項能夠運用低階微分項來計算，如二階微分式能夠表達為：對應(yīng)的差分式有：§5.1.1差分近似微分上式中三種不同表達形式依次是一階前向差分、一階后向差分和一階中心差分來近似表達微分。其中一階中心差分的精度較高。6/22/202562差分的Matlab實現(xiàn)在Matlab中，可用diff函數(shù)進行離散數(shù)據(jù)的近似求導(dǎo)調(diào)用形式：Y=diff(X,n)其中：X表達求導(dǎo)變量，能夠是向量或矩陣。如是矩陣形式則按各列作差分；n表達n階差分，即差分n次Y是X的差分成果注：用diff函數(shù)進行離散數(shù)據(jù)的近似求導(dǎo)與前向差分近似，但誤差較大。最佳將數(shù)據(jù)運用插值或擬合得到多項式，然后對近似多項式進行微分6/22/20257例5.1：丁二烯的氣相二聚反映以下：實驗在一定容器的反映器中進行，3260C時，測得物系中丁二烯的分壓(mmHg)與時間的關(guān)系如表5-2所示。。表5-2丁二烯二聚反映實驗數(shù)據(jù)t(min)(mmHg)t(min)(mmHg)0632.050362.05590.055348.010552.060336.015515.065325.020485.070314.025458.075304.030435.080294.035414.085284.040396.090274.045378.0用數(shù)值微分法計算所列時刻每一瞬間的反映速率6/22/20258解：程序以下：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);％求時間t的差分dpA=diff(pA);％求壓力的差分q=dpA./dt％q為數(shù)值微分成果執(zhí)行成果：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.00006/22/20259§5.1.2三次樣條插值函數(shù)求微分若三次樣條插值函數(shù)收斂于，那么導(dǎo)數(shù)收斂于，因此用樣條插值函數(shù)作為函數(shù)，不但彼此的函數(shù)值非常接近，而且導(dǎo)數(shù)值也很接近。

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

求導(dǎo)得：

（5-2）;式（5-1）和（5-2）不僅合用于求節(jié)點處的導(dǎo)數(shù)，并且可求非節(jié)點處的導(dǎo)數(shù)。（5-1）其中，6/22/2025112三次樣條插值函數(shù)求微分的Matlab函數(shù)Matlab求離散數(shù)據(jù)的三次樣條插值函數(shù)微分辦法分三個環(huán)節(jié)：Step1:對離散數(shù)據(jù)用csapi函數(shù)（或spline函數(shù)）得到其三次樣條插值函數(shù)調(diào)用形式pp=csapi(x,y)其中：x,y分別為離散數(shù)據(jù)對的自變量和因變量；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ù)在未知點處的導(dǎo)數(shù)值調(diào)用形式v=fnval(fprime,x)其中：fprime為三次樣條插值函數(shù)導(dǎo)函數(shù)；x為未知點處自變量值；v為未知點處的導(dǎo)數(shù)值。6/22/202512例5.2：某液體冷卻時，溫度隨時間的變化數(shù)據(jù)如表5-3所示：表5-3冷卻溫度隨時間的變化數(shù)據(jù)試分別計算t＝2，3，4min及t＝1.5，2.5，4.5min時的降溫速率。解：三次樣條插值函數(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);%計算導(dǎo)函數(shù)在t1處的導(dǎo)數(shù)值disp('對應(yīng)時間時的降溫速率:')disp([t1;dT])執(zhí)行成果：對應(yīng)時間時的降溫速率:2.00003.00004.00001.50002.50004.5000-5.3722-4.6722-3.8389-5.7972-4.9889-3.2222t(min)012345T(0C)92.085.379.574.570.267.0注：前者是計算節(jié)點處的一階導(dǎo)數(shù)，后者是計算非節(jié)點處的一階導(dǎo)數(shù)6/22/202513§5.1.3最小二乘法樣條擬合函數(shù)求微分

在實際化工應(yīng)用中，當(dāng)來自實驗觀察的離散數(shù)據(jù)不可避免地含有較大隨機誤差時，此時用插值公式求數(shù)值微分即使樣本點處誤差較小，但可能會使非樣本點處產(chǎn)生較大誤差為此，可采用最小二乘法樣條擬合實驗數(shù)據(jù)，獲得一種函數(shù)模型，然后再對其求導(dǎo)數(shù)。與樣條插值不同，樣條擬合不規(guī)定曲線通過全部的數(shù)據(jù)點，這樣解決求導(dǎo)成果會有很大改善6/22/2025141辦法概述

可用最小二乘法擬合成平滑B樣條曲線，即對于離散數(shù)據(jù)（所求的K次樣條擬合函數(shù)滿足：其中：再對擬合函數(shù)作平滑解決后求導(dǎo)，即可求出任意點處微分。）為權(quán)重系數(shù)，默認為16/22/202515Matlab求離散數(shù)據(jù)的最小二乘法平滑B樣條擬合函數(shù)求微分共三個環(huán)節(jié)：Step1:對離散數(shù)據(jù)用spaps函數(shù)得到最小二乘平滑B樣條擬合函數(shù)。調(diào)用格式：sp=spaps(x,y,tol)其中：x,y――要解決的離散數(shù)據(jù)（）tol――光滑時的允許精度，普通?。?0-2－10-4）Step2:可用fnder函數(shù)求最小二乘平滑B樣條擬合函數(shù)的導(dǎo)數(shù);Step3:可用fnval函數(shù)求導(dǎo)函數(shù)在未知點處的導(dǎo)數(shù)值。fnder()和fnval（）調(diào)用形式前文中已經(jīng)介紹過。2最小二乘法平滑B樣條擬合函數(shù)求微分的Matlab實現(xiàn)6/22/202516由離散數(shù)據(jù)求數(shù)值微分的四種辦法及有關(guān)MATLAB函數(shù)：（1）差分法用差分函數(shù)diff()近似計算導(dǎo)數(shù)，即dy=diff(y)./diff(x)。對于向量X，diff(X)表達了[X(2)-X(1)X(3)-X(2)...X(n)-X(n-1)].對于矩陣X，diff(X)表達了[X(2:n,:)-X(1:n-1,:)]小結(jié)注：此法用一階差分，精度較差，若改用二階差分，可大大提高精度，但編程麻煩些。（2）多項式擬合辦法向量p表達的多項式擬合函數(shù)導(dǎo)函數(shù)pppp在xi的導(dǎo)數(shù)值。其中：函數(shù)polyfit()和polyval()在前文中已介紹導(dǎo)函數(shù)polyder()的調(diào)用格式為：pp=polyder(p)

離散數(shù)據(jù)該函數(shù)對向量p表達的多項式函數(shù)進行求導(dǎo)，返回導(dǎo)函數(shù)為pp。6/22/202517（3）三次樣條插值辦法（4）樣條擬合辦法（最小二乘法）

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

求導(dǎo)函數(shù)fnder()的調(diào)用格式為：pp=fnder(f,dorder)該函數(shù)計算函數(shù)f的dorder階導(dǎo)函數(shù)，默認階數(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ù)值。

6/22/202518例5.3：反映物A在一等溫間歇反映器中發(fā)生的反映為：A測量得到的反映器中不同時間下反映物A的濃度表5-4間歇反映器動力學(xué)數(shù)據(jù)t(s)0204060120180300(mol/L)10865321系統(tǒng)的動力學(xué)模型為：，試根據(jù)表中數(shù)據(jù)擬定其反映速率方程。產(chǎn)物如表5-4所示。6/22/202519解：（1）系統(tǒng)的動力學(xué)模型為非線性形式，可將其線性化。對方程兩邊取對數(shù)：令則原模型變?yōu)椋海?）計算t=[0204060120180300];CA=[10865321];sp=spaps(t,CA,0.006);%生成光滑B樣條擬合函數(shù)pp=fnder(sp);%生成光滑B樣條擬合函數(shù)的導(dǎo)函數(shù)dCAdt=fnval(pp,t);%計算導(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')及得到擬合曲線的圖形程序以下：6/22/202520（3）線性擬合程序以下：y=log(-dCAdt);x=log(CA);p=polyfit(x,y,1);k=exp(p(2)),m=p(1)執(zhí)行成果：k=0.0059m=1.2904因此本例的反映速率方程為：6/22/2025211)被積函數(shù)以一組數(shù)據(jù)形式表達；2)被積函數(shù)過于特殊或原函數(shù)不能用初等函數(shù)表達，積分表中無法找到可沿用的現(xiàn)成公式；3)有的原函數(shù)十分復(fù)雜難以計算。對于積分：但是在工程技術(shù)和科學(xué)研究中,常會見到下列現(xiàn)象:0導(dǎo)入5.2數(shù)值積分6/22/202522 數(shù)值積分就是一種慣用的近似計算辦法。數(shù)值積分辦法不受以上幾個問題的限制，在化學(xué)化工領(lǐng)域應(yīng)用甚廣，如反映熱效應(yīng)計算、熱容計算、熵的計算、反映活化能的計算等。以上這些現(xiàn)象,Newton-Leibniz很難發(fā)揮作用，只能建立積分的近似計算辦法6/22/202523如：某反映器中進行的反映，計算出口物料的焓值：6/22/2025240.1數(shù)值積分的基本思路和辦法慣用的數(shù)值積分的基本思路來自于插值法，它通過構(gòu)造一種插值多項式作為的近似體現(xiàn)式，用的積分值作為的近似積分值。數(shù)值積分的辦法很豐富，慣用的插值型求積公式有兩類：一類是等距節(jié)點的牛頓－柯特斯求積公式；另一類是不等距節(jié)點的高斯型求積公式。6/22/202525

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

（5-3）6/22/202527在Newton-Cotes公式中,n=1,2,4時的公式是最慣用也是最重要三個公式,稱為低階公式（當(dāng)n=0時的公式為矩形公式）1）梯形(trapezia)公式Cotes系數(shù)為求積公式為6/22/202528上式稱為梯形求積公式,（5-4）6/22/2025292）Simpson公式Cotes系數(shù)為6/22/202530上式稱為Simpson求積公式，也稱三點公式或拋物線公式求積公式為式（5-4）和式（5-5）是化工領(lǐng)域慣用的兩個求積公式。與梯形法求積公式相比，Simpson法求積公式是一種較高精度的求積公式。用式（5-3）還可得到更高階數(shù)（5-5）6/22/202531Newton-Cotes公式當(dāng)n不不大于7時，公式的穩(wěn)定性將無法確保，因此,在實際應(yīng)用中普通不使用高階Newton-Cotes公式而是采用低階復(fù)合求積法在實際計算中為了確保計算的精度，往往首先用分點xk=a+kh,(k=0,1,…,n)將區(qū)間［a,b］分成n個相等的子區(qū)間，而后對每個子區(qū)間再應(yīng)用梯形公式或Simpson公式，分別得到：2復(fù)化法求積公式復(fù)化Simpson公式：6/22/202532盡管復(fù)化法求積公式含有很高的精度，但是它必須采用等步長辦法，從而限制了它的效率。這里我們介紹一種更加靈活選用步長的辦法，即自適應(yīng)步長法。（5-8）。（5-6）

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

上Simpson積分達到精度時，可認為區(qū)間取計算需滿足的精度為以Simpson積分法為例，某區(qū)間，記，考慮該區(qū)間上的Simpson積分和二等分以后的兩個Simpson積分和：3自適應(yīng)求積公式6/22/202533這樣重復(fù)下去，直至每個分段部分達到相應(yīng)精度（步長為時精度為）。這樣，不同段的步長可能是不一樣的，積分結(jié)果為每一小段的積分總和。自適應(yīng)步長Simpson法從開始按式（5-6）～（5-8）的方法檢驗。若滿足精度，則以為計算結(jié)果；若不滿足精度，則分成兩個小區(qū)間各自重復(fù)逐步上述過程，每個小區(qū)間精度為6/22/2025344Newton-Cotes求積公式的Matlab實現(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è)定X在區(qū)間[a,b]之間離散點的間隔；缺省參數(shù)X時，表達X被等分，每份寬為1。6/22/2025352）自適應(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ù)體現(xiàn)式中的必須使用點運算符號。a,b－分別是積分的下限和上限；q－積分成果。tol－默認誤差限，默認值為1.e-6.trace-取0表達不用圖形顯示積分過程，非0表達用圖形顯示積分過程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ù)庫（funfun）。6/22/202536例5.4：真實氣體的逸度可用下式計算：

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

6/22/202537解：本題是離散型數(shù)據(jù)，可用trapz函數(shù)求解數(shù)值積分。（1）計算數(shù)值積分，程序以下：%計算數(shù)值積分P=0:100:1000;a1=[15.4615.4615.4615.6115.8515.9316.0916.1316.1616.1616.14];a1=-a1;A=trapz(P,a1);%梯形法計算數(shù)值積分A=-A執(zhí)行成果：A=15865（2）計算逸度值，程序以下：%計算逸度lf=log10(1000)+A./(2.303.*82.06.*273.2);%lf代表f=10.^lf執(zhí)行成果：f=2.0290e+003因此f=2029.00atm6/22/202538例5.5：氯仿-苯雙組分精餾系統(tǒng)的氣液平衡數(shù)據(jù)如表5-6所示。規(guī)定進料和塔頂?shù)臉?gòu)成分別是，精餾段的回流比為精餾段理論板數(shù)的模型為表5-6精餾段氣液平衡數(shù)據(jù)0.1780.2750.3720.4560.6500.8440.2430.3820.5180.6160.7950.931，試用Matlab計算所需的精餾段理論板數(shù)。6/22/202539解：因模型中的的函數(shù)關(guān)系是以表格形式給出，我們?nèi)粢眯疗丈ǖ扔嬎爿^精確的精餾段理論板數(shù)的時候，需先采用擬正當(dāng)將離散數(shù)據(jù)（，）擬合成多項式，再將多項式代入被積函數(shù)求積。計算程序以下：functionli55clearall;xi=[0.1780.2750.3720.4560.6500.844];yi=[0.2430.3820.5180.6160.7950.931];sp=spline(xi,yi);%用spline（）擬合多項式%畫出擬合曲線，直觀比較擬合效果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);%計算精餾段理論板數(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í)行成果：N=5因此，精餾段理論板數(shù)為5塊。6/22/202540§5.2.2高斯－勒讓德公式 高斯法是一種精度較高的求積分法，它的普通公式是 式中ω(x)是一種權(quán)重函數(shù)，Aj為系數(shù)，xj為橫坐標上的節(jié)點 高斯－勒讓德多項式的權(quán)重函數(shù)ω(x)=1，因而，一種n點高斯－勒讓德求積公式含有以下形式：

1辦法概述6/22/202541右邊的f(xj)是函數(shù)f(x)在節(jié)點xj處的值，節(jié)點xj是勒讓德多項式Pn(x)的根。

Aj為系數(shù)，其值為式中P’n(x)是勒讓德多項式Pn(x)的一階導(dǎo)數(shù)6/22/202542pn(x)的前幾項體現(xiàn)式為

由高斯－勒讓德多項式得出的2～6點的根和系數(shù)見表5-76/22/202543表5-72～6點的高斯節(jié)點和高斯求積系數(shù)n2±0.577351.000003±0.774600.000000.555560.888894±0.86114±0.339980.347850.65