第5章

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

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

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

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

求導得：

（5-2）;式（5-1）和（5-2）不僅合用于求節點處的導數，并且可求非節點處的導數。（5-1）其中，6/22/2025112三次樣條插值函數求微分的Matlab函數Matlab求離散數據的三次樣條插值函數微分辦法分三個環節：Step1:對離散數據用csapi函數（或spline函數）得到其三次樣條插值函數調用形式pp=csapi(x,y)其中：x,y分別為離散數據對的自變量和因變量；pp為得到的三次樣條插值函數。Step2:用fnder函數求三次樣條插值函數的導數調用形式fprime=fnder(f,dorder)其中：f為三次樣條插值函數；dorder為三次樣條插值函數的求導階數；fprime為得到的三次樣條插值函數的導函數。Step3:可用fnval函數求導函數在未知點處的導數值調用形式v=fnval(fprime,x)其中：fprime為三次樣條插值函數導函數；x為未知點處自變量值；v為未知點處的導數值。6/22/202512例5.2：某液體冷卻時，溫度隨時間的變化數據如表5-3所示：表5-3冷卻溫度隨時間的變化數據試分別計算t＝2，3，4min及t＝1.5，2.5，4.5min時的降溫速率。解：三次樣條插值函數求數值微分的程序以下：t=[0:5];T=[92,85.3,79.5,74.5,70.2,67];cs=csapi(t,T);%生成三次樣條插值函數pp=fnder(cs);%生成三次樣條插值函數的導函數t1=[2,3,4,1.5,2.5,4.5];dT=fnval(pp,t1);%計算導函數在t1處的導數值disp('對應時間時的降溫速率:')disp([t1;dT])執行成果：對應時間時的降溫速率: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注：前者是計算節點處的一階導數，后者是計算非節點處的一階導數6/22/202513§5.1.3最小二乘法樣條擬合函數求微分

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

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

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

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

求導函數fnder()的調用格式為：pp=fnder(f,dorder)該函數計算函數f的dorder階導函數，默認階數dorder=1。離散數據三次樣條插值函數cscs的導數pppp在xi的導數值離散數據樣條擬合函數spsp的導數pppp在xi的導數值。

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

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

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

（5-7），令，當

上Simpson積分達到精度時，可認為區間取計算需滿足的精度為以Simpson積分法為例，某區間，記，考慮該區間上的Simpson積分和二等分以后的兩個Simpson積分和：3自適應求積公式6/22/202533這樣重復下去，直至每個分段部分達到相應精度（步長為時精度為）。這樣，不同段的步長可能是不一樣的，積分結果為每一小段的積分總和。自適應步長Simpson法從開始按式（5-6）～（5-8）的方法檢驗。若滿足精度，則以為計算結果；若不滿足精度，則分成兩個小區間各自重復逐步上述過程，每個小區間精度為6/22/2025344Newton-Cotes求積公式的Matlab實現

慣用的三種Newton-Cotes系列數值積分法的對應Matlab函數以下:1)復合梯形法數值積分：trapz()調用形式：Z=trapz(X,Y)其中:X,Y－分別代表數目相似的向量或數組，而Y與X的關系能夠是一函數型態（如y=sin(x)）或是不以函數描述的離散型態；Z－代表返回的積分值；注:離散型態數據用trapz函數時，還需設定X在區間[a,b]之間離散點的間隔；缺省參數X時，表達X被等分，每份寬為1。6/22/2025352）自適應Simpson法數值積分：quad（）基本調用格式：q=quad(fun,a,b)或q=quad(fun,a,b,tol,trace,p1,p2,…)其中:fun－被積函數。能夠是內置函數、m文獻或函數句柄,函數體現式中的必須使用點運算符號。a,b－分別是積分的下限和上限；q－積分成果。tol－默認誤差限，默認值為1.e-6.trace-取0表達不用圖形顯示積分過程，非0表達用圖形顯示積分過程p1,p2,…直接傳遞給函數fun的參數。3）自適應Lobatto法數值積分：quadl（）quadl是高階的自適應Newton-Cotes數值積分法函數，它比quad函數更有效，精度更高。其使用辦法和quad（）完全相似。需要理解更多的內容可查閱Matlab功效函數庫（funfun）。6/22/202536例5.4：真實氣體的逸度可用下式計算：

現測得00C下氫氣的有關數值如表5-5所示，試求1000atm下的逸度。表5-500C下氫氣的有關數值（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-氣體常數[]其中：-逸度；-壓力，atm；T-絕對溫度，K。

6/22/202537解：本題是離散型數據，可用trapz函數求解數值積分。（1）計算數值積分，程序以下：%計算數值積分P=0:100:1000;a1=[15.4615.4615.4615.6115.8515.9316.0916.1316.1616.1616.14];a1=-a1;A=trapz(P,a1);%梯形法計算數值積分A=-A執行成果：A=15865（2）計算逸度值，程序以下：%計算逸度lf=log10(1000)+A./(2.303.*82.06.*273.2);%lf代表f=10.^lf執行成果：f=2.0290e+003因此f=2029.00atm6/22/202538例5.5：氯仿-苯雙組分精餾系統的氣液平衡數據如表5-6所示。規定進料和塔頂的構成分別是，精餾段的回流比為精餾段理論板數的模型為表5-6精餾段氣液平衡數據0.1780.2750.3720.4560.6500.8440.2430.3820.5180.6160.7950.931，試用Matlab計算所需的精餾段理論板數。6/22/202539解：因模型中的的函數關系是以表格形式給出，我們若要用辛普森法等計算較精確的精餾段理論板數的時候，需先采用擬正當將離散數據（，）擬合成多項式，再將多項式代入被積函數求積。計算程序以下：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);%計算精餾段理論板數N=round(N+0.5)%數據取整functionf=func1(x,sp)y=fnval(sp,x);f=1./(y-x-(0.9-y)./5);執行成果：N=5因此，精餾段理論板數為5塊。6/22/202540§5.2.2高斯－勒讓德公式 高斯法是一種精度較高的求積分法，它的普通公式是 式中ω(x)是一種權重函數，Aj為系數，xj為橫坐標上的節點 高斯－勒讓德多項式的權重函數ω(x)=1，因而，一種n點高斯－勒讓德求積公式含有以下形式：

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

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

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