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

1、2021-11-151第第5章章 數(shù)值微分與數(shù)值積分?jǐn)?shù)值微分與數(shù)值積分22021-11-15n微積分是高等數(shù)學(xué)中的重要內(nèi)容，在化學(xué)微積分是高等數(shù)學(xué)中的重要內(nèi)容，在化學(xué)工程上有許多非常重要的應(yīng)用工程上有許多非常重要的應(yīng)用n微積分的數(shù)值方法，不同于高等數(shù)學(xué)中的微積分的數(shù)值方法，不同于高等數(shù)學(xué)中的解析方法，尤其適合求解沒(méi)有或很難求出解析方法，尤其適合求解沒(méi)有或很難求出微分或積分表達(dá)式的實(shí)際化工問(wèn)題的計(jì)算，微分或積分表達(dá)式的實(shí)際化工問(wèn)題的計(jì)算，例如：列表函數(shù)求微分或積分例如：列表函數(shù)求微分或積分 引言引言1-數(shù)值微積分方法不同于解析方法數(shù)值微積分方法不同于解析方法32021-11-15n數(shù)值微分和數(shù)值

2、積分與插值和擬合往往是密數(shù)值微分和數(shù)值積分與插值和擬合往往是密不可分的不可分的n如在進(jìn)行數(shù)值微分時(shí)，常針對(duì)離散的數(shù)據(jù)點(diǎn)，如在進(jìn)行數(shù)值微分時(shí)，常針對(duì)離散的數(shù)據(jù)點(diǎn)，利用插值和擬合可以減少誤差；而數(shù)值積分利用插值和擬合可以減少誤差；而數(shù)值積分的基本思路也來(lái)自于插值法。例如如果所積的基本思路也來(lái)自于插值法。例如如果所積函數(shù)的形式比較復(fù)雜或以表格形式給出，則函數(shù)的形式比較復(fù)雜或以表格形式給出，則可通過(guò)構(gòu)造一個(gè)插值多項(xiàng)式來(lái)代替原函數(shù)，可通過(guò)構(gòu)造一個(gè)插值多項(xiàng)式來(lái)代替原函數(shù)，從而使問(wèn)題簡(jiǎn)化從而使問(wèn)題簡(jiǎn)化引言引言2-數(shù)值微分和數(shù)值積分與插值和擬合數(shù)值微分和數(shù)值積分與插值和擬合關(guān)系密切關(guān)系密切42021-11-1

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

4、1t2tnpa1pa2pan0 0 導(dǎo)入導(dǎo)入52021-11-150.1 0.1 建立數(shù)值微分公式的三種思路建立數(shù)值微分公式的三種思路 常用三種思路建立數(shù)值微分公式：常用三種思路建立數(shù)值微分公式：1.從微分定義出發(fā)，通過(guò)近似處理，得到數(shù)值微從微分定義出發(fā)，通過(guò)近似處理，得到數(shù)值微分的近似公式分的近似公式2.從插值近似公式出發(fā)，對(duì)插值公式的近似求導(dǎo)從插值近似公式出發(fā)，對(duì)插值公式的近似求導(dǎo)可得到數(shù)值微分的近似公式可得到數(shù)值微分的近似公式3.先用最小二乘擬合方法根據(jù)已知數(shù)據(jù)得近似函先用最小二乘擬合方法根據(jù)已知數(shù)據(jù)得近似函數(shù)，再對(duì)此近似函數(shù)求微分可得到數(shù)值微分的數(shù)，再對(duì)此近似函數(shù)求微分可得到數(shù)值微分的

5、近似公式。然后對(duì)各方法數(shù)值微分后得到的多近似公式。然后對(duì)各方法數(shù)值微分后得到的多項(xiàng)式求值，即可求出任意點(diǎn)處的任意階微分項(xiàng)式求值，即可求出任意點(diǎn)處的任意階微分62021-11-15000()()( )()( )( )()22( )limlimlimhhhhhf xf xdf xf x hf xf xf x hf xdxhhh()()( )()( )( )()22( )hhf xf xdf xf xhf xf xf xhfxdxhhh( )( )dfxfxdx()()()( )( )()22( )hhfxfxfxhfxfxfxhfxhhh1 1 方法概述方法概述在微積分中，一階微分的計(jì)算可以在二相

6、鄰點(diǎn)x+h和x間函數(shù)取下列極限求得： 取其達(dá)到極限前的形式取其達(dá)到極限前的形式，就得到以下微分的差分近似式：，就得到以下微分的差分近似式： 注：注：高階微分項(xiàng)可以利用低階微分項(xiàng)來(lái)計(jì)算，如二階微分式可以表示為：高階微分項(xiàng)可以利用低階微分項(xiàng)來(lái)計(jì)算，如二階微分式可以表示為： 對(duì)應(yīng)的差分式有：對(duì)應(yīng)的差分式有：5.1.1 差分近似微分差分近似微分上式中三種不同表示形式依次是一階前向差分、一階后向差分和一階中上式中三種不同表示形式依次是一階前向差分、一階后向差分和一階中心差分來(lái)近似表示微分。其中一階中心差分的精度較高。心差分來(lái)近似表示微分。其中一階中心差分的精度較高。72021-11-152 2 差分的差

7、分的matlabmatlab實(shí)現(xiàn)實(shí)現(xiàn)n在在matlab中，可用中，可用diff函數(shù)進(jìn)行離散數(shù)據(jù)的近似函數(shù)進(jìn)行離散數(shù)據(jù)的近似求導(dǎo)求導(dǎo)n調(diào)用形式調(diào)用形式：y = diff(x,n)n其中：其中：x表示求導(dǎo)變量，可以是向量或矩陣。如表示求導(dǎo)變量，可以是向量或矩陣。如是矩陣形式則按各列作差分；是矩陣形式則按各列作差分； n表示表示n階差分，即差分階差分，即差分n次次 y是是x的差分結(jié)果的差分結(jié)果n注：注：用用diff函數(shù)進(jìn)行離散數(shù)據(jù)的近似求導(dǎo)與前向差函數(shù)進(jìn)行離散數(shù)據(jù)的近似求導(dǎo)與前向差分近似，分近似， 但誤差較大。最好將數(shù)據(jù)利用插值或擬但誤差較大。最好將數(shù)據(jù)利用插值或擬合得到多項(xiàng)式，然后對(duì)近似多項(xiàng)式進(jìn)行

8、微分合得到多項(xiàng)式，然后對(duì)近似多項(xiàng)式進(jìn)行微分82021-11-15464622()c hc h apadpdt例例5.1：丁二烯的氣相二聚反應(yīng)如下：丁二烯的氣相二聚反應(yīng)如下：實(shí)驗(yàn)在一定容器的反應(yīng)器中進(jìn)行，實(shí)驗(yàn)在一定容器的反應(yīng)器中進(jìn)行，3260c時(shí)，測(cè)得物系中丁二烯的分時(shí)，測(cè)得物系中丁二烯的分壓壓 (mmhg)與時(shí)間的關(guān)系如表與時(shí)間的關(guān)系如表5-2所示。所示。apap表表 5-2 丁二烯二聚反應(yīng)實(shí)驗(yàn)數(shù)據(jù)丁二烯二聚反應(yīng)實(shí)驗(yàn)數(shù)據(jù)t(min) (mmhg) t(min) (mmhg)0632.050362.05590.055348.010552.060336.015515.065325.020485.0

9、70314.025458.075304.030435.080294.035414.085284.040396.090274.045378.0用數(shù)值微分法計(jì)算所列時(shí)刻每一瞬間的反應(yīng)速率用數(shù)值微分法計(jì)算所列時(shí)刻每一瞬間的反應(yīng)速率 92021-11-15解：解：程序如下：t=0:5:90;pa=632.0 590.0 552.0 515.0 485.0 458.0 435.0 414.0 396.0 378.0 362.0 348.0 336.0 325.0 314.0 304.0 294.0 284.0 274.0;dt=diff(t); 求時(shí)間t的差分dpa=diff(pa); 求壓力的差分q=

10、dpa./dt q為數(shù)值微分結(jié)果執(zhí)行結(jié)果：執(zhí)行結(jié)果：q = columns 1 through 8 -8.4000 -7.6000 -7.4000 -6.0000 -5.4000 -4.6000 -4.2000 -3.6000 columns 9 through 16 -3.6000 -3.2000 -2.8000 -2.4000 -2.2000 -2.2000 -2.0000 -2.0000 columns 17 through 18 -2.0000 -2.0000102021-11-155.1.2 三次樣條插值函數(shù)求微分三次樣條插值函數(shù)求微分( )s x( )f x( )s x( )fx(

11、 )s x( )f x若三次樣條插值函數(shù)若三次樣條插值函數(shù)收斂于收斂于，那么導(dǎo)數(shù)，那么導(dǎo)數(shù)收斂于收斂于，因此用樣條插值函數(shù)，因此用樣條插值函數(shù)作為作為函數(shù)，不但彼此的函數(shù)值非常接近，而且導(dǎo)數(shù)值也很接近。函數(shù)，不但彼此的函數(shù)值非常接近，而且導(dǎo)數(shù)值也很接近。 用三次樣條插值函數(shù)求數(shù)值導(dǎo)數(shù)是可靠的，這是化工計(jì)用三次樣條插值函數(shù)求數(shù)值導(dǎo)數(shù)是可靠的，這是化工計(jì)算中求數(shù)值微分的有效方法。算中求數(shù)值微分的有效方法。的近似的近似112021-11-15( )( )fxs x1,2,in1,iixxx11()()( )( )iiiiiixxxxfxs xmmhh221111()()( )( )()226iiii

12、iiiiiiiiyyxxxxhfxs xmmmmhhh 1 1 方法概述方法概述用三次樣條插值函數(shù)建立的數(shù)值微分公式為：用三次樣條插值函數(shù)建立的數(shù)值微分公式為： 求導(dǎo)得：求導(dǎo)得： （5-2）;式（式（5-1）和（）和（5-2）不但適用于求節(jié)點(diǎn)處的導(dǎo)數(shù)，而且可求非節(jié)點(diǎn)處的導(dǎo)數(shù)。）不但適用于求節(jié)點(diǎn)處的導(dǎo)數(shù)，而且可求非節(jié)點(diǎn)處的導(dǎo)數(shù)。（5-1）其中，其中，122021-11-152 2 三次樣條插值函數(shù)求微分的三次樣條插值函數(shù)求微分的matlabmatlab函數(shù)函數(shù)matlab求離散數(shù)據(jù)的三次樣條插值函數(shù)微分方法分三個(gè)步驟：求離散數(shù)據(jù)的三次樣條插值函數(shù)微分方法分三個(gè)步驟：step 1:對(duì)離散數(shù)據(jù)用對(duì)離

13、散數(shù)據(jù)用csapi函數(shù)（或函數(shù)（或spline函數(shù)）得到其三次樣條插值函數(shù)函數(shù)）得到其三次樣條插值函數(shù)調(diào)用形式調(diào)用形式 pp = csapi(x,y)其中：其中： x,y分別為離散數(shù)據(jù)對(duì)的自變量和因變量；分別為離散數(shù)據(jù)對(duì)的自變量和因變量； pp為得到的三次樣條插值函數(shù)。為得到的三次樣條插值函數(shù)。step 2: 用用fnder函數(shù)求三次樣條插值函數(shù)的導(dǎo)數(shù)函數(shù)求三次樣條插值函數(shù)的導(dǎo)數(shù)調(diào)用形式調(diào)用形式 fprime = fnder(f,dorder)其中：其中： f為三次樣條插值函數(shù)；為三次樣條插值函數(shù)； dorder為三次樣條插值函數(shù)的求導(dǎo)階數(shù)；為三次樣條插值函數(shù)的求導(dǎo)階數(shù)； fprime為得到的

14、三次樣條插值函數(shù)的導(dǎo)函數(shù)。為得到的三次樣條插值函數(shù)的導(dǎo)函數(shù)。step3:可用可用fnval函數(shù)求導(dǎo)函數(shù)在未知點(diǎn)處的導(dǎo)數(shù)值函數(shù)求導(dǎo)函數(shù)在未知點(diǎn)處的導(dǎo)數(shù)值調(diào)用形式調(diào)用形式 v = fnval(fprime,x)其中：其中：fprime為三次樣條插值函數(shù)導(dǎo)函數(shù)；為三次樣條插值函數(shù)導(dǎo)函數(shù)； x為未知點(diǎn)處自變量值；為未知點(diǎn)處自變量值； v為未知點(diǎn)處的導(dǎo)數(shù)值。為未知點(diǎn)處的導(dǎo)數(shù)值。132021-11-15例例5.2 ：某液體冷卻時(shí)，溫度隨時(shí)間的變化數(shù)據(jù)如表某液體冷卻時(shí)，溫度隨時(shí)間的變化數(shù)據(jù)如表5-3所示：所示： 表表5-3 冷卻溫度隨時(shí)間的變化數(shù)據(jù)冷卻溫度隨時(shí)間的變化數(shù)據(jù)試分別計(jì)算試分別計(jì)算t2，3，4mi

15、n及及t1.5，2.5，4.5min時(shí)的降溫速率。時(shí)的降溫速率。解：解：三次樣條插值函數(shù)求數(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é)果：執(zhí)行結(jié)果：相應(yīng)時(shí)間時(shí)的降溫速率相應(yīng)時(shí)間時(shí)的降溫速率: 2.0000 3.0000 4.0000 1

16、.5000 2.5000 4.5000 -5.3722 -4.6722 -3.8389 -5.7972 -4.9889 -3.2222t (min)012345t(0c)92.085.379.574.570.267.0注：注：前者是計(jì)算節(jié)點(diǎn)處的一階導(dǎo)數(shù)，后者是計(jì)算非節(jié)點(diǎn)處的一階導(dǎo)數(shù)前者是計(jì)算節(jié)點(diǎn)處的一階導(dǎo)數(shù)，后者是計(jì)算非節(jié)點(diǎn)處的一階導(dǎo)數(shù)142021-11-155.1.3 最小二乘法樣條擬合函數(shù)求最小二乘法樣條擬合函數(shù)求微分微分n在實(shí)際化工應(yīng)用中，當(dāng)來(lái)自實(shí)驗(yàn)觀測(cè)的離在實(shí)際化工應(yīng)用中，當(dāng)來(lái)自實(shí)驗(yàn)觀測(cè)的離散數(shù)據(jù)不可避免地含有較大隨機(jī)誤差時(shí)，散數(shù)據(jù)不可避免地含有較大隨機(jī)誤差時(shí)，此時(shí)用插值公式求數(shù)值微分

17、雖然樣本點(diǎn)處此時(shí)用插值公式求數(shù)值微分雖然樣本點(diǎn)處誤差較小，但可能會(huì)使非樣本點(diǎn)處產(chǎn)生較誤差較小，但可能會(huì)使非樣本點(diǎn)處產(chǎn)生較大誤差大誤差n為此，可采用最小二乘法樣條擬合實(shí)驗(yàn)數(shù)為此，可采用最小二乘法樣條擬合實(shí)驗(yàn)數(shù)據(jù)，獲得一個(gè)函數(shù)模型，然后再對(duì)其求導(dǎo)據(jù)，獲得一個(gè)函數(shù)模型，然后再對(duì)其求導(dǎo)數(shù)。與樣條插值不同，樣條擬合不要求曲數(shù)。與樣條插值不同，樣條擬合不要求曲線經(jīng)過(guò)全部的數(shù)據(jù)點(diǎn)，這樣處理求導(dǎo)結(jié)果線經(jīng)過(guò)全部的數(shù)據(jù)點(diǎn)，這樣處理求導(dǎo)結(jié)果會(huì)有很大改善會(huì)有很大改善152021-11-15,iix y( )yf x2min( )iiiiew yf xiw1 1 方法概述方法概述 可用最小二乘法擬合成平滑可用最小二乘法

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

19、)其中：其中： x,y要處理的離散數(shù)據(jù)（）要處理的離散數(shù)據(jù)（） tol光滑時(shí)的允許精度，通常?。ü饣瑫r(shí)的允許精度，通常取（10-210-4）step 2:可用可用fnder函數(shù)求最小二乘平滑函數(shù)求最小二乘平滑b樣條擬合函數(shù)的導(dǎo)數(shù)樣條擬合函數(shù)的導(dǎo)數(shù);step3: 可用可用fnval函數(shù)求導(dǎo)函數(shù)在未知點(diǎn)處的導(dǎo)數(shù)值。函數(shù)求導(dǎo)函數(shù)在未知點(diǎn)處的導(dǎo)數(shù)值。 fnder( )和和fnval（）調(diào)用形式前文中已經(jīng)介紹過(guò)。（）調(diào)用形式前文中已經(jīng)介紹過(guò)。 2 2 最小二乘法平滑最小二乘法平滑b b樣條擬合函數(shù)求微分的樣條擬合函數(shù)求微分的matlabmatlab實(shí)現(xiàn)實(shí)現(xiàn)172021-11-15n由離散數(shù)據(jù)求數(shù)值微分的

20、四種方法及有關(guān)由離散數(shù)據(jù)求數(shù)值微分的四種方法及有關(guān)matlab函數(shù)：函數(shù)：（1）差分法）差分法 用差分函數(shù)用差分函數(shù)diff()近似計(jì)算導(dǎo)數(shù)，即近似計(jì)算導(dǎo)數(shù)，即dy=diff(y)./diff(x)。對(duì)于向量對(duì)于向量x，diff(x)表示了表示了x(2)-x(1) x(3)-x(2) . x(n)-x(n-1).對(duì)于矩陣對(duì)于矩陣x，diff(x)表示了表示了x(2:n,:) - x(1:n-1,:)小結(jié)小結(jié)注：注：此法用一階差分，精度較差，若改用二階差分，可大大提高精度，但編程麻煩些。此法用一階差分，精度較差，若改用二階差分，可大大提高精度，但編程麻煩些。()polyfit0(polyder(

21、)polyval（2）多項(xiàng)式擬合方法）多項(xiàng)式擬合方法向量向量p表示的多項(xiàng)式擬合函數(shù)表示的多項(xiàng)式擬合函數(shù)導(dǎo)函數(shù)導(dǎo)函數(shù)pppp在在xi的導(dǎo)數(shù)值。的導(dǎo)數(shù)值。其中：其中：函數(shù)函數(shù)polyfit( )和和polyval( )在前文中已介紹在前文中已介紹導(dǎo)函數(shù)導(dǎo)函數(shù)polyder( )的的調(diào)用格式調(diào)用格式為：為：pp=polyder(p) 離散數(shù)據(jù)離散數(shù)據(jù) 該函數(shù)對(duì)向量該函數(shù)對(duì)向量p表示的多項(xiàng)式函數(shù)進(jìn)行求導(dǎo)，返回導(dǎo)函數(shù)為表示的多項(xiàng)式函數(shù)進(jìn)行求導(dǎo)，返回導(dǎo)函數(shù)為pp。182021-11-15（3）三次樣條插值方法）三次樣條插值方法 （4）樣條擬合方法（最小二乘法）樣條擬合方法（最小二乘法） 其中：其中：函數(shù)函

22、數(shù)csaps()、spap2()、spaps()和和fnval()已在前已在前文中介紹，文中介紹， 求導(dǎo)函數(shù)求導(dǎo)函數(shù)fnder()的的調(diào)用格式調(diào)用格式為：為： pp=fnder(f,dorder) 該函數(shù)計(jì)算函數(shù)該函數(shù)計(jì)算函數(shù)f的的dorder階導(dǎo)函數(shù)，默認(rèn)階數(shù)階導(dǎo)函數(shù)，默認(rèn)階數(shù)dorder=1。 ()csapi()fnder()polyval離散數(shù)據(jù)離散數(shù)據(jù)三次樣條插值函數(shù)三次樣條插值函數(shù)cscs的導(dǎo)數(shù)的導(dǎo)數(shù)pppp在在xi的導(dǎo)數(shù)值的導(dǎo)數(shù)值()()2()spapsaspapcsaps或或()fnder ()fnval 離散數(shù)據(jù)離散數(shù)據(jù)樣條擬合函數(shù)樣條擬合函數(shù)spsp的導(dǎo)數(shù)的導(dǎo)數(shù)pppp在在x

23、i的導(dǎo)數(shù)值。的導(dǎo)數(shù)值。 192021-11-15 acacmaadck cdt例例5.3：反應(yīng)物反應(yīng)物a在一等溫間歇反應(yīng)器中發(fā)生的反應(yīng)為：在一等溫間歇反應(yīng)器中發(fā)生的反應(yīng)為：a測(cè)量得到的反應(yīng)器中不同時(shí)間下反應(yīng)物測(cè)量得到的反應(yīng)器中不同時(shí)間下反應(yīng)物a的濃度的濃度表表5-4 間歇反應(yīng)器動(dòng)力學(xué)數(shù)據(jù)間歇反應(yīng)器動(dòng)力學(xué)數(shù)據(jù)t(s)0204060120180300 (mol/l)10865321系統(tǒng)的動(dòng)力學(xué)模型為：系統(tǒng)的動(dòng)力學(xué)模型為：，試根據(jù)表中數(shù)據(jù)確定其反應(yīng)速率方程。，試根據(jù)表中數(shù)據(jù)確定其反應(yīng)速率方程。產(chǎn)物產(chǎn)物如表如表5-4所示。所示。202021-11-15ln()lnlnaadcmckdtln(),lna

24、adcyxcdtlnykmxadcdt解：解：（1）系統(tǒng)的動(dòng)力學(xué)模型為非線性形式，可將其線性化。）系統(tǒng)的動(dòng)力學(xué)模型為非線性形式，可將其線性化。對(duì)方程兩邊取對(duì)數(shù)：對(duì)方程兩邊取對(duì)數(shù)：令令 則原模型變?yōu)椋簞t原模型變?yōu)椋?（2）計(jì)算）計(jì)算t=0 20 40 60 120 180 300;ca=10 8 6 5 3 2 1;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(

25、sp,ti);plot(t,ca,bo,ti,cai,r-),xlabel(t),ylabel(ca)及得到擬合曲線的圖形程序如下：及得到擬合曲線的圖形程序如下：212021-11-15（3）線性擬合程序如下：）線性擬合程序如下：y=log(-dcadt);x=log(ca);p=polyfit(x,y,1);k=exp(p(2),m=p(1)執(zhí)行結(jié)果：執(zhí)行結(jié)果：k = 0.0059 m = 1.2904所以本例的反應(yīng)速率方程為：所以本例的反應(yīng)速率方程為：1.29040.0059aadccdt222021-11-151)被積函數(shù)以一組數(shù)據(jù)形式表示；被積函數(shù)以一組數(shù)據(jù)形式表示；2)被積函數(shù)過(guò)于特

26、殊或原函數(shù)不能用初等函數(shù)表示，積分表中無(wú)被積函數(shù)過(guò)于特殊或原函數(shù)不能用初等函數(shù)表示，積分表中無(wú) 法找到可沿用的現(xiàn)成公式；法找到可沿用的現(xiàn)成公式；3)有的原函數(shù)十分復(fù)雜難以計(jì)算。有的原函數(shù)十分復(fù)雜難以計(jì)算。badxxffi)()(對(duì)于積分：對(duì)于積分：公式有則由的原函數(shù)如果知道leibniznewtonxfxf),()(badxxf)()()()(afbfxfba但是在工程技術(shù)和科學(xué)研究中但是在工程技術(shù)和科學(xué)研究中,常會(huì)見(jiàn)到以下現(xiàn)象常會(huì)見(jiàn)到以下現(xiàn)象:0 導(dǎo)入導(dǎo)入5.2 數(shù)值積分?jǐn)?shù)值積分232021-11-15 數(shù)值積分?jǐn)?shù)值積分就是一種常用的近似計(jì)算方法。數(shù)就是一種常用的近似計(jì)算方法。數(shù)值積分方法不

27、受以上幾個(gè)問(wèn)題的限制，在化學(xué)化值積分方法不受以上幾個(gè)問(wèn)題的限制，在化學(xué)化工領(lǐng)域應(yīng)用甚廣，如反應(yīng)熱效應(yīng)計(jì)算、熱容計(jì)算、工領(lǐng)域應(yīng)用甚廣，如反應(yīng)熱效應(yīng)計(jì)算、熱容計(jì)算、熵的計(jì)算、反應(yīng)活化能的計(jì)算等。熵的計(jì)算、反應(yīng)活化能的計(jì)算等。 以上這些現(xiàn)象以上這些現(xiàn)象,newton-leibniz很難發(fā)揮作用，很難發(fā)揮作用，只能建立積分的近似計(jì)算方法只能建立積分的近似計(jì)算方法242021-11-15如：如：)(,h21tfcpdtcpntt252021-11-150.1 0.1 數(shù)值積分的基本思路和方法數(shù)值積分的基本思路和方法n 常用的數(shù)值積分的常用的數(shù)值積分的基本思路基本思路來(lái)自于插值法，它來(lái)自于插值法，它通過(guò)構(gòu)

28、造一個(gè)插值多項(xiàng)式通過(guò)構(gòu)造一個(gè)插值多項(xiàng)式 作為作為 的近似表的近似表達(dá)式，用達(dá)式，用 的積分值作為的積分值作為 的近似積分的近似積分值。值。n數(shù)值積分的數(shù)值積分的方法方法很豐富，常用的插值型求積公很豐富，常用的插值型求積公式有兩類(lèi)：一類(lèi)是等距節(jié)點(diǎn)的牛頓柯特斯求式有兩類(lèi)：一類(lèi)是等距節(jié)點(diǎn)的牛頓柯特斯求積公式；積公式； 另一類(lèi)是不等距節(jié)點(diǎn)的高斯型求積另一類(lèi)是不等距節(jié)點(diǎn)的高斯型求積公式。公式。( )np x( )f x( )np x( )f x262021-11-15 newton-cotes公式是指公式是指等距節(jié)點(diǎn)下等距節(jié)點(diǎn)下使用使用lagrange插值插值 多項(xiàng)式多項(xiàng)式建立的數(shù)值求積公式建立的數(shù)值求

29、積公式,)(bacxf設(shè)函數(shù)等份分割為將積分區(qū)間nba,nkkhaxk, 1 ,0,為步長(zhǎng)其中nabh各節(jié)點(diǎn)為各節(jié)點(diǎn)為插值多項(xiàng)式為的lagrangexf)(則：則：5.2.1牛頓柯特斯（牛頓柯特斯（newton-cotes）求積公式）求積公式1 1 方法概述方法概述272021-11-15這里這里 是既不依賴(lài)于被積函數(shù)，也不依賴(lài)于積分區(qū)間的常數(shù)，是既不依賴(lài)于被積函數(shù)，也不依賴(lài)于積分區(qū)間的常數(shù)，稱(chēng)為柯特斯系數(shù)。式（稱(chēng)為柯特斯系數(shù)。式（5-3）稱(chēng)為牛頓柯特斯求積公式。）稱(chēng)為牛頓柯特斯求積公式。00( )() ( )()nnbbkkkkaakkif x dxf x lx dxa f x0( )bbj

30、kkaajnkjjkxxalx dxdxxxxath/()iicaba00( 1)()! ()!nknijnjkctj dtn knkic其中：其中： 一般地：令一般地：令 ，得：得： （5-3）282021-11-15在在newton-cotes公式中公式中,n=1,2,4時(shí)的公式是最常用也是最重時(shí)的公式是最常用也是最重要三個(gè)公式要三個(gè)公式,稱(chēng)為低階公式稱(chēng)為低階公式（當(dāng)（當(dāng)n=0時(shí)的公式為矩形公式）時(shí)的公式為矩形公式）1）梯形）梯形(trapezia)公式公式abhbxaxn, 110則取dtt10)1()1(0ccotes系數(shù)為系數(shù)為21dtt10)1(1c21求積公式為求積公式為2920

31、21-11-151(1)1010( )()() ()()2kkkbaifbacf xf xf x)()(2bfafab)(1fi即上式稱(chēng)為上式稱(chēng)為梯形求積公式梯形求積公式,（5-4）302021-11-152）simpson公式公式2,2,2210abhbxabxaxn則取cotes系數(shù)為系數(shù)為dtttc20)2(0)2)(1(4161dtttc20)2(1)2(2164dtttc20)2(2)1(4161312021-11-15)(61)(64)(61)(210 xfxfxfab)()2(4)(6bfbafafab)(2fi上式稱(chēng)為上式稱(chēng)為simpson求積公式求積公式，也稱(chēng)，也稱(chēng)三點(diǎn)公式或

32、拋物線公式三點(diǎn)公式或拋物線公式求積公式為求積公式為2i20)2()()(kkkxfcab式（式（5-4）和式（）和式（5-5）是化工領(lǐng)域常用的兩個(gè)求積公式。與梯形法求）是化工領(lǐng)域常用的兩個(gè)求積公式。與梯形法求積公式相比，積公式相比，simpson法求積公式是一個(gè)較高精度的求積公式。用法求積公式是一個(gè)較高精度的求積公式。用式（式（5-3）還可得到更高階數(shù)）還可得到更高階數(shù)（5-5）322021-11-15n newton-cotes公式當(dāng)公式當(dāng)n大于大于7時(shí)，公式的穩(wěn)定性時(shí)，公式的穩(wěn)定性將無(wú)法保證，因此將無(wú)法保證，因此,在實(shí)際應(yīng)用中一般不使用高階在實(shí)際應(yīng)用中一般不使用高階newton-cotes

33、公式而是采用公式而是采用低階復(fù)合求積法低階復(fù)合求積法n在實(shí)際計(jì)算中為了保證計(jì)算的精度，往往首先用在實(shí)際計(jì)算中為了保證計(jì)算的精度，往往首先用分點(diǎn)分點(diǎn)xk=a+kh, (k=0,1,n)將區(qū)間將區(qū)間a,b分成分成n個(gè)個(gè)相等的子區(qū)間，而后對(duì)每個(gè)子區(qū)相等的子區(qū)間，而后對(duì)每個(gè)子區(qū) 間再應(yīng)用梯形公間再應(yīng)用梯形公式或式或simpson公式，分別得到：公式，分別得到： nabh11)(2)()(2 復(fù)化梯形公式nknkhafbfafht2 2 復(fù)化法求積公式復(fù)化法求積公式11102 ()4()()6nnkkkkhsf xf xf x復(fù)化復(fù)化simpson公式：公式：332021-11-15 盡管復(fù)化法求積公式

34、具有很高的精度，但是它必須采用等步長(zhǎng)方法，盡管復(fù)化法求積公式具有很高的精度，但是它必須采用等步長(zhǎng)方法，從而限制了它的效率。這里我們介紹一種更加靈活選取步長(zhǎng)的方法，即從而限制了它的效率。這里我們介紹一種更加靈活選取步長(zhǎng)的方法，即自適應(yīng)步長(zhǎng)法。自適應(yīng)步長(zhǎng)法。（5-8）,kka bkkkhba11 ()4 ()()62kkkkkhsf af ahf b2113 ()4 ()2 ()4 ()()12424kkkkkkkkkhsf af ahf ahf ahf b210.1()ss 210.1()ss （5-6） （5-7），令，令，當(dāng)，當(dāng) ,kka b2s上上simpson積分積分 達(dá)到精度達(dá)到精度時(shí)

35、，可認(rèn)為區(qū)間時(shí)，可認(rèn)為區(qū)間取計(jì)算需滿(mǎn)足的精度為取計(jì)算需滿(mǎn)足的精度為 以以simpson積分法為例，某區(qū)間積分法為例，某區(qū)間 ，記，記 ，考慮該區(qū)間，考慮該區(qū)間上的上的simpson積分和二等分以后的兩個(gè)積分和二等分以后的兩個(gè)simpson積分和：積分和：3 3 自適應(yīng)求積公式自適應(yīng)求積公式342021-11-15()/2kkhba/2k 這樣重復(fù)下去，直至每個(gè)分段部分達(dá)到相應(yīng)精度（步長(zhǎng)為這樣重復(fù)下去，直至每個(gè)分段部分達(dá)到相應(yīng)精度（步長(zhǎng)為時(shí)精度為時(shí)精度為 ）。這樣，不同段的步長(zhǎng)可能是不一樣的，積分結(jié)果為）。這樣，不同段的步長(zhǎng)可能是不一樣的，積分結(jié)果為每一小段的積分總和。每一小段的積分總和。,a

36、b2s/2 自適應(yīng)步長(zhǎng)自適應(yīng)步長(zhǎng)simpson法從法從 開(kāi)始按式（開(kāi)始按式（5-6）（）（5-8）的方法檢驗(yàn)）的方法檢驗(yàn)。若滿(mǎn)足精度，則以。若滿(mǎn)足精度，則以 為計(jì)算結(jié)果；若不滿(mǎn)足精度，則分成兩個(gè)小區(qū)為計(jì)算結(jié)果；若不滿(mǎn)足精度，則分成兩個(gè)小區(qū)間間各自重復(fù)逐步上述過(guò)程，每個(gè)小區(qū)間精度為各自重復(fù)逐步上述過(guò)程，每個(gè)小區(qū)間精度為352021-11-154 newton-cotes4 newton-cotes求積公式的求積公式的matlabmatlab實(shí)現(xiàn)實(shí)現(xiàn)n常用的三種常用的三種newton-cotes系列數(shù)值積分法的相應(yīng)系列數(shù)值積分法的相應(yīng)matlab函函數(shù)如下數(shù)如下:n1) 復(fù)合梯形法數(shù)值積分：復(fù)合梯

37、形法數(shù)值積分：trapz( )n調(diào)用形式：調(diào)用形式：z = trapz(x,y)n其中其中: x,y分別代表數(shù)目相同的向量或數(shù)組，而分別代表數(shù)目相同的向量或數(shù)組，而y 與與x 的的 關(guān)系可以是一關(guān)系可以是一 函數(shù)型態(tài)（如函數(shù)型態(tài)（如y=sin(x)）或是不以函數(shù)描）或是不以函數(shù)描述的離散型態(tài)；述的離散型態(tài)； z代表返回的積分值；代表返回的積分值；注注:離散型態(tài)數(shù)據(jù)用離散型態(tài)數(shù)據(jù)用trapz函數(shù)時(shí)，還需設(shè)定函數(shù)時(shí)，還需設(shè)定x在區(qū)間在區(qū)間 a,b 之間離之間離散點(diǎn)的間隔；缺省參數(shù)散點(diǎn)的間隔；缺省參數(shù)x時(shí)，表示時(shí)，表示x被等分，每份寬為被等分，每份寬為1。362021-11-152） 自適應(yīng)自適應(yīng)s

38、impson法數(shù)值積分：法數(shù)值積分：quad（）（）基本調(diào)用格式基本調(diào)用格式： q=quad(fun,a,b) 或或 q=quad(fun,a,b,tol,trace,p1,p2,)其中其中: fun被積函數(shù)?？梢允莾?nèi)置函數(shù)、被積函數(shù)?？梢允莾?nèi)置函數(shù)、m文件或函數(shù)句柄文件或函數(shù)句柄, 函數(shù)表達(dá)式函數(shù)表達(dá)式 中的必須使用點(diǎn)運(yùn)算符號(hào)。中的必須使用點(diǎn)運(yùn)算符號(hào)。 a, b分別是積分的下限和上限；分別是積分的下限和上限； q積分結(jié)果。積分結(jié)果。 tol默認(rèn)誤差限，默認(rèn)值為默認(rèn)誤差限，默認(rèn)值為1.e-6. trace-取取0表示不用圖形顯示積分過(guò)程，非表示不用圖形顯示積分過(guò)程，非0表示用圖形表示用圖形顯示

39、積分過(guò)程顯示積分過(guò)程 p1,p2,直接傳遞給函數(shù)直接傳遞給函數(shù)fun的參數(shù)。的參數(shù)。3）自適應(yīng)）自適應(yīng)lobatto法數(shù)值積分：法數(shù)值積分：quadl（）（） quadl是高階的自適應(yīng)是高階的自適應(yīng)newton-cotes數(shù)值積分法函數(shù)，它數(shù)值積分法函數(shù)，它比比quad函數(shù)更有效，精度更高。其使用方法和函數(shù)更有效，精度更高。其使用方法和quad（）完（）完全相同。全相同。需要了解更多的內(nèi)容可查閱需要了解更多的內(nèi)容可查閱matlab功能函數(shù)庫(kù)（功能函數(shù)庫(kù)（funfun）。）。372021-11-15flglg2.303afprt例例5.4：真實(shí)氣體的逸度真實(shí)氣體的逸度可用下式計(jì)算：可用下式計(jì)算：

40、 現(xiàn)測(cè)得現(xiàn)測(cè)得00c下氫氣的有關(guān)數(shù)值如表下氫氣的有關(guān)數(shù)值如表5-5所示，試求所示，試求1000 atm下的逸度。下的逸度。6310vmrtvpp6310vm表表5-5 00c下氫氣的相關(guān)數(shù)值下氫氣的相關(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-真實(shí)氣體的實(shí)測(cè)體積和按理想氣體真實(shí)氣體的實(shí)測(cè)體積和按理想氣體定律計(jì)算得到的體積之間的差值。定

41、律計(jì)算得到的體積之間的差值。r-氣體常數(shù)氣體常數(shù) 0,prtadpvpfp6382.06 10/()matm mol k其中：其中： -逸度；逸度； -壓力，壓力，atm；t-絕對(duì)溫度，絕對(duì)溫度，k。 r tvpfrtvpp382021-11-15解：解：本題是離散型數(shù)據(jù)，可用本題是離散型數(shù)據(jù)，可用trapz函數(shù)求解數(shù)值積分。函數(shù)求解數(shù)值積分。（1）計(jì)算數(shù)值積分，程序如下：）計(jì)算數(shù)值積分，程序如下：%計(jì)算數(shù)值積分p=0:100:1000;a1=15.46 15.46 15.46 15.61 15.85 15.93 16.09 16.13 16.16 16.16 16.14;a1=-a1;a=t

42、rapz(p,a1);% 梯形法計(jì)算數(shù)值積分a=-a執(zhí)行結(jié)果：執(zhí)行結(jié)果：a = 15865（2）計(jì)算逸度值，程序如下：）計(jì)算逸度值，程序如下：%計(jì)算逸度lf=log10(1000)+a./(2.303.*82.06.*273.2);% lf代表f=10.lf執(zhí)行結(jié)果：執(zhí)行結(jié)果：f = 2.0290e+003所以所以 f =2029.00 atmf392021-11-150.4,0.9fdxx5r () /dfxxddxnyxxyrxy例例5.5：氯仿：氯仿-苯雙組分精餾系統(tǒng)的氣液平衡數(shù)據(jù)如表苯雙組分精餾系統(tǒng)的氣液平衡數(shù)據(jù)如表5-6所示。規(guī)定所示。規(guī)定進(jìn)料和塔頂?shù)慕M成分別是進(jìn)料和塔頂?shù)慕M成分別是

43、 ，精餾段的回流比為，精餾段的回流比為精餾段理論板數(shù)的模型為精餾段理論板數(shù)的模型為表表5-6 精餾段氣液平衡數(shù)據(jù)精餾段氣液平衡數(shù)據(jù)0.1780.2750.3720.4560.6500.8440.2430.3820.5180.6160.7950.931，試用試用matlab計(jì)算所需的精餾段理論板數(shù)。計(jì)算所需的精餾段理論板數(shù)。402021-11-15n解：解：因模型中的的函數(shù)關(guān)系是以表格形式給出，我們?nèi)粢蚰Ｐ椭械牡暮瘮?shù)關(guān)系是以表格形式給出，我們?nèi)粢眯疗丈ǖ扔?jì)算較精確的精餾段理論板數(shù)的時(shí)候，需先用辛普森法等計(jì)算較精確的精餾段理論板數(shù)的時(shí)候，需先采用擬合法將離散數(shù)據(jù)（采用擬合法將離散數(shù)據(jù)（ ，

44、）擬合成多項(xiàng)式，再將多）擬合成多項(xiàng)式，再將多項(xiàng)式代入被積函數(shù)求積。項(xiàng)式代入被積函數(shù)求積。計(jì)算程序如下：function li55clear all;xi=0.178 0.275 0.372 0.456 0.650 0.844;yi=0.243 0.382 0.518 0.616 0.795 0.931;sp=spline(xi,yi);% 用spline（）擬合多項(xiàng)式% 畫(huà)出擬合曲線，直觀比較擬合效果xplot=linspace(xi(1),xi(end),100);yplot=fnval(sp,xplot);plot(xi,yi,o,xplot,yplot,-);n=quad(func1,0

45、.4,0.9, , ,sp); % 計(jì)算精餾段理論板數(shù)n=round(n+0.5) % 數(shù)據(jù)取整function f=func1(x,sp)y=fnval(sp,x);f=1./(y-x-(0.9-y)./5);執(zhí)行結(jié)果：執(zhí)行結(jié)果：n = 5所以，所以，精餾段理論板數(shù)為精餾段理論板數(shù)為5塊。塊。ixiy412021-11-15 5.2.2 高斯勒讓德公式高斯勒讓德公式高斯法是一種精度較高的求積分法，它的高斯法是一種精度較高的求積分法，它的一般公式是一般公式是式中式中(x)是一個(gè)權(quán)重函數(shù)，是一個(gè)權(quán)重函數(shù)，aj為系數(shù)，為系數(shù)，xj為橫坐標(biāo)上的節(jié)點(diǎn)為橫坐標(biāo)上的節(jié)點(diǎn) 高斯勒讓德多項(xiàng)式的權(quán)重函數(shù)高斯勒讓德多項(xiàng)式的權(quán)重函數(shù)(x)=1，因而，一個(gè)，因而，一個(gè)n點(diǎn)高斯勒讓德求積公點(diǎn)高斯勒讓德求積公式具有如下形式：式具有如下形式：banjjjxfadxxfx1)()()(111( )()njjjf x dxa f x1 方法概述方法概述 422021-11-15 右邊的右邊的f(xj)是函數(shù)是函數(shù)f(x)在節(jié)點(diǎn)在節(jié)點(diǎn)xj處的值，處的值， 節(jié)點(diǎn)節(jié)點(diǎn)xj是