第8章-matlab数值微分与积分课件

1、第8章 Matlab数值微分和积分化工学院 伍水生2019-11第8章 Matlab数值微分和积分化工学院 伍水生本章知识要点极限 (limit)数值微分(polyder, fnder)数值积分(quad, quadl, fnint)本章知识要点数值微分、数值积分在化工计算中的作用状态方程计算流体的焓和熵微分法反应动力学方程拟合等温活塞流反应器的设计计算微观离析反应器的计算数值微分、数值积分在化工计算中的作用状态方程计算流体的焓和单变量函数的极限格式1： L= limit( fun, x, x0)格式2： L= limit( fun, x, x0, left 或 right)8.1 数值微分8

2、.1.1 极限问题的解析解单变量函数的极限8.1 数值微分8.1.1 极限问题的解析例: 试求解极限问题 syms x a b; L=limit(x*(1+a/x)x*sin(b/x),x,inf) L = exp(a)*b例:求解单边极限问题 syms x; limit(exp(x3)-1)/(1-cos(sqrt(x-sin(x),x,0,right) ans =12例: 试求解极限问题多变量函数的极限：格式： L1=limit(limit(f,x,x0),y,y0) 或 L1=limit(limit(f,y,y0), x,x0) 如果x0 或y0不是确定的值，而是另一个变量的函数，如x-

3、g(y)，则上述的极限求取顺序不能交换。多变量函数的极限：例：求出二元函数极限值 syms x y a; f=exp(-1/(y2+x2) *sin(x)2/x2*(1+1/y2)(x+a2*y2); L=limit(limit(f,x,1/sqrt(y),y,inf)L =exp(a2)例：求出二元函数极限值8.2 函数导数的解析解函数的导数和高阶导数格式： y=diff(fun,x) %求导数(默认为1阶) y= diff(fun,x,n) %求n阶导数例： 求一阶导数： syms x; f1=diff(sin(x)/(x2+4*x+3),x)8.2 函数导数的解析解函数的导数和高阶导数s

4、yms x;f1=diff(sin(x)/x)得结果: f1=cos(x)/x-sin(x)/x2解: 输入命令: syms x; f1=diff(log(sin(x)得结果: f1=cos(x)/sin(x).在matlab中，函数lnx用log(x)表示，而log10(x)表示 lgx 。syms x;解: 输入命令: 在matlab中，函数l多元函数的偏导：格式： f=diff(diff(f,x,m),y,n) 或 f=diff(diff(f,y,n),x,m)例： 求偏导数 syms x y; z=(x2-2*x)*exp(-x2-y2-x*y); zx=diff(z,x)zx =(2

5、*x-2)*exp(-x2-y2-x*y)+(x2-2*x)*(-2*x-y)*exp(-x2-y2-x*y)多元函数的偏导： zy=diff(z,y)zy =(x2-2*x)*(-2*y-x)*exp(-x2-y2-x*y)直接绘制三维曲面 x,y=meshgrid(-3:.2:3,-2:.2:2); z=(x.2-2*x).*exp(-x.2-y.2-x.*y); surf(x,y,z), axis(-3 3 -2 2 -0.7 1.5) 第8章_matlab数值微分与积分课件隐函数的偏导数：格式：F=diff(f,xj)/diff(f,xi)隐函数的偏导数：例： z=f(x,y)=-x2

6、-y2-x*y; syms x y; f=x2-y2-x*y; diff(f,x)/diff(f,y)或者： syms x y; diff(x2-y2-x*y,x)/diff(x2-y2-x*y,y)ans =(-2*x-y)/(-2*y-x) 例： z=f(x,y)=-x2-y2-x*y;对于列表型函数往往需要用数值方法计算函数的微分数值微分的基本方法差分利用插值（拟合）多项式求微分利用三次样条插值（拟合）函数求微分数值微分可以放大误差，应谨慎使用对于列表型函数往往需要用数值方法计算函数的微分函数diff对于向量X，diff(X)表示了X(2)-X(1) X(3)-X(2) . X(n)-X

7、(n-1).对于矩阵X，diff(X)表示了X(2:n,:) - X(1:n-1,:) diff(x,n,dim)得到矩阵x在dim维上的n阶差值 diff(1:10).2)1 3 5 7 9 11 13 15 17 19x=1 3 8; 2 4 6diff(x,1,1)diff(x,1,2)diff(x,2,2)diff(x,3,2)1 1 -22 5;2 2Empty matrix: 2-by-030函数diff对于向量X，diff(X)表示了X(2)-X(Matlab数值微分实现方法有限差分法：用差分函数diff()近似计算导数多项式拟合方法：离散数据 多项式拟合函数 导函数pp 在xi

8、的导数值 三次样条插值方法样条拟合方法离散数据 样条插值函数cs导函数pp 在xi的导数值 离散数据 样条拟合函数sp导函数pp 在xi的导数值 Matlab数值微分实现方法有限差分法：用差分函数diff(例题某一液相反应浓度随时间变化的实验数据如下表：t/min00.20.61.02.05.010.0C(g/L)5.193.772.301.570.80.250.094试t=0.1,0.4min时的反应速度 x=0 0.2 0.6 1 2 5 10; C=5.19 3.77 2.30 1.57 0.8 0.25 0.094; f=spline(x,C); dC=fnder(f); dC1=fn

9、val(dC,0.1); dC2=fnval(dC,0.4); dC1, dC2三次样条插值方法:例题某一液相反应浓度随时间变化的实验数据如下表：t/min0 x=0 0.2 0.6 1 2 5 10; C=5.19 3.77 2.30 1.57 0.8 0.25 0.094; f=polyfit(x,C,3); dC=polyder(f); dC1=polyval(dC,0.1); dC2=polyval(dC,0.4);dC1, dC2多项式拟合方法： x=0 0.2 0.6 1 2 5 10;多项式拟8.2 数值积分数值积分在数值计算中有着重要作用,许多数值计算问题可以转化为数值积分问题

10、,如常微分方程初值问题等通常可用逼近多项式Pn(x)来代替被积函数f(x)，计算积分构造数值积分的方法很多，主要有Newton-Cotes系列数值积分法、Gauss积分法和Romberg积分法等8.2 数值积分数值积分在数值计算中有着重要作用,许多数值计数值积分MATLAB函数公式quad自适应Simpson求积公式（低阶）quadl自适应Lobatto求积公式；精度高，最常用trapz梯形求积公式；速度快，精度差cumtrapz梯形法求一个区间上的积分曲线cumsum等宽距法求一个区间上的积分曲线，精度很差fnint利用样条函数求不定积分；与spline，ppval配合使用，主要应用于表格“

11、函数”积分数值积分MATLAB函数公式quad自适应Simpson求积梯形法数值积分：trapz()调用格式：z=trapz(x,y) 用梯形求积方法计算y的积分近似值。对于向量y，trapz(y)返回y的积分；对于矩阵y，trapz(y)返回一行向量，向量中的各元素为矩阵y的对应列 向量的积分值；梯形法数值积分：trapz()调用格式：z=trapz(x格式： S=trapz(x,y)例： x1=0:pi/30:pi; y=sin(x1) cos(x1) sin(x1/2); S1=trapz(x1,y) S1 = 1.9982 0.0000 1.9995第8章_matlab数值微分与积分课

12、件等同于： x=0:pi/30:pi; y1=sin(x); y2=cos(x) ; y3=sin(x1/2); S1=trapz(x,y1) ; S2=trapz(x,y2) ; S3=trapz(x,y3) ;S1, S2, S3等同于：调用格式：q=quad(fun,a,b)q=quad(fun,a,b,tol)输入参数：fun被积函数。在定义fun时，被积函数表达式必须是向量形式，即表达式必须使用 点运算符(.*、./和.)以支持向量a,b即积分限a,btol 绝对误差限，默认值为1.e-6输出参数：q 积分结果自适应Lobatto法数值积分：quadl() 调用格式同quad自适应S

13、impson法数值积分：quad()调用格式：q=quad(fun,a,b)自适应Simp例：第一种，一般定义函数方法 y=quad(c3ffun,0,1.5)y = 0.9661第一种，一般定义函数方法 y=quad(c3ffun用 inline 函数定义被积函数： f=inline(2/sqrt(pi)*exp(-x.2),x); y=quad(f,0,1.5)y = 0.9661y=quad(2/sqrt(pi)*exp(-x.2), 0,1.5)y = 0.9661第8章_matlab数值微分与积分课件3.4.5 双重积分问题的数值解矩形区域上的二重积分的数值计算格式： 矩形区域的双重

14、积分： y=dblquad(Fun,xm,xM,ym,yM) 限定精度的双重积分： y=dblquad(Fun,xm,xM,ym,yM， )3.4.5 双重积分问题的数值解矩形区域上的二重积分的数值计例：求解 f=inline(exp(-x.2/2).*sin(x.2+y),x,y); y=dblquad(f,-2,2,-1,1)y = 1.57449318974494或者：function z=fun(x,y)z=exp(-x.2/2).*sin(x.2+y);y=dblquad(fun,-2,2,-1,1)或者： y=dblquad(exp(-x.2/2).*sin(x.2+y),-2,2

15、,-1,1)例：求解3.4.6 三重定积分的数值求解格式:I=triplequad(Fun,xm,xM,ym,yM,zm,zM,tol ,quadl) 其中quadl为具体求解一元积分的数值函数,也可选用quad或自编积分函数,但调用格式要与quadl一致。3.4.6 三重定积分的数值求解例： triplequad(inline(4*x.*z.*exp(-x.*x.*y-z.*z), x,y,z), 0, 1, 0, pi, 0, pi,1e-7,quadl)ans = 1.7328或者： triplequad(4*x.*z.*exp(-x.*x.*y-z.*z), 0, 1, 0, pi,

16、0, pi,1e-7,quadl)ans = 1.7328例：基于样条插值的数值积分运算格式： f=fnint(S)其中S为样条函数。例：考虑 中较稀疏的样本点,用样条积分的方式求出定积分及积分函数。 x=0,0.4,1 2,pi; y=sin(x); sp1=csapi(x,y); a=fnint(sp1,1); 建立三次样条函数 xx=fnval(a,0,pi); xx(2)-xx(1)ans = 2.0191基于样条插值的数值积分运算反应器停留时间分布的混合特性在t0的时刻，在一容器入口处突然向流进容器的流体脉冲注入一定量的示踪剂，同时在容器出口处测量流出物料中示踪剂浓度随时间的变化，实

17、验数据如下表：t/s020406080100120140160180200C(kmol/m3)10300000.45.516.211.11.70.10试计算流体在容器中的平均停留时间以及扩散准数。数学模型：平均停留时间方差： 扩散特征数为： 反应器停留时间分布的混合特性在t0的时刻，在一容器入口处突t = 0:20:200;C = 0, 0, 0, 0, 0.4, 5.5, 16.2, 11.1, 1.7, 0.1, 0; sp1= spline(t,C); %样条函数拟合 sp= spaps(t,C,1); %导函数 t0 = 0;tf = t(end); IC = quadl(Func,t

18、0,tf,sp); ICt = quadl(Func1,t0,tf,sp); ICt2 = quadl(Func2,t0,tf,sp); tm1 = ICt/IC; ss1 = ICt2/IC - tm12; DL2uL1 = ss1/tm12/2; tm1 , ss1 , DL2uL1function y = Func(x,sp) % f= Cy = fnval(sp,x); function y = Func1(x,sp) % f= Cty = fnval(sp,x).*x;function y = Func2(x,sp) % f= Ct2y = fnval(sp,x).*x.2;t = 0:20:200;C = 0, 0, 0,微分法进行动力学数据分析反应物A在一等温间歇反应器中发生反应为： ，测量得到反应器中不同时间下反应物A的浓度CA如下表所示。试根据表中数据确定其反应速率方程。数学模型：t/s0204060120180300CA(mol/L)10865321微分法进行动力学数据分析反应物A在一等温间歇反应器中发