《数值微分与数值积分.ppt》由会员分享，可在线阅读，更多相关《数值微分与数值积分.ppt（47页珍藏版）》请在三一办公上搜索。
1、2023/6/28,1,第5章,数值微分与数值积分,2,2023/6/28,微积分是高等数学中的重要内容，在化学工程上有许多非常重要的应用微积分的数值方法，不同于高等数学中的解析方法，尤其适合求解没有或很难求出微分或积分表达式的实际化工问题的计算，例如：列表函数求微分或积分,引言1-数值微积分方法不同于解析方法,3,2023/6/28,数值微分和数值积分与插值和拟合往往是密不可分的如在进行数值微分时，常针对离散的数据点，利用插值和拟合可以减少误差；而数值积分的基本思路也来自于插值法。例如如果所积函数的形式比较复杂或以表格形式给出，则可通过构造一个插值多项式来代替原函数，从而使问题简化,引言2-
2、数值微分和数值积分与插值和拟合关系密切,4,2023/6/28,5.1 数值微分,化工领域的实际问题中时常需要求列表函数在节点和非节点处的导数值，这是数值微分所要解决的问题。数值微分方法可近似求出某点的导数值例如在反应动力学的研究中，根据实验数据确定反应的动力学方程：这里实验测得一批离散点，要计算 只能借助数值微分求导解决 表5-1 反应动力学实验数据,0 导入,5,2023/6/28,0.1 建立数值微分公式的三种思路,常用三种思路建立数值微分公式：从微分定义出发，通过近似处理，得到数值微分的近似公式从插值近似公式出发，对插值公式的近似求导可得到数值微分的近似公式先用最小二乘拟合方法根据已知
3、数据得近似函数，再对此近似函数求微分可得到数值微分的近似公式。然后对各方法数值微分后得到的多项式求值，即可求出任意点处的任意阶微分,6,2023/6/28,1 方法概述在微积分中，一阶微分的计算可以在二相邻点x+h和x间函数取下列极限求得：,取其达到极限前的形式，就得到以下微分的差分近似式：,注：高阶微分项可以利用低阶微分项来计算，如二阶微分式可以表示为：,对应的差分式有：,5.1.1 差分近似微分,上式中三种不同表示形式依次是一阶前向差分、一阶后向差分和一阶中心差分来近似表示微分。其中一阶中心差分的精度较高。,7,2023/6/28,2 差分的Matlab实现,在Matlab中，可用diff
4、函数进行离散数据的近似求导调用形式：Y=diff(X,n)其中：X表示求导变量，可以是向量或矩阵。如是矩阵形式则按各列作差分；n表示n阶差分，即差分n次 Y是X的差分结果注：用diff函数进行离散数据的近似求导与前向差分近似，但误差较大。最好将数据利用插值或拟合得到多项式，然后对近似多项式进行微分,8,2023/6/28,例5.1：丁二烯的气相二聚反应如下：,实验在一定容器的反应器中进行，3260C时，测得物系中丁二烯的分压,(mmHg)与时间的关系如表5-2所示。,。,表 5-2 丁二烯二聚反应实验数据,用数值微分法计算所列时刻每一瞬间的反应速率,9,2023/6/28,解：程序如下：t=0
5、: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=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.80
6、00-2.4000-2.2000-2.2000-2.0000-2.0000 Columns 17 through 18-2.0000-2.0000,10,2023/6/28,5.1.2 三次样条插值函数求微分,11,2023/6/28,1 方法概述用三次样条插值函数建立的数值微分公式为：,求导得：,（5-2）,;,式（5-1）和（5-2）不但适用于求节点处的导数，而且可求非节点处的导数。,（5-1）,其中，,12,2023/6/28,2 三次样条插值函数求微分的Matlab函数Matlab求离散数据的三次样条插值函数微分方法分三个步骤：Step 1:对离散数据用csapi函数（或spline函
7、数）得到其三次样条插值函数调用形式 pp=csapi(x,y)其中：x,y分别为离散数据对的自变量和因变量；pp为得到的三次样条插值函数。Step 2:用fnder函数求三次样条插值函数的导数调用形式 fprime=fnder(f,dorder)其中：f为三次样条插值函数；dorder为三次样条插值函数的求导阶数；fprime为得到的三次样条插值函数的导函数。Step3:可用fnval函数求导函数在未知点处的导数值调用形式 v=fnval(fprime,x)其中：fprime为三次样条插值函数导函数；x为未知点处自变量值；v为未知点处的导数值。,13,2023/6/28,例5.2：某液体冷却时
8、，温度随时间的变化数据如表5-3所示：表5-3 冷却温度随时间的变化数据试分别计算t2，3，4min及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.5
9、000 2.5000 4.5000-5.3722-4.6722-3.8389-5.7972-4.9889-3.2222,注：前者是计算节点处的一阶导数，后者是计算非节点处的一阶导数,14,2023/6/28,5.1.3 最小二乘法样条拟合函数求微分,在实际化工应用中，当来自实验观测的离散数据不可避免地含有较大随机误差时，此时用插值公式求数值微分虽然样本点处误差较小，但可能会使非样本点处产生较大误差为此，可采用最小二乘法样条拟合实验数据，获得一个函数模型，然后再对其求导数。与样条插值不同，样条拟合不要求曲线经过全部的数据点，这样处理求导结果会有很大改善,15,2023/6/28,1 方法概述 可
10、用最小二乘法拟合成平滑B样条曲线，即对于离散数据（,所求的K次样条拟合函数,满足：,其中：,再对拟合函数作平滑处理后求导，即可求出任意点处微分。,）,为权重系数，默认为1,16,2023/6/28,Matlab求离散数据的最小二乘法平滑B样条拟合函数求微分共三个步骤：Step 1:对离散数据用spaps函数得到最小二乘平滑B样条拟合函数。调用格式：sp=spaps(x,y,tol)其中：x,y要处理的离散数据（）tol光滑时的允许精度，通常取（10-210-4）Step 2:可用fnder函数求最小二乘平滑B样条拟合函数的导数;Step3:可用fnval函数求导函数在未知点处的导数值。fnde
11、r()和fnval（）调用形式前文中已经介绍过。,2 最小二乘法平滑B样条拟合函数求微分的Matlab实现,17,2023/6/28,由离散数据求数值微分的四种方法及有关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表示的多项式拟合函数,导函数pp,pp在x
12、i的导数值。,其中：函数polyfit()和polyval()在前文中已介绍导函数polyder()的调用格式为：pp=polyder(p),离散数据,该函数对向量p表示的多项式函数进行求导，返回导函数为pp。,18,2023/6/28,（3）三次样条插值方法（4）样条拟合方法（最小二乘法）其中：函数csaps()、spap2()、spaps()和fnval()已在前文中介绍，求导函数fnder()的调用格式为：pp=fnder(f,dorder)该函数计算函数f的dorder阶导函数，默认阶数dorder=1。,离散数据,三次样条插值函数cs,cs的导数pp,pp在xi的导数值,离散数据,样
13、条拟合函数sp,sp的导数pp,pp在xi的导数值。,19,2023/6/28,例5.3：反应物A在一等温间歇反应器中发生的反应为：A,测量得到的反应器中不同时间下反应物A的浓度,表5-4 间歇反应器动力学数据,系统的动力学模型为：,，试根据表中数据确定其反应速率方程。,产物,如表5-4所示。,20,2023/6/28,解：（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(s
14、p);%生成光滑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),及得到拟合曲线的图形程序如下：,21,2023/6/28,（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所以本例的反应速率方程为：,22,2023/6/28,1)被积函
15、数以一组数据形式表示；2)被积函数过于特殊或原函数不能用初等函数表示，积分表中无 法找到可沿用的现成公式；3)有的原函数十分复杂难以计算。,对于积分：,但是在工程技术和科学研究中,常会见到以下现象:,0 导入,5.2 数值积分,23,2023/6/28,数值积分就是一种常用的近似计算方法。数值积分方法不受以上几个问题的限制，在化学化工领域应用甚广，如反应热效应计算、热容计算、熵的计算、反应活化能的计算等。,以上这些现象,Newton-Leibniz很难发挥作用，只能建立积分的近似计算方法,24,2023/6/28,如：,某反应器中进行的反应，计算出口物料的焓值：,25,2023/6/28,0.
16、1 数值积分的基本思路和方法,常用的数值积分的基本思路来自于插值法，它通过构造一个插值多项式 作为 的近似表达式，用 的积分值作为 的近似积分值。数值积分的方法很丰富，常用的插值型求积公式有两类：一类是等距节点的牛顿柯特斯求积公式；另一类是不等距节点的高斯型求积公式。,26,2023/6/28,Newton-Cotes公式是指等距节点下使用Lagrange插值 多项式建立的数值求积公式,各节点为,则：,牛顿柯特斯（Newton-Cotes）求积公式1 方法概述,27,2023/6/28,这里 是既不依赖于被积函数，也不依赖于积分区间的常数，称为柯特斯系数。式（5-3）称为牛顿柯特斯求积公式。,
17、其中：,一般地：令，,得：,（5-3）,28,2023/6/28,在Newton-Cotes公式中,n=1,2,4时的公式是最常用也是最重要三个公式,称为低阶公式（当n=0时的公式为矩形公式）,1）梯形(trapezia)公式,Cotes系数为,求积公式为,29,2023/6/28,上式称为梯形求积公式,（5-4）,30,2023/6/28,2）Simpson公式,Cotes系数为,31,2023/6/28,上式称为Simpson求积公式，也称三点公式或抛物线公式,求积公式为,式（5-4）和式（5-5）是化工领域常用的两个求积公式。与梯形法求积公式相比，Simpson法求积公式是一个较高精度的
18、求积公式。用式（5-3）还可得到更高阶数,（5-5）,32,2023/6/28,Newton-Cotes公式当n大于7时，公式的稳定性将无法保证，因此,在实际应用中一般不使用高阶Newton-Cotes公式而是采用低阶复合求积法在实际计算中为了保证计算的精度，往往首先用分点xk=a+kh,(k=0,1,n)将区间a,b分成n个相等的子区间，而后对每个子区 间再应用梯形公式或Simpson公式，分别得到：,2 复化法求积公式,复化Simpson公式：,33,2023/6/28,尽管复化法求积公式具有很高的精度，但是它必须采用等步长方法，从而限制了它的效率。这里我们介绍一种更加灵活选取步长的方法，
19、即自适应步长法。,（5-8）,。,3 自适应求积公式,34,2023/6/28,35,2023/6/28,4 Newton-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。,36,2023/
20、6/28,2）自适应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函数更有效，精度更高。
21、其使用方法和quad（）完全相同。,需要了解更多的内容可查阅Matlab功能函数库（funfun）。,37,2023/6/28,例5.4：真实气体的逸度,可用下式计算：,现测得00C下氢气的有关数值如表5-5所示，试求1000 atm下的逸度。,表5-5 00C下氢气的相关数值,-真实气体的实测体积和按理想气体定律计算得到的体积之间的差值。,38,2023/6/28,解：本题是离散型数据，可用trapz函数求解数值积分。（1）计算数值积分，程序如下：%计算数值积分P=0:100:1000;a1=15.46 15.46 15.46 15.61 15.85 15.93 16.09 16.13 16
22、.16 16.16 16.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.00 atm,39,2023/6/28,例5.5：氯仿-苯双组分精馏系统的气液平衡数据如表5-6所示。规定进料和塔顶的组成分别是,，精馏段的回流比为,精馏段理论板数的模型为,表5-6 精馏段气液平衡数据,，,试用Matlab计算所需的精馏段理论板数。,40,2023/6/28
23、,解：因模型中的的函数关系是以表格形式给出，我们若要用辛普森法等计算较精确的精馏段理论板数的时候，需先采用拟合法将离散数据（，）拟合成多项式，再将多项式代入被积函数求积。计算程序如下：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);plo
24、t(xi,yi,o,xplot,yplot,-);N=quad(func1,0.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块。,41,2023/6/28,5.2.2 高斯勒让德公式,高斯法是一种精度较高的求积分法，它的一般公式是式中(x)是一个权重函数，Aj为系数，xj为横坐标上的节点 高斯勒让德多项式的权重函数(x)=1，因而，一个n点高斯勒让德求积公式具有如下形式：,1 方法概述,42,2023/6/
25、28,右边的f(xj)是函数f(x)在节点xj处的值，节点xj是勒让德多项式Pn(x)的根。Aj 为系数，其值为,式中Pn(x)是勒让德多项式Pn(x)的一阶导数,43,2023/6/28,pn(x)的前几项表达式为,由高斯勒让德多项式得出的26点的根和系数见表5-7,44,2023/6/28,表5-7 26点的高斯节点和高斯求积系数,45,2023/6/28,区间a,b内的高斯勒让德求积公式 高斯勒让德求积公式也可转换成求区间a,b内的积分公式，得到如下结果：,补充,46,2023/6/28,2 高斯勒让德法的Matlab 实现,Matlab没有提供高斯勒让德法的有关计算函数，可由其数学原理
26、编程实现。参考程序如下：function q=gaussL(f,a,b,x,A)N=length(x);T=zeros(1,N);T=(a+b)/2+(b-a)/2)*x;q=(b-a)/2)*sum(A.*feval(f,T);程序说明：f被积函数，可编写m文件实现；a,b 分别是积分的上下限；q-所求积分值；x 和A 分别由查节点与系数表5-7得到。,47,2023/6/28,例5.6：等压过程中使乙炔体系温度由 加热到 所需的热量 可按下式计算：用3点高斯勒让德法计算从 加热到 所需的热量。解：（1）用M文件定义名为fun2.m的函数：function Q=fun2(t)Q=44.16-0.047.*t-0.00002.*t.2(2)查高斯节点和高斯求积系数表可知：x=-0.7746,0.7746,0 A=0.55556,0.55556,0.88889计算主程序为：function jiaocaili56QgaussL(fun2,25,100,-0.7746,0.7746,0,0.55556,0.55556,0.88889)执行结果：Q=3.0851e+003所以，等压过程中使乙炔体系温度从加热到所需的热量为3085.1焦耳。,
