1、5 数值微分与积分数值微分与积分5.1 数值微分数值微分 根据离散点上的函数值求取某点导数近似根据离散点上的函数值求取某点导数近似值的方法。值的方法。1）差商代替导数）差商代替导数 2）插值型数值求导）插值型数值求导 例例:已知自变量和函数值如下，已知自变量和函数值如下，求各点的导数值。求各点的导数值。x00.10.20.30.40.50.60.70.80.91.0f(x)0.48 0.38 0.31 0.33 0.36 0.41 0.51 0.43 0.35 0.29 0.28xxfxfdxdfiixxi )()(1向前差商向前差商 用差商近似代替导数，用差商近似代替导数，f（x）在）在xi
2、处的差商形式为处的差商形式为xxfxfdxdfiixxi )()(1向后差商向后差商xxfxfdxdfiixxi211)()(中心差商中心差商yiyi-1yi+1yxxi-1xixi+1向前向前向后向后5 数值微分与积分数值微分与积分 21211)(21)()()(2iixxiixxiixxdxfdxxdxdfxfxfii x=h(x)2=h2f(xi)在其相邻点在其相邻点 xi+1处可展开成泰勒级数处可展开成泰勒级数 考查差商代替导数的误差考查差商代替导数的误差 xdxfdxxfxfdxdfiixxiixx22121)()(整理后得向前差商整理后得向前差商右端除第一项外其后各项省略右端除第一
3、项外其后各项省略(O)x#)(O)()(1xxxfxfdxdfiixxi O(h)截断误差：步截断误差：步长长 h 的一次方量级（一的一次方量级（一阶精度）阶精度）数值微分数值微分2整理可得向后差商整理可得向后差商 21211)(21)()()(2iixxiixxiixxdxfdxxdxdfxfxfii在相邻点在相邻点 xi-1处可展开成泰勒级数处可展开成泰勒级数)(O)()(1xxxfxfdxdfiixxi O(h)一阶精度一阶精度)(O2)()(211xxxfxfdxdfiixxi xxfxfxfdxdfxx )(5.0)(5.1)(23121将将 xi+1和和xi-1处式处式泰勒级数展开
4、式泰勒级数展开式相减得中心差商相减得中心差商O(h2)二阶精度二阶精度 i=1和和 i=n（端点），二阶精度的数值微分式分别为（端点），二阶精度的数值微分式分别为xxfxfxfdxdfnnnxxn )(5.0)(2)(5.121数值微分数值微分2 2)插值型数值求导插值型数值求导 用插值多项式用插值多项式Ln(x)近似近似f(x)，对，对Ln(x)求导作为求导作为f(x)的导数的近似值，称为插值型求导公式。的导数的近似值，称为插值型求导公式。)()(xLxf )()()(11111 kkkkkkkkxfxxxxxfxxxxxL线性插值型公式对线性插值型公式对x求导求导hxfxfxLkk)()(
5、)(11 hxfxfxfkkk)()()(1 hxfxfxfkkk)()()(11 两点公式两点公式 两点的求导公式形式相同，两点公式的截断误差两点的求导公式形式相同，两点公式的截断误差可用插值余项公式可用插值余项公式(3-18)计算，为计算，为O(h)（一阶精度）（一阶精度）同理可由三点插值公式求导得三点公式同理可由三点插值公式求导得三点公式 数值微分数值微分2 上式求导可得三点求导公式上式求导可得三点求导公式)()()()()()()()()()(111111111111112 kkkkkkkkkkkkkkkkkkkkkxfxxxxxxxxxfxxxxxxxxxfxxxxxxxxxLhxf
6、xfxfxfkkkk2)()(4)(3)(111 hxfxfxfkkk2)()()(11 hxfxfxfxfkkkk2)(3)(4)()(111 三点求导公式的截断三点求导公式的截断误差为误差为O(h2)（二阶精（二阶精度）度）二阶导数的三点公式二阶导数的三点公式2111)()(2)()(hxfxfxfxfkkkk 211)()(2)()(hxfxfxfxfkkkk 2111)()(2)()(hxfxfxfxfkkkk 从截断误差的角度看，减小步长从截断误差的角度看，减小步长h会提高计算的会提高计算的精度，但另一方面，精度，但另一方面，h的减小会使舍入误差增大。的减小会使舍入误差增大。由于由于
7、数值微分计算对舍入误差非常敏感，因此步数值微分计算对舍入误差非常敏感，因此步长长h不宜太小不宜太小。例例5-1例题例题5-1 实验测得一组数据如下，求实验测得一组数据如下，求x=0.8处及各点处及各点的导数的导数 解解两种格式求得的微分值两种格式求得的微分值先求先求 x=0.8处处 导数导数 df/dxx00.10.20.30.40.50.60.70.80.91.0f(x)0.48 0.38 0.31 0.33 0.36 0.41 0.51 0.43 0.35 0.29 0.2806.01.035.029.08.09.0)8.0()9.0(8.0 ffdxdfx 二阶精度格式二阶精度格式)(O
8、2)()(211xxxfxfdxdfiixxi )(O)()(1xxxfxfdxdfiixxi 07.02.043.029.01.02)7.0()9.0(8.0 ffdxdfx 一阶精度格式一阶精度格式二阶格式二阶格式的精度更高。的精度更高。编程求导数。所给数据为等距节点，采用插值型编程求导数。所给数据为等距节点，采用插值型数值求导中的三点公式。数值求导中的三点公式。也可采用也可采用MATLAB中的数值微分函数中的数值微分函数diff求导数求导数 参考程序见教材参考程序见教材p.53程序代码程序代码 clear x=0:0.1:1;y=0.48 0.38 0.31 0.33 0.36 0.41
9、 0.51 0.43 0.35 0.29 0.28;h=0.1;%等距节点步长等距节点步长 for i=2:(length(x)-1);dy(i)=(y(i+1)-y(i-1)/2/h;%第第2点至点至第第n-1点导数值点导数值 end dy(1)=(2*y(2)-1.5*y(1)-0.5*y(3)/h;%第第1点导数值点导数值 dy(length(x)=(1.5*y(length(y)-2*y(length(y)-1)+0.5*y(length(y)-2)/h;%最末点导数值最末点导数值 fprintf(n%s tt%s tt%s,No,f(x),dy/dx)fprintf(n-)for i
10、=1:length(x)fprintf(n%d tt%.2f tt%.2f,i,y(i),dy(i)end数值微分结数值微分结果显示（二果显示（二阶精度）阶精度）No.xf(x)df(x)/dx1 00.48-1.150 2 0.10.38-0.8503 0.20.31-0.2504 0.30.330.2505 0.40.360.4006 0.50.410.7507 0.60.510.1008 0.70.43-0.8009 0.80.35-0.70010 0.90.29-0.35011 1.00.28 0.150diff(),polyder(),fnder()例例 已知函数已知函数y=cos(
11、x)，计算在，计算在x=0:pi/20:pi下对应下对应函数值向量函数值向量y，根据这些离散数据计算各点处函数的，根据这些离散数据计算各点处函数的一阶导数值一阶导数值%有限差分法：用差分函数有限差分法：用差分函数diff()近似计算导近似计算导数，即数，即dy=diff(y)./diff(x)%多项式拟合方法：先用多项式拟合方法：先用polyfit()根据离散根据离散数据拟合得到多项式插值函数数据拟合得到多项式插值函数p，%再用再用polyder()计算计算p的导数的导数pp,然后用然后用polyval()计算计算pp在在x的导数值的导数值%三次样条插值方法：先用三次样条插值方法：先用csap
12、i()根据离散根据离散数据生成三次样条插值函数数据生成三次样条插值函数cs，%再用再用fnder()计算计算cs的导数的导数pp,然后用然后用fnval()计算计算pp在在xi的导数值的导数值%用函数用函数y=cos(x)生成离散数据生成离散数据 x=0:pi/20:pi;y=cos(x);dy=-sin(x)%准确值准确值%有限差分法有限差分法 dy1=diff(y)./diff(x)x1=x(1:length(x)-1);clear all clc%多项式拟合方法多项式拟合方法 p=polyfit(x,y,5);%进行多项式拟合，得进行多项式拟合，得到多项式到多项式p pp=polyder
13、(p);%对拟合得到的多项式对拟合得到的多项式求导，得到导函数求导，得到导函数pp dy2=polyval(pp,x)%计算导函数计算导函数pp在在x的导数值的导数值%三次样条插值方法三次样条插值方法 cs=csapi(x,y);%生成三次样条插值函生成三次样条插值函数数cs pp=fnder(cs);%计算三次样条插值函计算三次样条插值函数数cs的导函数的导函数pp dy3=fnval(pp,x)%计算导函数计算导函数pp在在x的的导数值导数值%多项式拟合方法和三次样条插值方法的结多项式拟合方法和三次样条插值方法的结果与准确值之差果与准确值之差 disp(多项式拟合方法和三次样条插值方法的多
14、项式拟合方法和三次样条插值方法的结果与准确值之差结果与准确值之差dy-dy2;dy-dy3：)disp(dy-dy2;dy-dy3)%绘制几种方法的导函数曲线绘制几种方法的导函数曲线 plot(x,dy,ko,x1,dy1,k.-,x,dy2,b-,x,dy3,k-)xlabel(x)ylabel(y)title(函数函数cos(x)的导函数曲线的导函数曲线)legend(准确值准确值,有限差分法有限差分法,多项式多项式拟合方法拟合方法,三次样条插值方法三次样条插值方法)多项式拟合求导（多项式拟合求导（5 5阶）和三次样条插值求导阶）和三次样条插值求导所得曲线重合所得曲线重合 有限差分得到曲线
15、较差（只有有限差分得到曲线较差（只有1 1阶精度）阶精度）中心差分格式（中心差分格式（2 2阶）效果较好，但编程麻烦阶）效果较好，但编程麻烦 注：注：多项式拟合多项式拟合阶数敏感阶数敏感 离散数据由余弦函数生成，不包括噪声离散数据由余弦函数生成，不包括噪声（误差），三次样条插值结果较好，对（误差），三次样条插值结果较好，对于真实实验数据，必须改用于真实实验数据，必须改用三次样条拟三次样条拟合合 例例 微分法进行动力学数据分析微分法进行动力学数据分析 反应物反应物A A在一等温间歇反应器中发生反应为在一等温间歇反应器中发生反应为：A A产物产物 测量得到反应器中不同时间下反应物测量得到反应器中不
16、同时间下反应物A A的的浓度浓度c cA A如下如下表所示。试根据表中数据确定其反应速表所示。试根据表中数据确定其反应速率方程率方程 解解 数学模型数学模型 function Cha4demo8%动力学数据动力学数据 t=0 20 40 60 120 180 300;CA=10 8 6 5 3 2 1;%用最小二乘样条拟合法计算微分用最小二乘样条拟合法计算微分dCA/dt knots=3;K=3;%三次三次B样条样条 sp=spap2(knots,K,t,CA);sp=spap2(newknt(sp),K,t,CA);pp=fnder(sp);%计算计算B样条函数的导函数样条函数的导函数 dC
17、Adt=fnval(pp,t);%计算计算t处的导函数处的导函数值值%绘制图形绘制图形 ti=linspace(t(1),t(end),200);CAi=fnval(sp,ti);plot(t,CA,ro,ti,CAi,b-)，xlabel(t)，ylabel(C_A)figure fnplt(pp);%dCAdti=fnval(pp,ti)%plot(ti,dCAdti,-)xlabel(t)ylabel(dC/dt)%线性拟合线性拟合 rA=dCAdt;y=log(-rA);x=log(CA);p=polyfit(x,y,1);k=exp(p(2)n=p(1)化工常见数值积分化工常见数值积
18、分 双组分简单精馏（双组分简单精馏（LewisLewis法）法）精馏段理论塔板数精馏段理论塔板数 提馏段理论塔板数提馏段理论塔板数5.2 数值积分基础数值积分基础dfxxdRyxy-xxN/)(dfwxxwRxyy-xxN/)(d 填料吸收塔填料吸收塔 基于气相的总传质单元数（基于气相的总传质单元数（NTUNTU）为）为12yyOGyydyN)(*5.2 数值积分基础数值积分基础 定积分计算定积分计算 babaaFbFxFdxxf)()()()(4142444214 xabbaxbaxeeedxe 213dxx 并不是所有的积分都能求原函数和并不是所有的积分都能求原函数和/或积分值或积分值，原
19、因有：原因有：被积函数形式复杂，找不到积分公式求原函数被积函数形式复杂，找不到积分公式求原函数 原函数形式复杂，不便计算（如无穷级数）原函数形式复杂，不便计算（如无穷级数）被积函数以列表函数形式给出，无法用积分公被积函数以列表函数形式给出，无法用积分公式式问题：问题：能否不通过求原函数来求积分值能否不通过求原函数来求积分值?方法：方法：图解积分，图解积分，数值计算求积分数值计算求积分 数值计算求积分数值计算求积分构造一个比较简单的函数来构造一个比较简单的函数来近似被积函数（与插值法类似），然后通过计算简近似被积函数（与插值法类似），然后通过计算简单函数在区间上的积分值作为的积分近似值。单函数在
20、区间上的积分值作为的积分近似值。选择什么样的简单函数？选择什么样的简单函数？拉格朗日插值多项式拉格朗日插值多项式Ln(x):容易积分且能够近似大多数函数，通常被用来近容易积分且能够近似大多数函数，通常被用来近似被积函数似被积函数f(x)，称为插值型求积公式，称为插值型求积公式 对对Ln(x)积分就可求到的积分近似值。积分就可求到的积分近似值。下面讨论插值多项式下面讨论插值多项式Ln(x)中中n=1和和n=2时的积时的积分公式分公式()d()dbbnaaf xxL xx对插值多项式对插值多项式Ln(x)积分，积分，n不同时，可有不同的积不同时，可有不同的积分公式分公式 5.2.1 梯形公式（梯形
21、公式（n=1）线性插值公式线性插值公式)()()(1bfabaxafbaaxxL 1()d()d()()d()()2bbaabaf xxL xxxaxaf af bxabbabaf af b积分积分几何意义几何意义 几何意义几何意义：直线代替曲线，斜边梯形近似曲边梯形。：直线代替曲线，斜边梯形近似曲边梯形。)()(2)(2)(2bfafabbfabafabT 梯形公式梯形公式5.2.2 辛普森（辛普森（Simpson）公式（）公式（n=2）n=2时，插值多项式时，插值多项式L2(x)是二次插值式，在求积是二次插值式，在求积区间两端点，的中间加上一个中点区间两端点，的中间加上一个中点2/)(ab
22、c)()()()()()()()()()(2cbabcxaxbfbcacbxaxcfbacabxcxafxL 2()d()d()()()()()()()()()d()()()()()()()4()()6bbaabaf xxL xxxc xbxa xbxa xcf af cf cxac abca cbba bbbaf af cf b积分积分辛普生辛普生积分积分(三点公式三点公式)()(4)(6bfcfafabS cf(c)af(a)yf(b)b badxxLI)(2f(x)x 几何意义几何意义：用三点构造的抛物线：用三点构造的抛物线L2(x)代替被积函代替被积函数数f(x)构成曲边梯形面积近似原
23、曲边梯形面积。构成曲边梯形面积近似原曲边梯形面积。5.2.3 牛顿牛顿-科特斯（科特斯（Newton-Cotes）公式）公式（n等分）等分）积分得到五点公式积分得到五点公式柯特斯公式柯特斯公式对区间对区间n等分，过等分，过n+1个点时的拉格朗日个点时的拉格朗日(Lagrange)n次插值多项式为次插值多项式为 nknkjjkjkjnnkjyxxxxxL00,1,0,)()(牛顿牛顿-科特斯（科特斯（Newton-Cotes）公式）公式。使用该公式计算的积分值，使用该公式计算的积分值，关键是要求关键是要求Ak。00000()d()d()()d()d()()nnbbbjnkaaakjkjjknnn
24、bjkkkakkjkjjkxxfxxLxxfxxxxxxxfxAfxxx 求求Ak 积分后得积分后得由科特斯系数公式，给定由科特斯系数公式，给定n，就可计算，就可计算Cnk0dnbjkajkjj kxxAxxx00()d()()()nnbkknkkakkIf xxA f xbaCf x00(1)()d!()!n knnnkjj kCtjtnknkCnk科特斯（科特斯（Cotes）系数）系数7()32()12()32()7()90baCf af df cf ef b积分近似值的误差情况？积分近似值的误差情况？无论无论2点，点，3点，还是点，还是5点，实际上与准确值都有相点，实际上与准确值都有相当
25、的误差当的误差表表5-1列出了最常用的列出了最常用的n=16的科特斯（的科特斯（Cotes）系数，）系数，表中系数要乘第表中系数要乘第2列中的列中的K后才是后才是Cnk。如果如果n=4，即，即5个点，可直接写出牛顿个点，可直接写出牛顿-科特斯公式为科特斯公式为 5.3 复化求积公式复化求积公式实际计算时，为提高计算精度，通常是将积分区间实际计算时，为提高计算精度，通常是将积分区间划分为若干个小区间，对每个小区间分别求积分，划分为若干个小区间，对每个小区间分别求积分，再将积分结果累加得到整个区间的积分值再将积分结果累加得到整个区间的积分值复化复化求积方法求积方法 5.3.1 复化梯形公式复化梯形
26、公式 将积分区间将积分区间(a,b)划分为划分为n等分，步长等分，步长:h=(b-a)/n，kxakh0,1,kn将全部小区间积分值累加，得将全部小区间积分值累加，得复化梯形公式复化梯形公式111()d()()2kkxkkkkxxxf xxf xf x对小区间运用梯形公式求积分对小区间运用梯形公式求积分)()(2)(2)()(2111101bfxfafhxfxfxxTnkkknkkkkn5.3.2 复化辛普生（复化辛普生（Simpson）公式）公式辛普森公式要用到辛普森公式要用到3个节点，所以区间的等分数个节点，所以区间的等分数n应应为偶数，节点个数为奇数。为计算方便，对为偶数，节点个数为奇数
27、。为计算方便，对3个节点个节点2等分的小区间运用辛普森公式求积分等分的小区间运用辛普森公式求积分)()(4)(6d)(2122kkkkkxxxfxfxfxxxxfkk n等分的复化辛普森公式等分的复化辛普森公式1,3,124220212)(4)(2)()(3)()(4)(6nkknkknkkkkkknxfxfbfafhxfxfxfxxS，5.3.3 变步长辛普生求积变步长辛普生求积 定步长求积分定步长求积分积分区间分成若干个小区间，积分区间分成若干个小区间，区间数一定（每个区间长度即步长一定）区间数一定（每个区间长度即步长一定）给定步长给定步长（区间个数）（区间个数）复化求积公式复化求积公式计
28、算积分值计算积分值变步长求积分变步长求积分积分过程中将步长逐次对积分过程中将步长逐次对分，反复利用复化求积公式进行计算分，反复利用复化求积公式进行计算计算计算a,b、区间积区间积分值分值步长减半步长减半(区区间数加倍间数加倍)后后计算积分值计算积分值两次积分两次积分值之差是值之差是否否精度？精度？结束结束计算计算TF5.4 龙贝格（龙贝格（Romberg)积分法积分法 梯形公式算法简单，但精度低梯形公式算法简单，但精度低,误差大，能否利用误差大，能否利用简单的算法得到高精度的结果？简单的算法得到高精度的结果？可证明，梯形法积分值可证明，梯形法积分值Tn的截断误差近似与步长的截断误差近似与步长
29、h2 成正比。因此，当步长二分后，截断误差减至二分成正比。因此，当步长二分后，截断误差减至二分前的前的1/4，即：，即：41*2 nnTITI整理得整理得 nnnTTTI 2231*T2n的误差的误差 由上式，由上式，T2n 的误差约为的误差约为 1/3(T2n-Tn)，如将这个误，如将这个误差值作为差值作为T2n的补偿，应能得到精度更高的结果。的补偿，应能得到精度更高的结果。将梯形公式二分前后的积分值将梯形公式二分前后的积分值 Tn 和和T2n 按上式组合，按上式组合，可证明得到的结果与辛普生公式相同。可证明得到的结果与辛普生公式相同。nnnnnTTTTTT313431222 5.4 龙贝格
30、（龙贝格（Romberg)积分法积分法 通过梯形公式的线性组合得到的辛普生积分结果：通过梯形公式的线性组合得到的辛普生积分结果：辛普生公式积分值辛普生公式积分值 Sn 的截断误差近似与的截断误差近似与 h4 成正成正比比(可证明可证明)，步长二分后，截断误差减至原有误差，步长二分后，截断误差减至原有误差的的 1/16，有：，有：nnnTTST31342 161*2 nnSISI整理得整理得nnSSI1511516*2 将辛普生公式二分前后的积分值将辛普生公式二分前后的积分值 Sn 和和 S2n 按上式按上式线性组合，所得到的结果为柯特斯公式的积分值，线性组合，所得到的结果为柯特斯公式的积分值，
31、即：即：nnnSSC15115162 柯特斯公式积分值的精度高于辛普生公式。柯特斯公式积分值的精度高于辛普生公式。5.4 龙贝格（龙贝格（Romberg)积分法积分法nnnCCR63163642 龙贝格求积公式龙贝格求积公式 用梯形公式逐次减半步长得到的积分值序列：用梯形公式逐次减半步长得到的积分值序列：同理，柯特斯公式的误差与同理，柯特斯公式的误差与 h6 成正比，通过柯特成正比，通过柯特斯公式的线性组合同样可达到精度更高的结果：斯公式的线性组合同样可达到精度更高的结果：144 1iiiTTS212441iiiSSC313441iiiCCR 龙贝格求积法龙贝格求积法从单一梯形公式出发，逐次从
32、单一梯形公式出发，逐次减半步长，运用误差估计对近似值进行修正，得减半步长，运用误差估计对近似值进行修正，得到加快计算速度，提高精度的效果。到加快计算速度，提高精度的效果。优点：算法简单，精度高，速度快，计算量小。优点：算法简单，精度高，速度快，计算量小。收敛判断：两个相邻的积分近似值之差满足精收敛判断：两个相邻的积分近似值之差满足精度要求度要求。MATLAB 龙贝格积分函数龙贝格积分函数 quadl例题例题5-2：分别用定步长辛普生积分、变步长辛普生分别用定步长辛普生积分、变步长辛普生积分、龙贝格求积法求定积分的积分、龙贝格求积法求定积分的 值，定步长值，定步长求积时取求积时取N=6。解解 先
33、用梯形公式先用梯形公式(2点点)和辛普生公式和辛普生公式(3点点)计算计算10dxe x )()(2bfafabT 01201ee 85914.1171828.221 用梯形公式（取两点）：用梯形公式（取两点）：计算程序见教材计算程序见教材p.58。辛普生积分和龙贝格积分。辛普生积分和龙贝格积分直接调用函数。直接调用函数。计算结果：计算结果：S=1.718282 用辛普生公式（取三点）用辛普生公式（取三点）)()(4)(6bfcfafabS 05.014601eee 71886.1164872.1471828.261%例例5-2(1)a=0;b=1;n=7;h=(b-a)/n;%定积分定积分上
34、下限和步长上下限和步长 S=(b-a)*(exp(a)-exp(b)/2/(3*n);for i=1:n S=S+(2*exp(a+(i-1/2)*h)+exp(a+i*h)*(b-a)/(3*n);%复化辛普生求积复化辛普生求积 end fprintf(tS=%.6f,S)%输出结果输出结果 clear S=quad(exp,0,1);%调用辛普生积调用辛普生积分函数分函数 fprintf(tS=%.6f,S)clear S=quadl(exp,0,1);%调用龙贝格积调用龙贝格积分函数分函数 fprintf(tS=%.6f,S)例题例题 5-3 (选自选自化工原理化工原理)用 水 吸 收
35、混 合 气 中 的 丙 酮。进 气 丙 酮 含 量用 水 吸 收 混 合 气 中 的 丙 酮。进 气 丙 酮 含 量0.0 5(m o l%)，塔 内 吸 收，塔 内 吸 收 9 8%，吸 收 系 数，吸 收 系 数Kv=1.111110-4，常压和，常压和20时气体流量时气体流量Vt=0.5556 m3/s。相平衡关系式为。相平衡关系式为Y=MX，其中，其中M=1.68。求：。求：(1）出塔吸收液的最大浓度和吸收剂最小用量；）出塔吸收液的最大浓度和吸收剂最小用量；（2）当吸收液的浓度为最大浓度的）当吸收液的浓度为最大浓度的60时的吸收时的吸收剂用量；剂用量；（3）若填料比表面积）若填料比表面
36、积 a0=190 m2/m3，空塔气速，空塔气速 u=0.8m/s，计算塔径和填料层高度。，计算塔径和填料层高度。y1=0.05(mol%)解：解：（1）主要计算公式）主要计算公式(a)计算计算 Y1，Y2，X2进气丙酮比摩尔分率进气丙酮比摩尔分率 Y1=5/95 出气丙酮比摩尔分率出气丙酮比摩尔分率 Y2=(1-0.98)Y1 进塔吸收剂丙酮浓度进塔吸收剂丙酮浓度 X2=0 x2=0 x1=？y2=？吸收吸收98%例题例题 5-3(本例选自本例选自化工原理化工原理)(b)计算计算X1m，Lmin，X1，L惰性气体流量惰性气体流量(kmol/s)V=273(1-0.05)/29322.4Vt吸
37、收速率吸收速率(kmol/s)GA=V(Y1-Y2)吸收液出口最大浓度吸收液出口最大浓度 Xlm=Y1/M (与与Y1平衡的液体浓度）平衡的液体浓度）吸收剂最小用量吸收剂最小用量(kmol/s)Lmin=GA/(X1m-X2)吸收液出口实际浓度吸收液出口实际浓度 X1=0.6 Xlm 吸收剂实际用量吸收剂实际用量(kmol/s)L=GA(X1-X2)y1=0.05(mol%)x2=0 x1=？y2=？吸收吸收98%Y=MXLmin=?L=？V=？常压和常压和20时气体流量时气体流量)/(uVDt4(c)计算塔径（计算塔径（m）例题例题 5-3(本例选自本例选自化工原理化工原理)计算方法计算方法
38、()dbaSf xxdYMXYYY 121此处自变量为此处自变量为Y，故应将，故应将X用用Y表示表示(d)计算传质单元数计算传质单元数 NOG数值积分法数值积分法F(Y)12*YYOGYYdYN22()V YYL XX塔的任一截面与塔顶间物料平衡塔的任一截面与塔顶间物料平衡1212()()V YYL XX全塔物料平衡式全塔物料平衡式212122()()()YYXXXYYX联立上两式求出联立上两式求出上式代入上式代入*1YYF2121221()()/()FYYYXXYYXM(e)计算传质单元高度和填料高度计算传质单元高度和填料高度HOG0OGvVHKa（m）OGOGZHN（m）未知量未知量X用用
39、Y表示表示X（2）计算程序及结果）计算程序及结果参考程序见教材参考程序见教材p.59。计算结果计算结果出塔吸收液最大浓度出塔吸收液最大浓度 Xmax=0.0313吸收剂最小用量吸收剂最小用量 Lmin=0.0361 kmol/s出塔吸收液实际浓度出塔吸收液实际浓度 x1=0.0188吸收剂实际用量吸收剂实际用量 L=0.0602 kmol/s填料塔径填料塔径 D=0.9404填料层高度填料层高度 Z=11.5691m function y=f(x,Y1,Y2,M,x1,x2)y=1./(x-(x-Y2).*(x1-x2)./(Y1-Y2)+x2).*M);%被积函数被积函数F的表达式的表达式
40、必须是向量形式必须是向量形式 clear KV=1.1111*10(-4);Vt=0.5556;M=1.68;a0=190;u=0.8;Y1=5/95;%入塔气体丙酮比摩尔分率入塔气体丙酮比摩尔分率 Y2=(1-0.98)*Y1;%出塔气丙酮比摩尔分出塔气丙酮比摩尔分率率 x2=0;%吸收剂入口浓度吸收剂入口浓度 V=273*(1-0.05)*Vt/293/22.4;%惰性惰性气体流量气体流量 GA=V*(Y1-Y2);%吸收速率吸收速率 X1m=Y1/M;%吸收液出口最大浓度吸收液出口最大浓度 Lmin=GA/(X1m-x2);%吸收剂最小用量吸收剂最小用量 x1=X1m*0.6;%吸收液出
41、口实际浓度吸收液出口实际浓度 L=GA/(x1-x2);%吸收剂实际用量吸收剂实际用量 D=sqrt(4*Vt/(pi*u);%塔径计算塔径计算 Nog=quad(f,Y2,Y1,Y1,Y2,M,x1,x2);%调用辛普生积分函数计算传质单元数调用辛普生积分函数计算传质单元数 Hog=V/(KV*a0*(pi*D2/4);%计算传计算传质单元高度质单元高度 Z=Nog*Hog;%计算填料层高度计算填料层高度 string=kmol/s;fprintf(tXmax=%.4ftLmin=%.4ft%sntx1=%.4fttL=%.4ftt%s,X1m,Lmin,string,x1,L,string
42、)fprintf(ntD=%.4fttZ=%.4ftt%s,D,Z,m)被积函数为列表函数时的积分问题如何计算？被积函数为列表函数时的积分问题如何计算？x1.01.11.21.31.41.51.61.71.8f(x)1.543 1.668 1.8111.9712.151 2.3522.5772.8283.107dxxf8.10.1)(积分：积分：步长分别取步长分别取(a)h=0.1,(b)h=0.2，(c)h=0.4。(a)h=0.1，用复化梯形公式，用复化梯形公式101)()(2nkkkxfxfhT)()(2)(211bfxfafhnkk7683.1)828.2577.2352.2151.2
43、971.1811.1668.1(2107.3543.121.0 T辛普生公式辛普生公式1312，4，2423nkknkknxfxfbfafhS,)()()()(7668.1)828.2352.2971.1668.1(4)577.2151.2811.1(2107.3543.131.0Sdxxf8.10.1)(积分：积分：x1.01.11.21.31.41.51.61.71.8f(x)1.543 1.6681.8111.9712.1512.3522.5772.8283.107 (b)h=0.2，（，（5个点）个点）7728.1)577.2151.2811.1(2107.3543.122.0 T辛普
44、生公式辛普生公式 76693.1)577.2811.1(4)151.2(2107.3543.132.0S (c)h=0.4,（3个点）个点）7904.1151.22107.3543.124.0 T梯形公式梯形公式 梯形公式梯形公式7672.1)151.2(4107.3543.134.0S辛普生公式辛普生公式 计算结果比较计算结果比较步长步长h=0.1h=0.2h=0.4梯形公式梯形公式1.76831.77281.7904辛普生公式辛普生公式1.76681.766931.7672柯特斯公式柯特斯公式1.76679 步长相同时，辛普生公式步长相同时，辛普生公式 与用柯特斯公式的计算与用柯特斯公式的
45、计算精度高于梯形公式精度高于梯形公式 例题例题 5-4 用脉冲法测量某反应器的停留时间分用脉冲法测量某反应器的停留时间分布，由实验数据计算寿命分布密度与时间布，由实验数据计算寿命分布密度与时间t的的关系如下表关系如下表t10152025303540E(t)103 01.929.020.9531.10 32.20 30.60求分布密度积分求分布密度积分4010()dE tt解：解：根据所给数据，小区间分为根据所给数据，小区间分为6等分，步长等分，步长h=5，节点个数为节点个数为7，用辛普生复化求积公式计算。，用辛普生复化求积公式计算。0)10()(faf31060.30)40()(fbf403105()d030.62(9.031.1)4(1.9220.9532.2)100.55183E tt1,3,1242)(4)(2)()(3nkknkknxfxfbfafhS，本章要求本章要求 清楚数值积分的原理清楚数值积分的原理 掌握数值积分的梯形公式和辛普生公式、复化掌握数值积分的梯形公式和辛普生公式、复化梯形公式和复化辛普生公式梯形公式和复化辛普生公式 了解龙贝格积分的主要思想了解龙贝格积分的主要思想 能用数值积分求解实际问题能用数值积分求解实际问题分别用定步长和变步长辛普森积分计算分别用定步长和变步长辛普森积分计算102d14xx800dcos.xx
