05数值微分与数值积分

1、2021-11-151第第5章章 数值微分与数值积分数值微分与数值积分22021-11-15n微积分是高等数学中的重要内容，在化学微积分是高等数学中的重要内容，在化学工程上有许多非常重要的应用工程上有许多非常重要的应用n微积分的数值方法，不同于高等数学中的微积分的数值方法，不同于高等数学中的解析方法，尤其适合求解没有或很难求出解析方法，尤其适合求解没有或很难求出微分或积分表达式的实际化工问题的计算，微分或积分表达式的实际化工问题的计算，例如：列表函数求微分或积分例如：列表函数求微分或积分 引言引言1-数值微积分方法不同于解析方法数值微积分方法不同于解析方法32021-11-15n数值微分和数值

2、积分与插值和拟合往往是密数值微分和数值积分与插值和拟合往往是密不可分的不可分的n如在进行数值微分时，常针对离散的数据点，如在进行数值微分时，常针对离散的数据点，利用插值和拟合可以减少误差；而数值积分利用插值和拟合可以减少误差；而数值积分的基本思路也来自于插值法。例如如果所积的基本思路也来自于插值法。例如如果所积函数的形式比较复杂或以表格形式给出，则函数的形式比较复杂或以表格形式给出，则可通过构造一个插值多项式来代替原函数，可通过构造一个插值多项式来代替原函数，从而使问题简化从而使问题简化引言引言2-数值微分和数值积分与插值和拟合数值微分和数值积分与插值和拟合关系密切关系密切42021-11-1

3、55.1 数值微分数值微分n化工领域的实际问题中时常需要求列表函数在节化工领域的实际问题中时常需要求列表函数在节点和非节点处的导数值，这是数值微分所要解决点和非节点处的导数值，这是数值微分所要解决的问题。数值微分方法可近似求出某点的导数值的问题。数值微分方法可近似求出某点的导数值n例如在反应动力学的研究中，根据实验数据确定反应的动例如在反应动力学的研究中，根据实验数据确定反应的动力学方程力学方程 ：n这里实验测得一批离散点，要计算这里实验测得一批离散点，要计算 只能借助数值只能借助数值微分求导解决微分求导解决 表表5-1 反应动力学实验数据反应动力学实验数据 mapadpk pdtadpdtt

4、1t2tnpa1pa2pan0 0 导入导入52021-11-150.1 0.1 建立数值微分公式的三种思路建立数值微分公式的三种思路 常用三种思路建立数值微分公式：常用三种思路建立数值微分公式：1.从微分定义出发，通过近似处理，得到数值微从微分定义出发，通过近似处理，得到数值微分的近似公式分的近似公式2.从插值近似公式出发，对插值公式的近似求导从插值近似公式出发，对插值公式的近似求导可得到数值微分的近似公式可得到数值微分的近似公式3.先用最小二乘拟合方法根据已知数据得近似函先用最小二乘拟合方法根据已知数据得近似函数，再对此近似函数求微分可得到数值微分的数，再对此近似函数求微分可得到数值微分的

5、近似公式。然后对各方法数值微分后得到的多近似公式。然后对各方法数值微分后得到的多项式求值，即可求出任意点处的任意阶微分项式求值，即可求出任意点处的任意阶微分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 方法概述方法概述在微积分中，一阶微分的计算可以在二相

6、邻点x+h和x间函数取下列极限求得： 取其达到极限前的形式取其达到极限前的形式，就得到以下微分的差分近似式：，就得到以下微分的差分近似式： 注：注：高阶微分项可以利用低阶微分项来计算，如二阶微分式可以表示为：高阶微分项可以利用低阶微分项来计算，如二阶微分式可以表示为： 对应的差分式有：对应的差分式有：5.1.1 差分近似微分差分近似微分上式中三种不同表示形式依次是一阶前向差分、一阶后向差分和一阶中上式中三种不同表示形式依次是一阶前向差分、一阶后向差分和一阶中心差分来近似表示微分。其中一阶中心差分的精度较高。心差分来近似表示微分。其中一阶中心差分的精度较高。72021-11-152 2 差分的差

7、分的matlabmatlab实现实现n在在matlab中，可用中，可用diff函数进行离散数据的近似函数进行离散数据的近似求导求导n调用形式调用形式：y = diff(x,n)n其中：其中：x表示求导变量，可以是向量或矩阵。如表示求导变量，可以是向量或矩阵。如是矩阵形式则按各列作差分；是矩阵形式则按各列作差分； n表示表示n阶差分，即差分阶差分，即差分n次次 y是是x的差分结果的差分结果n注：注：用用diff函数进行离散数据的近似求导与前向差函数进行离散数据的近似求导与前向差分近似，分近似， 但误差较大。最好将数据利用插值或拟但误差较大。最好将数据利用插值或拟合得到多项式，然后对近似多项式进行

8、微分合得到多项式，然后对近似多项式进行微分82021-11-15464622()c hc h apadpdt例例5.1：丁二烯的气相二聚反应如下：丁二烯的气相二聚反应如下：实验在一定容器的反应器中进行，实验在一定容器的反应器中进行，3260c时，测得物系中丁二烯的分时，测得物系中丁二烯的分压压 (mmhg)与时间的关系如表与时间的关系如表5-2所示。所示。apap表表 5-2 丁二烯二聚反应实验数据丁二烯二聚反应实验数据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用数值微分法计算所列时刻每一瞬间的反应速率用数值微分法计算所列时刻每一瞬间的反应速率 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); 求时间t的差分dpa=diff(pa); 求压力的差分q=

10、dpa./dt q为数值微分结果执行结果：执行结果：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 三次样条插值函数求微分三次样条插值函数求微分( )s x( )f x( )s x( )fx(

11、 )s x( )f x若三次样条插值函数若三次样条插值函数收敛于收敛于，那么导数，那么导数收敛于收敛于，因此用样条插值函数，因此用样条插值函数作为作为函数，不但彼此的函数值非常接近，而且导数值也很接近。函数，不但彼此的函数值非常接近，而且导数值也很接近。 用三次样条插值函数求数值导数是可靠的，这是化工计用三次样条插值函数求数值导数是可靠的，这是化工计算中求数值微分的有效方法。算中求数值微分的有效方法。的近似的近似112021-11-15( )( )fxs x1,2,in1,iixxx11()()( )( )iiiiiixxxxfxs xmmhh221111()()( )( )()226iiii

12、iiiiiiiiyyxxxxhfxs xmmmmhhh 1 1 方法概述方法概述用三次样条插值函数建立的数值微分公式为：用三次样条插值函数建立的数值微分公式为： 求导得：求导得： （5-2）;式（式（5-1）和（）和（5-2）不但适用于求节点处的导数，而且可求非节点处的导数。）不但适用于求节点处的导数，而且可求非节点处的导数。（5-1）其中，其中，122021-11-152 2 三次样条插值函数求微分的三次样条插值函数求微分的matlabmatlab函数函数matlab求离散数据的三次样条插值函数微分方法分三个步骤：求离散数据的三次样条插值函数微分方法分三个步骤：step 1:对离散数据用对离

13、散数据用csapi函数（或函数（或spline函数）得到其三次样条插值函数函数）得到其三次样条插值函数调用形式调用形式 pp = csapi(x,y)其中：其中： x,y分别为离散数据对的自变量和因变量；分别为离散数据对的自变量和因变量； pp为得到的三次样条插值函数。为得到的三次样条插值函数。step 2: 用用fnder函数求三次样条插值函数的导数函数求三次样条插值函数的导数调用形式调用形式 fprime = fnder(f,dorder)其中：其中： f为三次样条插值函数；为三次样条插值函数； dorder为三次样条插值函数的求导阶数；为三次样条插值函数的求导阶数； fprime为得到的

14、三次样条插值函数的导函数。为得到的三次样条插值函数的导函数。step3:可用可用fnval函数求导函数在未知点处的导数值函数求导函数在未知点处的导数值调用形式调用形式 v = fnval(fprime,x)其中：其中：fprime为三次样条插值函数导函数；为三次样条插值函数导函数； x为未知点处自变量值；为未知点处自变量值； v为未知点处的导数值。为未知点处的导数值。132021-11-15例例5.2 ：某液体冷却时，温度随时间的变化数据如表某液体冷却时，温度随时间的变化数据如表5-3所示：所示： 表表5-3 冷却温度随时间的变化数据冷却温度随时间的变化数据试分别计算试分别计算t2，3，4mi

15、n及及t1.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.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注：注：前者是计算节点处的一阶导数，后者是计算非节点处的一阶导数前者是计算节点处的一阶导数，后者是计算非节点处的一阶导数142021-11-155.1.3 最小二乘法样条拟合函数求最小二乘法样条拟合函数求微分微分n在实际化工应用中，当来自实验观测的离在实际化工应用中，当来自实验观测的离散数据不可避免地含有较大随机误差时，散数据不可避免地含有较大随机误差时，此时用插值公式求数值微分

17、虽然样本点处此时用插值公式求数值微分虽然样本点处误差较小，但可能会使非样本点处产生较误差较小，但可能会使非样本点处产生较大误差大误差n为此，可采用最小二乘法样条拟合实验数为此，可采用最小二乘法样条拟合实验数据，获得一个函数模型，然后再对其求导据，获得一个函数模型，然后再对其求导数。与样条插值不同，样条拟合不要求曲数。与样条插值不同，样条拟合不要求曲线经过全部的数据点，这样处理求导结果线经过全部的数据点，这样处理求导结果会有很大改善会有很大改善152021-11-15,iix y( )yf x2min( )iiiiew yf xiw1 1 方法概述方法概述 可用最小二乘法拟合成平滑可用最小二乘法

18、拟合成平滑b样条曲线，即对于离散数据（样条曲线，即对于离散数据（所求的所求的k次样条拟合函数次样条拟合函数满足：满足：其中：其中：再对拟合函数作平滑处理后求导，即可求出任意点处微分。再对拟合函数作平滑处理后求导，即可求出任意点处微分。）为权重系数，默认为为权重系数，默认为1162021-11-15nmatlab求离散数据的最小二乘法平滑求离散数据的最小二乘法平滑b样条拟合函数求微分样条拟合函数求微分共三个步骤：共三个步骤：step 1:对离散数据用对离散数据用spaps函数得到最小二乘平滑函数得到最小二乘平滑b样条拟合样条拟合函数。函数。调用格式：调用格式：sp = spaps(x,y,tol

19、)其中：其中： x,y要处理的离散数据（）要处理的离散数据（） tol光滑时的允许精度，通常取（光滑时的允许精度，通常取（10-210-4）step 2:可用可用fnder函数求最小二乘平滑函数求最小二乘平滑b样条拟合函数的导数样条拟合函数的导数;step3: 可用可用fnval函数求导函数在未知点处的导数值。函数求导函数在未知点处的导数值。 fnder( )和和fnval（）调用形式前文中已经介绍过。（）调用形式前文中已经介绍过。 2 2 最小二乘法平滑最小二乘法平滑b b样条拟合函数求微分的样条拟合函数求微分的matlabmatlab实现实现172021-11-15n由离散数据求数值微分的

20、四种方法及有关由离散数据求数值微分的四种方法及有关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,:)小结小结注：注：此法用一阶差分，精度较差，若改用二阶差分，可大大提高精度，但编程麻烦些。此法用一阶差分，精度较差，若改用二阶差分，可大大提高精度，但编程麻烦些。()polyfit0(polyder(

21、)polyval（2）多项式拟合方法）多项式拟合方法向量向量p表示的多项式拟合函数表示的多项式拟合函数导函数导函数pppp在在xi的导数值。的导数值。其中：其中：函数函数polyfit( )和和polyval( )在前文中已介绍在前文中已介绍导函数导函数polyder( )的的调用格式调用格式为：为：pp=polyder(p) 离散数据离散数据 该函数对向量该函数对向量p表示的多项式函数进行求导，返回导函数为表示的多项式函数进行求导，返回导函数为pp。182021-11-15（3）三次样条插值方法）三次样条插值方法 （4）样条拟合方法（最小二乘法）样条拟合方法（最小二乘法） 其中：其中：函数函

22、数csaps()、spap2()、spaps()和和fnval()已在前已在前文中介绍，文中介绍， 求导函数求导函数fnder()的的调用格式调用格式为：为： pp=fnder(f,dorder) 该函数计算函数该函数计算函数f的的dorder阶导函数，默认阶数阶导函数，默认阶数dorder=1。 ()csapi()fnder()polyval离散数据离散数据三次样条插值函数三次样条插值函数cscs的导数的导数pppp在在xi的导数值的导数值()()2()spapsaspapcsaps或或()fnder ()fnval 离散数据离散数据样条拟合函数样条拟合函数spsp的导数的导数pppp在在x

23、i的导数值。的导数值。 192021-11-15 acacmaadck cdt例例5.3：反应物反应物a在一等温间歇反应器中发生的反应为：在一等温间歇反应器中发生的反应为：a测量得到的反应器中不同时间下反应物测量得到的反应器中不同时间下反应物a的浓度的浓度表表5-4 间歇反应器动力学数据间歇反应器动力学数据t(s)0204060120180300 (mol/l)10865321系统的动力学模型为：系统的动力学模型为：，试根据表中数据确定其反应速率方程。，试根据表中数据确定其反应速率方程。产物产物如表如表5-4所示。所示。202021-11-15ln()lnlnaadcmckdtln(),lna

24、adcyxcdtlnykmxadcdt解：解：（1）系统的动力学模型为非线性形式，可将其线性化。）系统的动力学模型为非线性形式，可将其线性化。对方程两边取对数：对方程两边取对数：令令 则原模型变为：则原模型变为： （2）计算）计算t=0 20 40 60 120 180 300;ca=10 8 6 5 3 2 1;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(

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)执行结果：执行结果：k = 0.0059 m = 1.2904所以本例的反应速率方程为：所以本例的反应速率方程为：1.29040.0059aadccdt222021-11-151)被积函数以一组数据形式表示；被积函数以一组数据形式表示；2)被积函数过于特

26、殊或原函数不能用初等函数表示，积分表中无被积函数过于特殊或原函数不能用初等函数表示，积分表中无 法找到可沿用的现成公式；法找到可沿用的现成公式；3)有的原函数十分复杂难以计算。有的原函数十分复杂难以计算。badxxffi)()(对于积分：对于积分：公式有则由的原函数如果知道leibniznewtonxfxf),()(badxxf)()()()(afbfxfba但是在工程技术和科学研究中但是在工程技术和科学研究中,常会见到以下现象常会见到以下现象:0 导入导入5.2 数值积分数值积分232021-11-15 数值积分数值积分就是一种常用的近似计算方法。数就是一种常用的近似计算方法。数值积分方法不

27、受以上几个问题的限制，在化学化值积分方法不受以上几个问题的限制，在化学化工领域应用甚广，如反应热效应计算、热容计算、工领域应用甚广，如反应热效应计算、热容计算、熵的计算、反应活化能的计算等。熵的计算、反应活化能的计算等。 以上这些现象以上这些现象,newton-leibniz很难发挥作用，很难发挥作用，只能建立积分的近似计算方法只能建立积分的近似计算方法242021-11-15如：如：)(,h21tfcpdtcpntt252021-11-150.1 0.1 数值积分的基本思路和方法数值积分的基本思路和方法n 常用的数值积分的常用的数值积分的基本思路基本思路来自于插值法，它来自于插值法，它通过构

28、造一个插值多项式通过构造一个插值多项式 作为作为 的近似表的近似表达式，用达式，用 的积分值作为的积分值作为 的近似积分的近似积分值。值。n数值积分的数值积分的方法方法很丰富，常用的插值型求积公很丰富，常用的插值型求积公式有两类：一类是等距节点的牛顿柯特斯求式有两类：一类是等距节点的牛顿柯特斯求积公式；积公式； 另一类是不等距节点的高斯型求积另一类是不等距节点的高斯型求积公式。公式。( )np x( )f x( )np x( )f x262021-11-15 newton-cotes公式是指公式是指等距节点下等距节点下使用使用lagrange插值插值 多项式多项式建立的数值求积公式建立的数值求

29、积公式,)(bacxf设函数等份分割为将积分区间nba,nkkhaxk, 1 ,0,为步长其中nabh各节点为各节点为插值多项式为的lagrangexf)(则：则：5.2.1牛顿柯特斯（牛顿柯特斯（newton-cotes）求积公式）求积公式1 1 方法概述方法概述272021-11-15这里这里 是既不依赖于被积函数，也不依赖于积分区间的常数，是既不依赖于被积函数，也不依赖于积分区间的常数，称为柯特斯系数。式（称为柯特斯系数。式（5-3）称为牛顿柯特斯求积公式。）称为牛顿柯特斯求积公式。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时的公式是最常用也是最重时的公式是最常用也是最重要三个公式要三个公式,称为低阶公式称为低阶公式（当（当n=0时的公式为矩形公式）时的公式为矩形公式）1）梯形）梯形(trapezia)公式公式abhbxaxn, 110则取dtt10)1()1(0ccotes系数为系数为21dtt10)1(1c21求积公式为求积公式为2920

31、21-11-151(1)1010( )()() ()()2kkkbaifbacf xf xf x)()(2bfafab)(1fi即上式称为上式称为梯形求积公式梯形求积公式,（5-4）302021-11-152）simpson公式公式2,2,2210abhbxabxaxn则取cotes系数为系数为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上式称为上式称为simpson求积公式求积公式，也称，也称三点公式或

32、抛物线公式三点公式或抛物线公式求积公式为求积公式为2i20)2()()(kkkxfcab式（式（5-4）和式（）和式（5-5）是化工领域常用的两个求积公式。与梯形法求）是化工领域常用的两个求积公式。与梯形法求积公式相比，积公式相比，simpson法求积公式是一个较高精度的求积公式。用法求积公式是一个较高精度的求积公式。用式（式（5-3）还可得到更高阶数）还可得到更高阶数（5-5）322021-11-15n newton-cotes公式当公式当n大于大于7时，公式的稳定性时，公式的稳定性将无法保证，因此将无法保证，因此,在实际应用中一般不使用高阶在实际应用中一般不使用高阶newton-cotes

33、公式而是采用公式而是采用低阶复合求积法低阶复合求积法n在实际计算中为了保证计算的精度，往往首先用在实际计算中为了保证计算的精度，往往首先用分点分点xk=a+kh, (k=0,1,n)将区间将区间a,b分成分成n个个相等的子区间，而后对每个子区相等的子区间，而后对每个子区 间再应用梯形公间再应用梯形公式或式或simpson公式，分别得到：公式，分别得到： nabh11)(2)()(2 复化梯形公式nknkhafbfafht2 2 复化法求积公式复化法求积公式11102 ()4()()6nnkkkkhsf xf xf x复化复化simpson公式：公式：332021-11-15 尽管复化法求积公式

34、具有很高的精度，但是它必须采用等步长方法，尽管复化法求积公式具有很高的精度，但是它必须采用等步长方法，从而限制了它的效率。这里我们介绍一种更加灵活选取步长的方法，即从而限制了它的效率。这里我们介绍一种更加灵活选取步长的方法，即自适应步长法。自适应步长法。（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），令，令，当，当 ,kka b2s上上simpson积分积分 达到精度达到精度时

35、，可认为区间时，可认为区间取计算需满足的精度为取计算需满足的精度为 以以simpson积分法为例，某区间积分法为例，某区间 ，记，记 ，考虑该区间，考虑该区间上的上的simpson积分和二等分以后的两个积分和二等分以后的两个simpson积分和：积分和：3 3 自适应求积公式自适应求积公式342021-11-15()/2kkhba/2k 这样重复下去，直至每个分段部分达到相应精度（步长为这样重复下去，直至每个分段部分达到相应精度（步长为时精度为时精度为 ）。这样，不同段的步长可能是不一样的，积分结果为）。这样，不同段的步长可能是不一样的，积分结果为每一小段的积分总和。每一小段的积分总和。,a

36、b2s/2 自适应步长自适应步长simpson法从法从 开始按式（开始按式（5-6）（）（5-8）的方法检验）的方法检验。若满足精度，则以。若满足精度，则以 为计算结果；若不满足精度，则分成两个小区为计算结果；若不满足精度，则分成两个小区间间各自重复逐步上述过程，每个小区间精度为各自重复逐步上述过程，每个小区间精度为352021-11-154 newton-cotes4 newton-cotes求积公式的求积公式的matlabmatlab实现实现n常用的三种常用的三种newton-cotes系列数值积分法的相应系列数值积分法的相应matlab函函数如下数如下:n1) 复合梯形法数值积分：复合梯

37、形法数值积分：trapz( )n调用形式：调用形式：z = trapz(x,y)n其中其中: x,y分别代表数目相同的向量或数组，而分别代表数目相同的向量或数组，而y 与与x 的的 关系可以是一关系可以是一 函数型态（如函数型态（如y=sin(x)）或是不以函数描）或是不以函数描述的离散型态；述的离散型态； z代表返回的积分值；代表返回的积分值；注注:离散型态数据用离散型态数据用trapz函数时，还需设定函数时，还需设定x在区间在区间 a,b 之间离之间离散点的间隔；缺省参数散点的间隔；缺省参数x时，表示时，表示x被等分，每份宽为被等分，每份宽为1。362021-11-152） 自适应自适应s

38、impson法数值积分：法数值积分：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表示用图形表示用图形显示

39、积分过程显示积分过程 p1,p2,直接传递给函数直接传递给函数fun的参数。的参数。3）自适应）自适应lobatto法数值积分：法数值积分：quadl（）（） quadl是高阶的自适应是高阶的自适应newton-cotes数值积分法函数，它数值积分法函数，它比比quad函数更有效，精度更高。其使用方法和函数更有效，精度更高。其使用方法和quad（）完（）完全相同。全相同。需要了解更多的内容可查阅需要了解更多的内容可查阅matlab功能函数库（功能函数库（funfun）。）。372021-11-15flglg2.303afprt例例5.4：真实气体的逸度真实气体的逸度可用下式计算：可用下式计算：

40、 现测得现测得00c下氢气的有关数值如表下氢气的有关数值如表5-5所示，试求所示，试求1000 atm下的逸度。下的逸度。6310vmrtvpp6310vm表表5-5 00c下氢气的相关数值下氢气的相关数值 （atm）（atm）015.4660053.4316.09100239.5115.4670048.1416.13200127.4915.4680044.1716.1630090.2915.6190041.0616.1640071.8615.85100038.5516.1450060.7615.93-真实气体的实测体积和按理想气体真实气体的实测体积和按理想气体定律计算得到的体积之间的差值。定

41、律计算得到的体积之间的差值。r-气体常数气体常数 0,prtadpvpfp6382.06 10/()matm mol k其中：其中： -逸度；逸度； -压力，压力，atm；t-绝对温度，绝对温度，k。 r tvpfrtvpp382021-11-15解：解：本题是离散型数据，可用本题是离散型数据，可用trapz函数求解数值积分。函数求解数值积分。（1）计算数值积分，程序如下：）计算数值积分，程序如下：%计算数值积分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);% 梯形法计算数值积分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.00 atmf392021-11-150.4,0.9fdxx5r () /dfxxddxnyxxyrxy例例5.5：氯仿：氯仿-苯双组分精馏系统的气液平衡数据如表苯双组分精馏系统的气液平衡数据如表5-6所示。规定所示。规定进料和塔顶的组成分别是进料和塔顶的组成分别是

43、 ，精馏段的回流比为，精馏段的回流比为精馏段理论板数的模型为精馏段理论板数的模型为表表5-6 精馏段气液平衡数据精馏段气液平衡数据0.1780.2750.3720.4560.6500.8440.2430.3820.5180.6160.7950.931，试用试用matlab计算所需的精馏段理论板数。计算所需的精馏段理论板数。402021-11-15n解：解：因模型中的的函数关系是以表格形式给出，我们若要因模型中的的函数关系是以表格形式给出，我们若要用辛普森法等计算较精确的精馏段理论板数的时候，需先用辛普森法等计算较精确的精馏段理论板数的时候，需先采用拟合法将离散数据（采用拟合法将离散数据（ ，

44、）拟合成多项式，再将多）拟合成多项式，再将多项式代入被积函数求积。项式代入被积函数求积。计算程序如下：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（）拟合多项式% 画出拟合曲线，直观比较拟合效果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); % 计算精馏段理论板数n=round(n+0.5) % 数据取整function f=func1(x,sp)y=fnval(sp,x);f=1./(y-x-(0.9-y)./5);执行结果：执行结果：n = 5所以，所以，精馏段理论板数为精馏段理论板数为5块。块。ixiy412021-11-15 5.2.2 高斯勒让德公式高斯勒让德公式高斯法是一种精度较高的求积分法，它的高斯法是一种精度较高的求积分法，它的一般公式是一般公式是式中式中(x)是一个权重函数，是一个权重函数，aj为系数，为系数，xj为横坐标上的节点为横坐标上的节点 高斯勒让德多项式的权重函数高斯勒让德多项式的权重函数(x)=1，因而，一个，因而，一个n点高斯勒让德求积公点高斯勒让德求积公式具有如下形式：式具有如下形式：banjjjxfadxxfx1)()()(111( )()njjjf x dxa f x1 方法概述方法概述 422021-11-15 右边的右边的f(xj)是函数是函数f(x)在节点在节点xj处的值，处的值， 节点节点xj是