数值微分与数值积分课件

1、2022/9/261第5章 数值微分与数值积分2022/9/241第5章 数值微分与数值积分22022/9/26微积分是高等数学中的重要内容，在化学工程上有许多非常重要的应用微积分的数值方法，不同于高等数学中的解析方法，尤其适合求解没有或很难求出微分或积分表达式的实际化工问题的计算，例如：列表函数求微分或积分 引言1-数值微积分方法不同于解析方法22022/9/24微积分是高等数学中的重要内容，在化学工程32022/9/26数值微分和数值积分与插值和拟合往往是密不可分的如在进行数值微分时，常针对离散的数据点，利用插值和拟合可以减少误差；而数值积分的基本思路也来自于插值法。例如如果所积函数的形式

2、比较复杂或以表格形式给出，则可通过构造一个插值多项式来代替原函数，从而使问题简化引言2-数值微分和数值积分与插值和拟合关系密切32022/9/24数值微分和数值积分与插值和拟合往往是密不42022/9/265.1 数值微分化工领域的实际问题中时常需要求列表函数在节点和非节点处的导数值，这是数值微分所要解决的问题。数值微分方法可近似求出某点的导数值例如在反应动力学的研究中，根据实验数据确定反应的动力学方程 ：这里实验测得一批离散点，要计算 只能借助数值微分求导解决 表5-1 反应动力学实验数据 t1t2tnpA1pA2PAn0 导入42022/9/245.1 数值微分化工领域的实际问题52022

3、/9/260.1 建立数值微分公式的三种思路 常用三种思路建立数值微分公式：从微分定义出发，通过近似处理，得到数值微分的近似公式从插值近似公式出发，对插值公式的近似求导可得到数值微分的近似公式先用最小二乘拟合方法根据已知数据得近似函数，再对此近似函数求微分可得到数值微分的近似公式。然后对各方法数值微分后得到的多项式求值，即可求出任意点处的任意阶微分52022/9/240.1 建立数值微分公式的三种思路 62022/9/261 方法概述在微积分中，一阶微分的计算可以在二相邻点x+h和x间函数取下列极限求得： 取其达到极限前的形式，就得到以下微分的差分近似式： 注：高阶微分项可以利用低阶微分项来计

4、算，如二阶微分式可以表示为： 对应的差分式有：5.1.1 差分近似微分上式中三种不同表示形式依次是一阶前向差分、一阶后向差分和一阶中心差分来近似表示微分。其中一阶中心差分的精度较高。62022/9/241 方法概述取其达到极限前的形式，就得到72022/9/262 差分的Matlab实现在Matlab中，可用diff函数进行离散数据的近似求导调用形式：Y = diff(X,n)其中：X表示求导变量，可以是向量或矩阵。如是矩阵形式则按各列作差分； n表示n阶差分，即差分n次 Y是X的差分结果注：用diff函数进行离散数据的近似求导与前向差分近似， 但误差较大。最好将数据利用插值或拟合得到多项式，

5、然后对近似多项式进行微分72022/9/242 差分的Matlab实现在Matla82022/9/26例5.1：丁二烯的气相二聚反应如下：实验在一定容器的反应器中进行，3260C时，测得物系中丁二烯的分压 (mmHg)与时间的关系如表5-2所示。表 5-2 丁二烯二聚反应实验数据t(min) (mmHg)t(min) (mmHg)0632.050362.05590.055348.010552.060336.015515.065325.020485.070314.025458.075304.030435.080294.035414.085284.040396.090274.045378.0用数值

6、微分法计算所列时刻每一瞬间的反应速率 82022/9/24例5.1：丁二烯的气相二聚反应如下：实验92022/9/26解：程序如下：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=dpA./dt q为数值微分结果执行结果：q = Columns 1 through 8 -8.4000 -7.6000 -7.400

7、0 -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.000092022/9/24解：程序如下：102022/9/265.1.2 三次样条插值函数求微分若三次样条插值函数收敛于，那么导数收敛于，因此用样条插值函数作为函数，不但彼此的函数值非常接近，而且导数值也很接近。 用三次样条插值函数求数值导数是可靠的，这是化工计算中求数值微分的

8、有效方法。的近似102022/9/245.1.2 三次样条插值函数求微112022/9/261 方法概述用三次样条插值函数建立的数值微分公式为： 求导得： （5-2）;式（5-1）和（5-2）不但适用于求节点处的导数，而且可求非节点处的导数。（5-1）其中，112022/9/241 方法概述求导得： （5-2）;式122022/9/262 三次样条插值函数求微分的Matlab函数Matlab求离散数据的三次样条插值函数微分方法分三个步骤：Step 1:对离散数据用csapi函数（或spline函数）得到其三次样条插值函数调用形式 pp = csapi(x,y)其中： x,y分别为离散数据对的自

9、变量和因变量； pp为得到的三次样条插值函数。Step 2: 用fnder函数求三次样条插值函数的导数调用形式 fprime = fnder(f,dorder)其中： f为三次样条插值函数； dorder为三次样条插值函数的求导阶数； fprime为得到的三次样条插值函数的导函数。Step3:可用fnval函数求导函数在未知点处的导数值调用形式 v = fnval(fprime,x)其中：fprime为三次样条插值函数导函数； x为未知点处自变量值； v为未知点处的导数值。122022/9/242 三次样条插值函数求微分的Matl132022/9/26例5.2 ：某液体冷却时，温度随时间的变化

10、数据如表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.50

11、00 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注：前者是计算节点处的一阶导数，后者是计算非节点处的一阶导数132022/9/24例5.2 ：某液体冷却时，温度随时间142022/9/265.1.3 最小二乘法样条拟合函数求微分在实际化工应用中，当来自实验观测的离散数据不可避免地含有较大随机误差时，此时用插值公式求数值微分虽然样本点处误差较小，但可能会使非样本点处产生较大误差为此，可采用最小二乘法样条拟合实验数据，获得一个函数模型，

12、然后再对其求导数。与样条插值不同，样条拟合不要求曲线经过全部的数据点，这样处理求导结果会有很大改善142022/9/245.1.3 最小二乘法样条拟合函数求152022/9/261 方法概述 可用最小二乘法拟合成平滑B样条曲线，即对于离散数据（所求的K次样条拟合函数满足：其中：再对拟合函数作平滑处理后求导，即可求出任意点处微分。）为权重系数，默认为1152022/9/241 方法概述所求的K次样条拟合函数满足162022/9/26Matlab求离散数据的最小二乘法平滑B样条拟合函数求微分共三个步骤：Step 1:对离散数据用spaps函数得到最小二乘平滑B样条拟合函数。调用格式：sp = sp

13、aps(x,y,tol)其中： x,y要处理的离散数据（） tol光滑时的允许精度，通常取（10-210-4）Step 2:可用fnder函数求最小二乘平滑B样条拟合函数的导数;Step3: 可用fnval函数求导函数在未知点处的导数值。 fnder( )和fnval（）调用形式前文中已经介绍过。 2 最小二乘法平滑B样条拟合函数求微分的Matlab实现162022/9/24Matlab求离散数据的最小二乘法平滑172022/9/26由离散数据求数值微分的四种方法及有关MATLAB函数：（1）差分法 用差分函数diff()近似计算导数，即dy=diff(y)./diff(x)。对于向量X，di

14、ff(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表示的多项式拟合函数导函数pppp在xi的导数值。其中：函数polyfit( )和polyval( )在前文中已介绍导函数polyder( )的调用格式为：pp=polyder(p) 离散数据 该函数对向量p表示的多项式函数进行求导，返回导函数为pp。172022/9/24由离散数据求数值微分的四种方法及有关M182022/9/26（

15、3）三次样条插值方法 （4）样条拟合方法（最小二乘法） 其中：函数csaps()、spap2()、spaps()和fnval()已在前文中介绍， 求导函数fnder()的调用格式为： pp=fnder(f,dorder) 该函数计算函数f的dorder阶导函数，默认阶数dorder=1。离散数据三次样条插值函数cscs的导数pppp在xi的导数值 离散数据样条拟合函数spsp的导数pppp在xi的导数值。 182022/9/24（3）三次样条插值方法 离散数据三次192022/9/26例5.3：反应物A在一等温间歇反应器中发生的反应为：A测量得到的反应器中不同时间下反应物A的浓度表5-4 间歇

16、反应器动力学数据t(s)0204060120180300 (mol/L)10865321系统的动力学模型为：，试根据表中数据确定其反应速率方程。产物如表5-4所示。192022/9/24例5.3：反应物A在一等温间歇反应器中202022/9/26解：（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);

17、%计算导函数在t处的数值微分% 绘制图形ti=linspace(t(1),t(end),200);CAi=fnval(sp,ti);plot(t,CA,bo,ti,CAi,r-),xlabel(t),ylabel(CA)及得到拟合曲线的图形程序如下：202022/9/24解：（1）系统的动力学模型为非线性形式212022/9/26（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所以本例的反应速率方程为：212022/9/24（3）线性拟合程序如下：22

18、2022/9/261)被积函数以一组数据形式表示；2)被积函数过于特殊或原函数不能用初等函数表示，积分表中无 法找到可沿用的现成公式；3)有的原函数十分复杂难以计算。对于积分：但是在工程技术和科学研究中,常会见到以下现象:0 导入5.2 数值积分222022/9/241)被积函数以一组数据形式表示；对于积232022/9/26 数值积分就是一种常用的近似计算方法。数值积分方法不受以上几个问题的限制，在化学化工领域应用甚广，如反应热效应计算、热容计算、熵的计算、反应活化能的计算等。 以上这些现象,Newton-Leibniz很难发挥作用，只能建立积分的近似计算方法232022/9/24 数值积分

19、就是一种常用的近似计算方242022/9/26如：某反应器中进行的反应，计算出口物料的焓值：242022/9/24如：某反应器中进行的反应，计算出口物料么么么么方面Sds绝对是假的么么么么方面Sds绝对是假的262022/9/260.1 数值积分的基本思路和方法 常用的数值积分的基本思路来自于插值法，它通过构造一个插值多项式 作为 的近似表达式，用 的积分值作为 的近似积分值。数值积分的方法很丰富，常用的插值型求积公式有两类：一类是等距节点的牛顿柯特斯求积公式； 另一类是不等距节点的高斯型求积公式。262022/9/240.1 数值积分的基本思路和方法 272022/9/26 Newton-C

20、otes公式是指等距节点下使用Lagrange插值 多项式建立的数值求积公式各节点为则：5.2.1牛顿柯特斯（Newton-Cotes）求积公式1 方法概述272022/9/24 Newton-Cotes公式是指282022/9/26这里 是既不依赖于被积函数，也不依赖于积分区间的常数，称为柯特斯系数。式（5-3）称为牛顿柯特斯求积公式。其中： 一般地：令 ，得： （5-3）282022/9/24这里 是既不依赖于被积函292022/9/26在Newton-Cotes公式中,n=1,2,4时的公式是最常用也是最重要三个公式,称为低阶公式（当n=0时的公式为矩形公式）1）梯形(trapezia)

21、公式Cotes系数为求积公式为292022/9/24在Newton-Cotes公式中,n=302022/9/26上式称为梯形求积公式,（5-4）302022/9/24上式称为梯形求积公式,（5-4）312022/9/262）Simpson公式Cotes系数为312022/9/242）Simpson公式Cotes系数为322022/9/26上式称为Simpson求积公式，也称三点公式或抛物线公式求积公式为式（5-4）和式（5-5）是化工领域常用的两个求积公式。与梯形法求积公式相比，Simpson法求积公式是一个较高精度的求积公式。用式（5-3）还可得到更高阶数（5-5）322022/9/24上式

22、称为Simpson求积公式，也称三332022/9/26 Newton-Cotes公式当n大于7时，公式的稳定性将无法保证，因此,在实际应用中一般不使用高阶Newton-Cotes公式而是采用低阶复合求积法在实际计算中为了保证计算的精度，往往首先用分点xk=a+kh, (k=0,1,n)将区间a,b分成n个相等的子区间，而后对每个子区 间再应用梯形公式或Simpson公式，分别得到： 2 复化法求积公式复化Simpson公式：332022/9/24 Newton-Cotes公式当n大于342022/9/26 尽管复化法求积公式具有很高的精度，但是它必须采用等步长方法，从而限制了它的效率。这里我

23、们介绍一种更加灵活选取步长的方法，即自适应步长法。（5-8）。（5-6） （5-7），令，当 上Simpson积分 达到精度时，可认为区间取计算需满足的精度为 以Simpson积分法为例，某区间 ，记 ，考虑该区间上的Simpson积分和二等分以后的两个Simpson积分和：3 自适应求积公式342022/9/24 尽管复化法求积公式具有很高的精352022/9/26 这样重复下去，直至每个分段部分达到相应精度（步长为时精度为 ）。这样，不同段的步长可能是不一样的，积分结果为每一小段的积分总和。 自适应步长Simpson法从 开始按式（5-6）（5-8）的方法检验。若满足精度，则以 为计算结果

24、；若不满足精度，则分成两个小区间各自重复逐步上述过程，每个小区间精度为352022/9/24 这样重复下去，直362022/9/264 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被等

25、分，每份宽为1。362022/9/244 Newton-Cotes求积公372022/9/262） 自适应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法数值积分：quad

26、l（） quadl是高阶的自适应Newton-Cotes数值积分法函数，它比quad函数更有效，精度更高。其使用方法和quad（）完全相同。需要了解更多的内容可查阅Matlab功能函数库（funfun）。372022/9/242） 自适应Simpson法数值积分：382022/9/26例5.4：真实气体的逸度可用下式计算： 现测得00C下氢气的有关数值如表5-5所示，试求1000 atm下的逸度。表5-5 00C下氢气的相关数值 （atm）（atm）015.4660053.4316.09100239.5115.4670048.1416.13200127.4915.4680044.1716.16

27、30090.2915.6190041.0616.1640071.8615.85100038.5516.1450060.7615.93-真实气体的实测体积和按理想气体定律计算得到的体积之间的差值。R-气体常数 其中： -逸度； -压力，atm；T-绝对温度，K。 382022/9/24例5.4：真实气体的逸度可用下式计算：392022/9/26解：本题是离散型数据，可用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.

28、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 atm392022/9/24解：本题是离散型数据，可用trapz函402022/9/26例5.5：氯仿-苯双组分精馏系统的气液平衡数据如表5-6所示。规定进料和塔顶的组成分别是 ，精馏段的回流比为精馏段理论板数的模型为表5-6 精馏段气液平衡数据0.1780.2750.3720.4

29、560.6500.8440.2430.3820.5180.6160.7950.931，试用Matlab计算所需的精馏段理论板数。402022/9/24例5.5：氯仿-苯双组分精馏系统的气液412022/9/26解：因模型中的的函数关系是以表格形式给出，我们若要用辛普森法等计算较精确的精馏段理论板数的时候，需先采用拟合法将离散数据（ ， ）拟合成多项式，再将多项式代入被积函数求积。计算程序如下：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;s

30、p=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,sp)y=fnval(sp,x);f=1./(y-x-(0.9-y)./5);执行结果：N = 5所以，精馏段理论板数为5块。412022/9/24解：因模型中的的函数关系

31、是以表格形式给422022/9/26 5.2.2 高斯勒让德公式高斯法是一种精度较高的求积分法，它的一般公式是式中(x)是一个权重函数，Aj为系数，xj为横坐标上的节点 高斯勒让德多项式的权重函数(x)=1，因而，一个n点高斯勒让德求积公式具有如下形式：1 方法概述 422022/9/24 5.2432022/9/26 右边的f(xj)是函数f(x)在节点xj处的值， 节点xj是勒让德多项式Pn(x)的根。 Aj 为系数，其值为式中Pn(x)是勒让德多项式Pn(x)的一阶导数432022/9/24式中Pn(x)是勒让德多项式Pn(x442022/9/26pn(x)的前几项表达式为 由高斯勒让德