1、第8章 Matlab数值微分和积分化工学院化工学院 伍水生伍水生2019-11本章知识要点 极限(limit)数值微分(polyder,fnder)数值积分(quad,quadl,fnint)数值微分、数值积分在化工计算中的作用状态方程计算流体的焓和熵微分法反应动力学方程拟合等温活塞流反应器的设计计算微观离析反应器的计算 单变量函数的极限 格式1：L=limit(fun,x,x0)格式2：L=limit(fun,x,x0,left 或 right)8.1 数值微分8.1.1 极限问题的解析解 例:试求解极限问题 syms x a b;L=limit(x*(1+a/x)x*sin(b/x),x,
2、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-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
3、(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)的导数。求例xxysin1syms x;syms x;f1=diff(sin(x)/x)f1=diff(sin(x)/x)得结果得结果:f1=cos(x)/x-sin(x)/x2的导数。求例xysinln2解解:输入命令输入命令:syms x;f1=diff(log(sin(x)得结果得结果:f1=cos(x)/sin(x).在
4、在matlab中，函数中，函数lnx用用log(x)表示，而表示，而log10(x)表示表示 lgx。多元函数的偏导：格式：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*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
5、: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)隐函数的偏导数：格式：F=diff(f,xj)/diff(f,xi)例：z=f(x,y)=-x2-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)对于列表型函数往往需要用数值方法计算函数的微分 数值微分的基本方法1.差分2.利用插值（拟合）多项式求微分3.利用三
6、次样条插值（拟合）函数求微分 数值微分可以放大误差，应谨慎使用函数diff对于向量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,:)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-030Matlab数值微分实现方法有限差分法：用差分函
7、数diff()近似计算导数多项式拟合方法：离散数据()polyfit多项式拟合函数()polyder导函数pp()polyval在xi的导数值 三次样条插值方法样条拟合方法离散数据 样条插值函数cs导函数pp ()fnval在xi的导数值()spline()fnder离散数据 样条拟合函数sp导函数pp 在xi的导数值 ()fnder()()2()spapsaspapcsaps或或()fnval例题某一液相反应浓度随时间变化的实验数据如下表：t/min00.20.61.02.05.010.0C(g/L)5.193.772.301.570.80.250.094试t=0.1,0.4min时的反应速
8、度 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=fnval(dC,0.1);dC2=fnval(dC,0.4);dC1,dC2三次样条插值方法: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多项式拟合方法：8.2 数值积分数值积分在数值计算中
9、有着重要作用,许多数值计算问题可以转化为数值积分问题,如常微分方程初值问题等通常可用逼近多项式Pn(x)来代替被积函数f(x)，计算积分构造数值积分的方法很多，主要有Newton-Cotes系列数值积分法、Gauss积分法和Romberg积分法等bandxxP)(数值积分数值积分MATLAB函数函数公式公式quad自适应自适应Simpson求积公式（低阶）求积公式（低阶）quadl自适应自适应Lobatto求积公式；精度高，最求积公式；精度高，最常用常用trapz梯形求积公式；速度快，精度差梯形求积公式；速度快，精度差cumtrapz梯形法求一个区间上的积分曲线cumsum等宽距法求一个区间上
10、的积分曲线，精度很差fnint利用样条函数求不定积分；与spline，ppval配合使用，主要应用于表格“函数”积分梯形法数值积分：梯形法数值积分：trapz()调用格式：z=trapz(x,y)用梯形求积方法计算y的积分近似值。对于向量y，trapz(y)返回y的积分；对于矩阵y，trapz(y)返回一行向量，向量中的各元素为矩阵y的对应列 向量的积分值；格式：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 等同于：x=0:pi/30:pi;y1=sin(x
11、);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被积函数。在定义在定义funfun时，被积函数表达时，被积函数表达式必须是向量形式，即表达式必须使用式必须是向量形式，即表达式必须使用 点运算符点运算符(.(.*、././和和.).)以支持向量以支持向量a,b即积分限a,btol 绝对误差限，默认值为1.e-6输出参数：q 积分结果自适应Lobatto法数值积分：quadl()调用格式同quad自
12、适应Simpson法数值积分：quad()例：第一种，一般定义函数方法 y=quad(c3ffun,0,1.5)y=0.9661 用 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.96613.4.5 双重积分问题的数值解 矩形区域上的二重积分的数值计算 格式：矩形区域的双重积分：y=dblquad(Fun,xm,xM,ym,yM)限定精度的双重积分：y=dblquad(Fun,xm,xM,ym,yM，)例：求解 f=inl
13、ine(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,-1,1)3.4.6 三重定积分的数值求解 格式:I=triplequad(Fun,xm,xM,ym,yM,zm,zM,tol ,quadl)其中quadl为具体求解一元积分的数值函数,也可选用quad或自编积分函数,但调用格
14、式要与quadl一致。例：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,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);建
15、立三次样条函数 xx=fnval(a,0,pi);xx(2)-xx(1)ans=2.01910sin xdx反应器停留时间分布的混合特性在t0的时刻，在一容器入口处突然向流进容器的流体脉冲注入一定量的示踪剂，同时在容器出口处测量流出物料中示踪剂浓度随时间的变化，实验数据如下表：t/s020406080100120140160180200C(kmol/m3)10300000.45.516.211.11.70.10试计算流体在容器中的平均停留时间以及扩散准数。数学模型：00CdtCtdttm平均停留时间方差：20022mtCdtdtCt扩散特征数为：222mLtuLDt=0:20:200;C=0,
16、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,t0,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=f
17、nval(sp,x).*x;function y=Func2(x,sp)%f=Ct2y=fnval(sp,x).*x.2;微分法进行动力学数据分析反应物A在一等温间歇反应器中发生反应为：，测量得到反应器中不同时间下反应物A的浓度CA如下表所示。试根据表中数据确定其反应速率方程。数学模型：产物At/s0204060120180300CA(mol/L)10865321nAAkCr kCnrAAlnln)ln(%动力学数据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);pp=fnder(sp);%计算B样条函数的导函数 dCAdt=fnval(pp,t);%计算t处的导函数值 rA=dCAdt;y=log(-rA);x=log(CA);p=polyfit(x,y,1);k=exp(p(2)n=p(1)
