1、2022-12-111第第5 5章章 数值微分与数值积分数值微分与数值积分22022-12-11n微积分是高等数学中的重要内容，在化学微积分是高等数学中的重要内容，在化学工程上有许多非常重要的应用工程上有许多非常重要的应用n微积分的数值方法，不同于高等数学中的微积分的数值方法，不同于高等数学中的解析方法，尤其适合求解没有或很难求出解析方法，尤其适合求解没有或很难求出微分或积分表达式的实际化工问题的计算，微分或积分表达式的实际化工问题的计算，例如：列表函数求微分或积分例如：列表函数求微分或积分 引言引言1-1-数值微积分方法不同于解析方数值微积分方法不同于解析方法法32022-12-11n数值微
2、分和数值积分与插值和拟合往往是密数值微分和数值积分与插值和拟合往往是密不可分的不可分的n如在进行数值微分时，常针对离散的数据点，如在进行数值微分时，常针对离散的数据点，利用插值和拟合可以减少误差；而数值积分利用插值和拟合可以减少误差；而数值积分的基本思路也来自于插值法。例如如果所积的基本思路也来自于插值法。例如如果所积函数的形式比较复杂或以表格形式给出，则函数的形式比较复杂或以表格形式给出，则可通过构造一个插值多项式来代替原函数，可通过构造一个插值多项式来代替原函数，从而使问题简化从而使问题简化引言引言2-2-数值微分和数值积分与插值和拟数值微分和数值积分与插值和拟合关系密切合关系密切4202
3、2-12-115.1 5.1 数值微分数值微分n化工领域的实际问题中时常需要求列表函数在节化工领域的实际问题中时常需要求列表函数在节点和非节点处的导数值，这是数值微分所要解决点和非节点处的导数值，这是数值微分所要解决的问题。数值微分方法可近似求出某点的导数值的问题。数值微分方法可近似求出某点的导数值n例如在反应动力学的研究中，根据实验数据确定反应的动例如在反应动力学的研究中，根据实验数据确定反应的动力学方程力学方程 ：n这里实验测得一批离散点，要计算这里实验测得一批离散点，要计算 只能借助数值只能借助数值微分求导解决微分求导解决 表表5-1 5-1 反应动力学实验数据反应动力学实验数据 mAp
4、Adpk pdtAdpdtt1t2tnpA1pA2PAn0 0 导入导入52022-12-110.1 0.1 建立数值微分公式的三种思路建立数值微分公式的三种思路 常用三种思路建立数值微分公式：常用三种思路建立数值微分公式：1.1.从微分定义出发，通过近似处理，得到数值微从微分定义出发，通过近似处理，得到数值微分的近似公式分的近似公式2.2.从插值近似公式出发，对插值公式的近似求导从插值近似公式出发，对插值公式的近似求导可得到数值微分的近似公式可得到数值微分的近似公式3.3.先用最小二乘拟合方法根据已知数据得近似函先用最小二乘拟合方法根据已知数据得近似函数，再对此近似函数求微分可得到数值微分的
5、数，再对此近似函数求微分可得到数值微分的近似公式。然后对各方法数值微分后得到的多近似公式。然后对各方法数值微分后得到的多项式求值，即可求出任意点处的任意阶微分项式求值，即可求出任意点处的任意阶微分62022-12-11000()()()()()()()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 5.1.1 差分近似微分差分近似微分上式中三种不同表示形式依次是一阶前向差分、一阶后向差分和一阶中上式中三种不同表示形式依次是一阶前向差分、一阶后向差分和一阶中心差分来近似表示微分。其中一阶中心差分的精度较高。心差分来近似表示微分。其中一阶中心差分的精度较高。72022-12-
7、112 2 差分的差分的MatlabMatlab实现实现n在在MatlabMatlab中，可用中，可用diffdiff函数进行离散数据的近似函数进行离散数据的近似求导求导n调用形式调用形式：Y=diff(X,n)Y=diff(X,n)n其中：其中：X X表示求导变量，可以是向量或矩阵。如表示求导变量，可以是向量或矩阵。如是矩阵形式则按各列作差分；是矩阵形式则按各列作差分；n n表示表示n n阶差分，即差分阶差分，即差分n n次次 Y Y是是X X的差分结果的差分结果n注：注：用用diffdiff函数进行离散数据的近似求导与前向函数进行离散数据的近似求导与前向差分近似，差分近似，但误差较大。最好
8、将数据利用插值或但误差较大。最好将数据利用插值或拟合得到多项式，然后对近似多项式进行微分拟合得到多项式，然后对近似多项式进行微分82022-12-11464622()C HC H ApAdpdt例例5.5.1 1：丁二烯的气相二聚反应如下：丁二烯的气相二聚反应如下：实验在一定容器的反应器中进行，实验在一定容器的反应器中进行，3263260 0C C时，测得物系中丁二烯的时，测得物系中丁二烯的分压分压(mmHg)(mmHg)与时间的关系如表与时间的关系如表5-25-2所示。所示。ApAp表表 5-2 5-2 丁二烯二聚反应实验数据丁二烯二聚反应实验数据t(min)(mmHg)t(min)(mmH
9、g)0632.050362.05590.055348.010552.060336.015515.065325.020485.070314.025458.075304.030435.080294.035414.085284.040396.090274.045378.0用数值微分法计算所列时刻每一瞬间的反应速率用数值微分法计算所列时刻每一瞬间的反应速率 92022-12-11解：解：程序如下：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 30
10、4.0 294.0 284.0 274.0;dt=diff(t);求时间t的差分dpA=diff(pA);求压力的差分q=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.000010202
11、2-12-115.1.2 5.1.2 三次样条插值函数求微分三次样条插值函数求微分()S x()f x()S x()fx()S x()f x若三次样条插值函数若三次样条插值函数收敛于收敛于，那么导数，那么导数收敛于收敛于，因此用样条插值函数，因此用样条插值函数作为作为函数，不但彼此的函数值非常接近，而且导数值也很接近。函数，不但彼此的函数值非常接近，而且导数值也很接近。用三次样条插值函数求数值导数是可靠的，这是化工计用三次样条插值函数求数值导数是可靠的，这是化工计算中求数值微分的有效方法。算中求数值微分的有效方法。的近似的近似112022-12-11()()fxS x1,2,in1,iixxx
12、11()()()()iiiiiixxxxfxS xMMhh221111()()()()()226iiiiiiiiiiiiyyxxxxhfxS xMMMMhhh 1 1 方法概述方法概述用三次样条插值函数建立的数值微分公式为：用三次样条插值函数建立的数值微分公式为：求导得：求导得：（5-25-2）;式（式（5-15-1）和（）和（5-25-2）不但适用于求节点处的导数，而且可求非节点处的导数。）不但适用于求节点处的导数，而且可求非节点处的导数。（5-15-1）其中，其中，122022-12-112 2 三次样条插值函数求微分的三次样条插值函数求微分的MatlabMatlab函数函数MatlabM
13、atlab求离散数据的三次样条插值函数微分方法分三个步骤：求离散数据的三次样条插值函数微分方法分三个步骤：Step 1:Step 1:对离散数据用对离散数据用csapicsapi函数（或函数（或splinespline函数）得到其三次样条插值函数函数）得到其三次样条插值函数调用形式调用形式 pp=csapi(x,y)pp=csapi(x,y)其中：其中：x,yx,y分别为离散数据对的自变量和因变量；分别为离散数据对的自变量和因变量；pppp为得到的三次样条插值函数。为得到的三次样条插值函数。Step 2:Step 2:用用fnderfnder函数求三次样条插值函数的导数函数求三次样条插值函数的
14、导数调用形式调用形式 fprime=fnder(f,dorder)fprime=fnder(f,dorder)其中：其中：f f为三次样条插值函数；为三次样条插值函数；dorderdorder为三次样条插值函数的求导阶数；为三次样条插值函数的求导阶数；fprimefprime为得到的三次样条插值函数的导函数。为得到的三次样条插值函数的导函数。Step3:Step3:可用可用fnvalfnval函数求导函数在未知点处的导数值函数求导函数在未知点处的导数值调用形式调用形式 v=fnval(fprime,x)v=fnval(fprime,x)其中：其中：fprimefprime为三次样条插值函数导函
15、数；为三次样条插值函数导函数；x x为未知点处自变量值；为未知点处自变量值；v v为未知点处的导数值。为未知点处的导数值。132022-12-11例例5.2 5.2 ：某液体冷却时，温度随时间的变化数据如表某液体冷却时，温度随时间的变化数据如表5-35-3所示：所示：表表5-3 5-3 冷却温度随时间的变化数据冷却温度随时间的变化数据试分别计算试分别计算t t2 2，3 3，4min4min及及t t1.51.5，2.52.5，4.5min4.5min时的降温速率。时的降温速率。解：解：三次样条插值函数求数值微分的程序如下：三次样条插值函数求数值微分的程序如下：t=0:5;T=92,85.3,
16、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.5000 2.5000 4.5000 2.0000 3.0000 4.0000 1.5000 2.5000 4.5000 -5.3722 -4.6722 -3.8389 -5.7972 -4.
17、9889 -3.2222 -5.3722 -4.6722 -3.8389 -5.7972 -4.9889 -3.2222t(min)012345T(0C)92.085.379.574.570.267.0注：注：前者是计算节点处的一阶导数，后者是计算非节点处的一阶导数前者是计算节点处的一阶导数，后者是计算非节点处的一阶导数142022-12-115.1.3 5.1.3 最小二乘法样条拟合函数求最小二乘法样条拟合函数求微分微分n在实际化工应用中，当来自实验观测的离在实际化工应用中，当来自实验观测的离散数据不可避免地含有较大随机误差时，散数据不可避免地含有较大随机误差时，此时用插值公式求数值微分虽然
18、样本点处此时用插值公式求数值微分虽然样本点处误差较小，但可能会使非样本点处产生较误差较小，但可能会使非样本点处产生较大误差大误差n为此，可采用最小二乘法样条拟合实验数为此，可采用最小二乘法样条拟合实验数据，获得一个函数模型，然后再对其求导据，获得一个函数模型，然后再对其求导数。与样条插值不同，样条拟合不要求曲数。与样条插值不同，样条拟合不要求曲线经过全部的数据点，这样处理求导结果线经过全部的数据点，这样处理求导结果会有很大改善会有很大改善152022-12-11,iix y()yf x2min()iiiiEw yf xiw1 1 方法概述方法概述 可用最小二乘法拟合成平滑可用最小二乘法拟合成平
19、滑B B样条曲线，即对于离散数据（样条曲线，即对于离散数据（所求的所求的K K次样条拟合函数次样条拟合函数满足：满足：其中：其中：再对拟合函数作平滑处理后求导，即可求出任意点处微分。再对拟合函数作平滑处理后求导，即可求出任意点处微分。）为权重系数，默认为为权重系数，默认为1 1162022-12-11nMatlabMatlab求离散数据的最小二乘法平滑求离散数据的最小二乘法平滑B B样条拟合函数求微分样条拟合函数求微分共三个步骤：共三个步骤：Step 1:Step 1:对离散数据用对离散数据用spapsspaps函数得到最小二乘平滑函数得到最小二乘平滑B B样条拟合样条拟合函数。函数。调用格式
20、：调用格式：sp=spaps(x,y,tol)sp=spaps(x,y,tol)其中：其中：x,yx,y要处理的离散数据（）要处理的离散数据（）toltol光滑时的允许精度，通常取（光滑时的允许精度，通常取（10-210-210-410-4）Step 2:Step 2:可用可用fnderfnder函数求最小二乘平滑函数求最小二乘平滑B B样条拟合函数的导数样条拟合函数的导数;Step3:Step3:可用可用fnvalfnval函数求导函数在未知点处的导数值。函数求导函数在未知点处的导数值。fnder()fnder()和和fnvalfnval（）调用形式前文中已经介绍过。（）调用形式前文中已经介
21、绍过。2 2 最小二乘法平滑最小二乘法平滑B B样条拟合函数求微分的样条拟合函数求微分的MatlabMatlab实现实现172022-12-11n由离散数据求数值微分的四种方法及有关由离散数据求数值微分的四种方法及有关MATLABMATLAB函数：函数：（1 1）差分法）差分法 用差分函数用差分函数diff()diff()近似计算导数，即近似计算导数，即dy=diff(y)./diff(x)dy=diff(y)./diff(x)。对于向量对于向量X X，diff(X)diff(X)表示了表示了X(2)-X(1)X(3)-X(2).X(n)-X(n-1).X(2)-X(1)X(3)-X(2).X
22、(n)-X(n-1).对于矩阵对于矩阵X X，diff(X)diff(X)表示了表示了X(2:n,:)-X(1:n-1,:)X(2:n,:)-X(1:n-1,:)小结小结注：注：此法用一阶差分，精度较差，若改用二阶差分，可大大提高精度，但编程麻烦些。此法用一阶差分，精度较差，若改用二阶差分，可大大提高精度，但编程麻烦些。()polyfit0(polyder()polyval（2 2）多项式拟合方法）多项式拟合方法向量向量p p表示的多项式拟合函数表示的多项式拟合函数导函数导函数pppppppp在在xi xi的导数值。的导数值。其中：其中：函数函数polyfit()polyfit()和和poly
23、val()polyval()在前文中已介绍在前文中已介绍导函数导函数polyder()polyder()的的调用格式调用格式为：为：pp=polyder(p)pp=polyder(p)离散数据离散数据 该函数对向量该函数对向量p p表示的多项式函数进行求导，返回导函数为表示的多项式函数进行求导，返回导函数为pppp。182022-12-11（3 3）三次样条插值方法）三次样条插值方法 （4 4）样条拟合方法（最小二乘法）样条拟合方法（最小二乘法）其中：其中：函数函数csaps()csaps()、spap2()spap2()、spaps()spaps()和和fnval()fnval()已在前已在
24、前文中介绍，文中介绍，求导函数求导函数fnder()fnder()的的调用格式调用格式为：为：pp=fnder(f,dorder)pp=fnder(f,dorder)该函数计算函数该函数计算函数f f的的dorderdorder阶导函数，默认阶数阶导函数，默认阶数dorder=1dorder=1。()csapi()fnder()polyval离散数据离散数据三次样条插值函数三次样条插值函数cscscscs的导数的导数pppppppp在在xi xi的导数值的导数值()()2()spapsaspapcsaps或或()fnder()fnval 离散数据离散数据样条拟合函数样条拟合函数spspspsp
25、的导数的导数pppppppp在在xi xi的导数值。的导数值。192022-12-11 ACACmAAdCk Cdt例例5.5.3 3：反应物反应物A A在一等温间歇反应器中发生的反应为：在一等温间歇反应器中发生的反应为：A A测量得到的反应器中不同时间下反应物测量得到的反应器中不同时间下反应物A A的浓度的浓度表表5-4 5-4 间歇反应器动力学数据间歇反应器动力学数据t(s)0204060120180300 (mol/L)10865321系统的动力学模型为：系统的动力学模型为：，试根据表中数据确定其反应速率方程。，试根据表中数据确定其反应速率方程。产物产物如表如表5-45-4所示。所示。2
26、02022-12-11ln()lnlnAAdCmCkdtln(),lnAAdCyxCdtlnykmxAdCdt解：解：（1 1）系统的动力学模型为非线性形式，可将其线性化。）系统的动力学模型为非线性形式，可将其线性化。对方程两边取对数：对方程两边取对数：令令 则原模型变为：则原模型变为：（2 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=l
27、inspace(t(1),t(end),200);CAi=fnval(sp,ti);plot(t,CA,bo,ti,CAi,r-),xlabel(t),ylabel(CA)及得到拟合曲线的图形程序如下：及得到拟合曲线的图形程序如下：212022-12-11（3 3）线性拟合程序如下：）线性拟合程序如下：y=log(-dCAdt);x=log(CA);p=polyfit(x,y,1);k=exp(p(2),m=p(1)执行结果：执行结果：k=k=0.0059 0.0059 m=m=1.2904 1.2904所以本例的反应速率方程为：所以本例的反应速率方程为：1.29040.0059AAdCCdt
28、222022-12-111)1)被积函数以一组数据形式表示；被积函数以一组数据形式表示；2)2)被积函数过于特殊或原函数不能用初等函数表示，积分表中无被积函数过于特殊或原函数不能用初等函数表示，积分表中无 法找到可沿用的现成公式；法找到可沿用的现成公式；3)3)有的原函数十分复杂难以计算。有的原函数十分复杂难以计算。badxxffI)()(对于积分：对于积分：公式有则由的原函数如果知道LeibnizNewtonxFxf),()(badxxf)()()()(aFbFxFba但是在工程技术和科学研究中但是在工程技术和科学研究中,常会见到以下现象常会见到以下现象:0 0 导入导入5.2 5.2 数值
29、积分数值积分232022-12-11 数值积分数值积分就是一种常用的近似计算方法。数就是一种常用的近似计算方法。数值积分方法不受以上几个问题的限制，在化学化值积分方法不受以上几个问题的限制，在化学化工领域应用甚广，如反应热效应计算、热容计算、工领域应用甚广，如反应热效应计算、热容计算、熵的计算、反应活化能的计算等。熵的计算、反应活化能的计算等。以上这些现象以上这些现象,Newton-Leibniz,Newton-Leibniz很难发挥作用，很难发挥作用，只能建立积分的近似计算方法只能建立积分的近似计算方法242022-12-11如：如：)(,H21TfCpdTCpntt么么么么方面nSds绝对
30、是假的262022-12-110.1 0.1 数值积分的基本思路和方法数值积分的基本思路和方法n 常用的数值积分的常用的数值积分的基本思路基本思路来自于插值法，它来自于插值法，它通过构造一个插值多项式通过构造一个插值多项式 作为作为 的近似表的近似表达式，用达式，用 的积分值作为的积分值作为 的近似积分的近似积分值。值。n数值积分的数值积分的方法方法很丰富，常用的插值型求积公很丰富，常用的插值型求积公式有两类：一类是等距节点的牛顿柯特斯求式有两类：一类是等距节点的牛顿柯特斯求积公式；积公式；另一类是不等距节点的高斯型求积另一类是不等距节点的高斯型求积公式。公式。()nP x()f x()nP
31、x()f x272022-12-11 Newton-CotesNewton-Cotes公式是指公式是指等距节点下等距节点下使用使用LagrangeLagrange插值插值 多项式多项式建立的数值求积公式建立的数值求积公式,)(baCxf设函数等份分割为将积分区间nba,nkkhaxk,1,0,为步长其中nabh各节点为各节点为插值多项式为的Lagrangexf)(则：则：5.2.15.2.1牛顿柯特斯（牛顿柯特斯（Newton-CotesNewton-Cotes）求积公式）求积公式1 1 方法概述方法概述282022-12-11这里这里 是既不依赖于被积函数，也不依赖于积分区间的常数，是既不依
32、赖于被积函数，也不依赖于积分区间的常数，称为柯特斯系数。式（称为柯特斯系数。式（5-35-3）称为牛顿柯特斯求积公式。）称为牛顿柯特斯求积公式。00()()()()nnbbkkkkaakkIf x dxf x lx dxA f x0()bbjkkaajnkjjkxxAlx dxdxxxxath/()iiCAba00(1)()!()!nknijnjkCtj dtn knkiC其中：其中：一般地：令一般地：令 ，得：得：（5-35-3）292022-12-11在在Newton-CotesNewton-Cotes公式中公式中,n=1,2,4,n=1,2,4时的公式是最常用也是最时的公式是最常用也是最
33、重要三个公式重要三个公式,称为低阶公式称为低阶公式（当（当n=0n=0时的公式为矩形公式）时的公式为矩形公式）1 1）梯形）梯形(trapezia)(trapezia)公式公式abhbxaxn,110则取dtt10)1()1(0CCotesCotes系数为系数为21dtt10)1(1C21求积公式为求积公式为302022-12-111(1)1010()()()()()2kkkbaIfbaCf xf xf x)()(2bfafab)(1fI即上式称为上式称为梯形求积公式梯形求积公式,（5-45-4）312022-12-112 2）SimpsonSimpson公式公式2,2,2210abhbxab
34、xaxn则取CotesCotes系数为系数为dtttC20)2(0)2)(1(4161dtttC20)2(1)2(2164dtttC20)2(2)1(4161322022-12-11)(61)(64)(61)(210 xfxfxfab)()2(4)(6bfbafafab)(2fI上式称为上式称为SimpsonSimpson求积公式求积公式，也称，也称三点公式或抛物线公式三点公式或抛物线公式求积公式为求积公式为2I20)2()()(kkkxfCab式（式（5-45-4）和式（）和式（5-55-5）是化工领域常用的两个求积公式。与梯形法）是化工领域常用的两个求积公式。与梯形法求积公式相比，求积公式
35、相比，SimpsonSimpson法求积公式是一个较高精度的求积公式。法求积公式是一个较高精度的求积公式。用式（用式（5-35-3）还可得到更高阶数）还可得到更高阶数（5-55-5）332022-12-11n Newton-CotesNewton-Cotes公式当公式当n n大于大于7 7时，公式的稳定性时，公式的稳定性将无法保证，因此将无法保证，因此,在实际应用中一般不使用高阶在实际应用中一般不使用高阶Newton-CotesNewton-Cotes公式而是采用公式而是采用低阶复合求积法低阶复合求积法n在实际计算中为了保证计算的精度，往往首先用在实际计算中为了保证计算的精度，往往首先用分点分
36、点x xk k=a a+khkh,(k k=0,1,=0,1,n n)将区间将区间a a,b b分成分成n n个个相等的子区间，而后对每个子区相等的子区间，而后对每个子区 间再应用梯形公间再应用梯形公式或式或SimpsonSimpson公式，分别得到：公式，分别得到：nabh11)(2)()(2 复化梯形公式nknkhafbfafhT2 2 复化法求积公式复化法求积公式11102 ()4()()6nnkkkkhSf xf xf x复化复化SimpsonSimpson公式：公式：342022-12-11 尽管复化法求积公式具有很高的精度，但是它必须采用等步长方法，尽管复化法求积公式具有很高的精度
37、，但是它必须采用等步长方法，从而限制了它的效率。这里我们介绍一种更加灵活选取步长的方法，即从而限制了它的效率。这里我们介绍一种更加灵活选取步长的方法，即自适应步长法。自适应步长法。（5-85-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积分积分 达到精度达到精度时，可认为区间时，可认为区间取计算需满足的精度为取计算需满足的精度为 以以Simpso
38、n积分法为例，某区间积分法为例，某区间 ，记，记 ，考虑该区间，考虑该区间上的上的Simpson积分和二等分以后的两个积分和二等分以后的两个Simpson积分和：积分和：3 3 自适应求积公式自适应求积公式352022-12-11()/2kkhba/2k 这样重复下去，直至每个分段部分达到相应精度（步长为这样重复下去，直至每个分段部分达到相应精度（步长为时精度为时精度为 ）。这样，不同段的步长可能是不一样的，积分结果为）。这样，不同段的步长可能是不一样的，积分结果为每一小段的积分总和。每一小段的积分总和。,a b2S/2 自适应步长自适应步长Simpson法从法从 开始按式（开始按式（5-6）
39、（）（5-8）的方法检验）的方法检验。若满足精度，则以。若满足精度，则以 为计算结果；若不满足精度，则分成两个小区为计算结果；若不满足精度，则分成两个小区间间各自重复逐步上述过程，每个小区间精度为各自重复逐步上述过程，每个小区间精度为362022-12-114 Newton-Cotes4 Newton-Cotes求积公式的求积公式的MatlabMatlab实现实现n常用的三种常用的三种Newton-CotesNewton-Cotes系列数值积分法的相应系列数值积分法的相应MatlabMatlab函数如下函数如下:n1)1)复合梯形法数值积分：复合梯形法数值积分：trapz()trapz()n调
40、用形式：调用形式：Z=trapz(X,Y)Z=trapz(X,Y)n其中其中:X,YX,Y分别代表数目相同的向量或数组，而分别代表数目相同的向量或数组，而Y Y 与与X X 的的 关系可以是一关系可以是一 函数型态（如函数型态（如y y=sin(=sin(x x)）或是不以函数描述）或是不以函数描述的离散型态；的离散型态；Z Z代表返回的积分值；代表返回的积分值；注注:离散型态数据用离散型态数据用trapztrapz函数时，还需设定函数时，还需设定X X在区间在区间 a,ba,b 之间离之间离散点的间隔；缺省参数散点的间隔；缺省参数X X时，表示时，表示X X被等分，每份宽为被等分，每份宽为1
41、 1。372022-12-112 2）自适应自适应SimpsonSimpson法数值积分：法数值积分：quadquad（）（）基本调用格式基本调用格式：q=quad(fun,a,b)q=quad(fun,a,b)或或 q=quad(fun,a,b,tol,trace,p1,p2,)q=quad(fun,a,b,tol,trace,p1,p2,)其中其中:fun fun被积函数。可以是内置函数、被积函数。可以是内置函数、mm文件或函数句柄文件或函数句柄,函数表达式函数表达式 中的必须使用点运算符号。中的必须使用点运算符号。a,ba,b分别是积分的下限和上限；分别是积分的下限和上限；q q积分结果
42、。积分结果。toltol默认误差限，默认值为默认误差限，默认值为1.e-6.1.e-6.trace-trace-取取0 0表示不用图形显示积分过程，非表示不用图形显示积分过程，非0 0表示用图形表示用图形显示积分过程显示积分过程 p1,p2,p1,p2,直接传递给函数直接传递给函数funfun的参数。的参数。3 3）自适应）自适应LobattoLobatto法数值积分：法数值积分：quadlquadl（）（）quadlquadl是高阶的自适应是高阶的自适应Newton-CotesNewton-Cotes数值积分法函数，它数值积分法函数，它比比quadquad函数更有效，精度更高。其使用方法和函
43、数更有效，精度更高。其使用方法和quadquad（）完（）完全相同。全相同。需要了解更多的内容可查阅需要了解更多的内容可查阅MatlabMatlab功能函数库（功能函数库（funfunfunfun）。）。382022-12-11flglg2.303AfPRT例例5.5.4 4：真实气体的逸度真实气体的逸度可用下式计算：可用下式计算：现测得现测得0 00 0C C下氢气的有关数值如表下氢气的有关数值如表5-55-5所示，试求所示，试求1000 atm1000 atm下的逸度。下的逸度。6310VmRTVPP6310Vm表表5-5 05-5 00 0C C下氢气的相关数值下氢气的相关数值 （atm
44、）（atm）015.4660053.4316.09100239.5115.4670048.1416.13200127.4915.4680044.1716.1630090.2915.6190041.0616.1640071.8615.85100038.5516.1450060.7615.93-真实气体的实测体积和按理想气体真实气体的实测体积和按理想气体定律计算得到的体积之间的差值。定律计算得到的体积之间的差值。R-气体常数气体常数 0,pRTAdPVPfP6382.06 10/()matm mol K其中：其中：-逸度；逸度；-压力，压力，atm；T-绝对温度，绝对温度，K。R TVPfRTVP
45、P392022-12-11解：解：本题是离散型数据，可用本题是离散型数据，可用trapztrapz函数求解数值积分。函数求解数值积分。（1 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=trapz(P,a1);%梯形法计算数值积分A=-A执行结果：执行结果：A=A=15865 15865（2 2）计算逸度值，程序如下：）计算逸度值，程序如下：%计算逸度lf=log10(1000)+A./(2.
46、303.*82.06.*273.2);%lf代表f=10.lf执行结果：执行结果：f=f=2.0290e+003 2.0290e+003所以所以 f =2029.00 atm f =2029.00 atmf402022-12-110.4,0.9fdxx5R()/dfxxddxNyxxyRxy例例5.5.5 5：氯仿：氯仿-苯双组分精馏系统的气液平衡数据如表苯双组分精馏系统的气液平衡数据如表5-65-6所示。规所示。规定进料和塔顶的组成分别是定进料和塔顶的组成分别是 ，精馏段的回流比为，精馏段的回流比为精馏段理论板数的模型为精馏段理论板数的模型为表表5-6 5-6 精馏段气液平衡数据精馏段气液平
47、衡数据0.1780.2750.3720.4560.6500.8440.2430.3820.5180.6160.7950.931，试用试用MatlabMatlab计算所需的精馏段理论板数。计算所需的精馏段理论板数。412022-12-11n解：解：因模型中的的函数关系是以表格形式给出，我们若要因模型中的的函数关系是以表格形式给出，我们若要用辛普森法等计算较精确的精馏段理论板数的时候，需先用辛普森法等计算较精确的精馏段理论板数的时候，需先采用拟合法将离散数据（采用拟合法将离散数据（，）拟合成多项式，再将多）拟合成多项式，再将多项式代入被积函数求积。项式代入被积函数求积。计算程序如下：functio
48、n 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.4,0.9,sp);%计算精馏段理论板数N=round(N+0.5)%数据取整function f=func1(x,s
49、p)y=fnval(sp,x);f=1./(y-x-(0.9-y)./5);执行结果：执行结果：N=N=5 5所以，所以，精馏段理论板数为精馏段理论板数为5 5块。块。ixiy422022-12-11 5.2.2 5.2.2 高斯勒让德公式高斯勒让德公式高斯法是一种精度较高的求积分法，它的高斯法是一种精度较高的求积分法，它的一般公式是一般公式是式中式中(x)(x)是一个权重函数，是一个权重函数，A Aj j为系数，为系数，x xj j为为横坐标上的节点横坐标上的节点 高斯勒让德多项式的权重函数高斯勒让德多项式的权重函数(x)=1(x)=1，因而，一个，因而，一个n n点高斯勒让德求积公点高斯勒
50、让德求积公式具有如下形式：式具有如下形式：banjjjxfAdxxfx1)()()(111()()njjjf x dxA f x1 1 方法概述方法概述 432022-12-11 右边的右边的f(xf(xj j)是函数是函数f(x)f(x)在节点在节点x xj j处的值，处的值，节点节点x xj j是勒让德多项式是勒让德多项式P Pn n(x)(x)的根。的根。A Aj j 为系数，为系数，其值为其值为)1()(222jjnjxxpA111)()(njjjxfAdxxf式中式中P Pn n(x x)是勒让德多项式是勒让德多项式P Pn n(x x)的一阶导数的一阶导数442022-12-11p
