1、第四部分 插值、拟合与数值微分和积分,本章知识要点,插值方法(interp,spline)拟合方法(polyfit,csaps)数值微分(polyder, fnder)数值积分(quad, quadl, fnint),插值、拟合、数值微分、数值积分在化工计算中的作用,表格式物性数据的内插离散实验数据点的处理状态方程计算流体的焓和熵微分法反应动力学方程拟合等温活塞流反应器的设计计算微观离析反应器的计算,插值简介,插值的数学问题可以描述为：已知n+1个数对xi, f(xi)，其中i0,1n，（xi互不相同，称之为节点），求取函数g(xi)=f(xi)。当xi, f(xi)有相当的精确度，但它们的函
2、数关系难以确定或难以计算时，则可利用这些数据点来构造一个较简单的函数来近似表达原函数关系。根据逼近函数的不同，常见的插值方法： Lagrange多项式插值（线性插值）分段插值三次样条插值三角插值有理式插值,插值方法的选择,已知熔盐在423473K的密度和粘度如下表所示，估计450K时的密度和粘度。,Matlab的插值（Interpolation）函数,一维插值interp1,调用格式：yiinterp1(T,P,Ti， spline) 已知数据向量（x,y），计算并返回在插值向量xi处的函数值yi=interp1(x,y,xi, method)yi=interp1(x,y,xi, method
3、, extrap)method用于指定插值算法，其值可以是：nearest最近插值linear线性插值（默认值）spline分段三次样条插值pchip分段三次Hermite插值cubic与pchip相同,注意：向量x为单调。若y为矩阵，则对y的每一列进行插值向量xi中有元素不在x的范围内，则对应yi值为NaN extrap用于指定当向量xi中有元素不在x的范围内时，采用method所指定的插值算法进行外插计算与之对应的yi值,已知x=0:10，y=0 0.8415 0.9093 0.1411 -0.7568 -0.9589 -0.2794 0.6570 0.9894 0.4121 -0.544
4、0;(ysinx)，比较一维线性、线性最近、立方和三次样条插值所得xi0,0.15,0.30,0.45,10处的值yi。如果初始数据点为x0,2,4,10，ysinx，以上方法插值效果。,一维插值方法比较,x=0:1:10;y =0 0.8415 0.9093 0.1411 -0.7568 -0.9589 -0.2794 0.657 0.9894 0.4121 -0.5440;plot(x,y,bo),hold onfplot(sin,0 10) xi=0:0.15:10;yi=interp1(x,y,xi);plot(xi,yi,r+),text(0.7028,0.4649,线性插值righ
5、tarrow)yi2=interp1(x,y,xi,nearst);plot(xi,yi2,g*),text(3.537,0.1374,leftarrow最近插值)yi3=interp1(x,y,xi,cubic);plot(xi,yi3,mx),text(2.408,0.8333,leftarrow三次插值)yi4=interp1(x,y,xi,spline);plot(xi,yi4,k:),text(4.62,0.8158,三次样条插值rightarrow),初始数据对于插值的影响,x=0:2:10;y=sin(x);plot(x,y,go),hold onezplot(sin,0 10)
6、xi=0:0.15:10;yi=sin(xi);xj=0:10;yi1=interp1(x,y,xj,spline);plot(xj,yi1,r+),yj2=interp1(xi,yi,xj,spline);plot(xj,yj2,k*),spline与pchip,Spline()的调用格式为：yi=spline(x,y,xi) 此函数等同于yi=interp1(x,y,xi, spline)pp=spline(x,y) 返回三次样条插值的分段多项式形式的向量spline函数可以保证插值函数的三阶导数连续,pchip()的调用格式为：yi=pchip(x,y,xi) 此函数等同于yi=inte
7、rp1(x,y,xi, pchip)pp=pchip(x,y) 返回三次样条插值的分段多项式形式的向量,pchip与spline的区别,利用点(x=sin(k/6),y=cos(k/6)，其中k=0 1 2 3来逼近单位圆的前四分之一圆周。比较pchip与spline的差别。,pchip与spline的区别：三次样条在相邻的节点上并不保证单调性；而Hermite分段三次样条则可保证插值的局部单调性,t=linspace(0,pi/2,4)x=cos(t);y=sin(t);xx=linspace(0,1,40);plot(x,y,s,xx,pchip(x,y,xx);spline(x,y,xx
8、)grid on, axis equallegend(Orignal data,pchip,spline),其它样条插值函数csape,csape()的调用格式为：pp=csape(x,y)pp=csape(x,y,conds),conds可为以下字符串：complete指定端点处一阶导数second指定端点处二阶导数periodic左右端点处一、二阶导数相等not-a-knot第二个和倒数第二个节点处三阶导数连续variational自然边界条件，即端点二阶导数为零,除csape外，Matlab的样条工具箱还提供了其它样条插值函数，如csapi，spapi等,csape的使用,设f(x)为区
9、间0,3上的函数，剖分节点为xi=0 1 2 3；节点上的函数值为yi=0 0.5 2 1.5，左右端点处的一阶导数值分别为0.2和-1。试求区间0,3上满足上述条件的三次样条插值函数,程序：x=0 1 2 3;y=0.2 0 0.5 2.0 1.5 -1;pp=csape(x,y,complete),几个有用的函数,breaks coff k m n=unmkpp(pp)命令可以看到样条函数的具体信息其中coff是一个矩阵，其第i行是第i个三次多项式的系数，多项式形式如下：fnval或ppval可以计算样条函数在指定点的函数值fnder和fnint可分别计算样条函数的导数和积分fnplt可用
10、于绘制样条函数的图形,二维插值：interp2,调用格式：zi=interp2(x,y,z,xi,yi,method)method算法属性值可以是；nearest最近插值linear线性插值（默认）spline三次样条插值(spline)cubic立方插值,二维插值函数的使用,假设有一组分度系数的“海底深度测量数据”，由以下一段程序生成：randn(state,2);x=-5:5;y=-5:5;X,Y=meshgrid(x,y);Z=-500+1.2*exp(-(X-1).2+(Y-2).2)-0.7*exp(-(exp(X+2).2+(Y+1).2);surf(X,Y,Z)，view(-25
11、,25)试由插值方式绘制海底形状图。,xi=linspace(-5,5,50);yi=linspace(-5,5,50);XI,YI=meshgrid(xi,yi);ZI=interp2(X,Y,Z,XI,YI,*cubic);surf(XI,YI,ZI)view(-25,25),拟合简介,拟合和插值的区别在于：拟合时，所得函数不需要过所有数值点插值函数不宜外推，拟合函数在某些情况则可以拟合方法中最常用的是最小二乘曲线拟合最小二乘法的基本思路是使拟合因变量y在给定点xi上使残差平方和最小用于拟合的函数可以是多样的，本讲只介绍多项式函数拟合和样条函数拟合,最小二乘多项式拟合：polyfit,调用
12、格式如下：p=polyfit(x,y,n)p,s=polyfit(x,y,n),说明：输入参数：(x,y)为已知数据向量，n为多项式阶数；输出参数p为拟合生成的多项式的系数向量（长度为n+1），s为结构参数，供函数polyval()调用以获得误差估计值。函数polyval()常常与polyfit()联合使用，其调用格式为：y=polyval(p,x)。 y,deltay=polyval(p,xi,s), 返回xi处的拟合函数值，deltay是50%置信区间的y的误差。,多项式次数对拟合效果的影响,x=0.5:0.5:3;y=1.75 2.45 3.81 4.80 8.00 8.60;xi=0.
13、5:0.05:3;a2,S2=polyfit(x,y,2);y2,delta2=polyval(a2,xi,S2);plot(x,y,ro,xi,y2,m-),text(1.04,4.67,二次多项式拟合leftarrow)hold on,a4,S4=polyfit(x,y,4);y4,delta4=polyval(a4,xi,S4);plot(x,y,ro,xi,y4,b-),hold on,text(1.563,6.775,四次多项式拟合leftarrow)a7,S7=polyfit(x,y,7);y7,delta7=polyval(a7,xi,S7);plot(x,y,ro,xi,y7,
14、k-),hold on,text(1.931,9.532,七次多项式拟合leftarrow),hold off,已知下表数据，用polyfit进行多项式拟合,最小二乘法拟合生成样条曲线,拟合得到曲线函数sp以后，可利用fnval()计算任意自变量下的函数值。,功能：平滑生成三次样条函数，即对于数据(xi,yi)，所求的三次样条函数y=f(x)满足 ：,函数csaps()的用法,调用格式：sp=csaps(x,y,p) ys=csaps(x,y,p,xx,w)输入参数：x,y 要处理的离散数据(xi,yi) p 平滑参数，取值区间为0,1。当p=0时，相当于最小二乘直线拟 合；当p=1时，相当于
15、“自然的”三次样条拟合 xx 用于指定在给定点xx上计算其三次样条函数值(ys) w 权值（权重），默认为1输出参数：sp 拟合得到的样条函数 ys 在给定点xx上的三次样条函数值,功能：用最小二乘法拟合生成B样条曲线，即对于数据(xi,yi)，所求k次样条函数y=f(x)满足 ：,函数spap2 ()的用法,调用格式：sp=spap2(knots,k,x,y) sp=spap2(knots,k,x,y,w)输入参数：knots节点序数(knot sequence)k 样条函数的阶次，一般取k=3，有时取k=4x,y要处理的离散数据(xi,yi)w权值（权重），默认为1输出参数：sp 拟合得到
16、的样条函数,功能：平滑生成B样条函数,函数spaps ()的用法,调用格式：sp=spaps(x,y,tol) sp,ys=spaps(x,y,tol,m,w)输入参数：x,y要处理的离散数据(xi,yi) tol光滑时的允许精度 m 默认值是2，即平滑生成三次B样条曲线w权值（权重），默认为1输出参数：sp拟合得到的样条函数 ys在x上经平滑处理的B样条函数值,采用样条拟合函数csaps进行拟合，比较p的取值对于拟合效果的影响,函数csaps()的用法,x=0.5:0.5:3;y=1.75 2.45 3.81 4.80 8.00 8.60;xi=0.5:0.05:3;plot(x,y,ro)
17、hold on,C=0 0 0;1 0 0;0 0 1; 0 1 1; 0.5 0.5 0.5;j=1;for i=0:0.25:1.0 Cj=C(j,:); ys(j,:)=csaps(x,y,i,xi); plot(xi,ys(j,:),color,Cj),pause j=j+1;end,对于微小的数据扰动，多项式拟合通常比样条函数更为敏感,数值微分,对于列表型函数往往需要用数值方法计算函数的微分数值微分的基本方法差分利用插值（拟合）多项式求微分利用三次样条插值（拟合）函数求微分数值微分可以放大误差，应谨慎使用,Matlab数值微分实现方法,有限差分法：用差分函数diff()近似计算导数,
18、多项式拟合方法：,离散数据,多项式拟合函数,导函数pp,在xi的导数值,三次样条插值方法,样条拟合方法,离散数据,样条插值函数cs,导函数pp,在xi的导数值,离散数据,样条拟合函数sp,导函数pp,在xi的导数值,函数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 19,利用diff函数求sin(x)在0,10上的导数值,x
19、=1 3 8; 2 4 6diff(x,1,1)diff(x,1,2)diff(x,2,2)diff(x,3,2),1 1 -2,2 5;2 2,Empty matrix: 2-by-0,30,fplot(cos,1 10)hold onh=1;x=1:h:10;hh=0.01;xx=1:hh:10;diffy=diff(sin(x)/h;diffyy=diff(sin(xx)/hh;plot(1:h:10-h,diffy,r:,1:hh:10-hh,diffyy,k-o),函数gradient,FX = gradient(F)FX,FY = gradient(F)Fx,Fy,Fz = gra
20、dient(F). = gradient(F,h). = gradient(F,h1,h2,.),其中Fx=dF/dx，FY=dF/dy，Fz=dF/dz；当F为向量时，dF=gradient(F)是一维梯度h1, h2，是步长，缺省值为1如果h1，h2，为向量,其长度必须匹配F的维数,例题,某一液相反应浓度随时间变化的实验数据如下表：,试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;plot(x,C,bo),hold onpp=csaps(x,C);fnplt(pp)dC=fnder(
21、pp);dC1=ppval(dC,0.1)dC2=fnval(dC,0.4)fprintf(The reaction rate at 0.1min and 0.4min are%5.4f,%5.4f.(g/(L min)respectivelyn,dC1,dC2),数值积分,数值积分在数值计算中有着重要作用,许多数值计算问题可以转化为数值积分问题,如常微分方程初值问题等通常可用逼近多项式Pn(x)来代替被积函数f(x)，计算积分构造数值积分的方法很多，主要有Newton-Cotes系列数值积分法、Gauss积分法和Romberg积分法等,数值积分,梯形法数值积分：trapz(),调用格式：z=
22、trapz(y) 用梯形求积方法计算y的积分近似值。对于向量y，trapz(y)返回y的积分；对于矩阵y，trapz(y)返回一行向量，向量中的各元素为矩阵y的对应列 向量的积分值；,调用格式：q=quad(fun,a,b)q=quad(fun,a,b,tol)q=quad(fun,a,b,tol,trace,p1,p2,)输入参数：fun被积函数。在定义fun时，被积函数表达式必须是向量形式，即表达式必须使用点运算符(.*、./和.)以支持向量a,b即积分限a,btol 绝对误差限，默认值为1.e-6p1,p2 直接传递给函数fun的已知参数输出参数：q 积分结果自适应Lobatto法数值积
23、分：quadl() 调用格式同quad,自适应Simpson法数值积分：quad(),不同积分函数的比较,求积分：,比较cumsum，trapz，quad，quadl的积分精度，该积分的精确解为0.9008407878,function Cha4demo5format longd=pi/1000;t=0:d:3*pi;nt=length(t);y=fun(t);sc=cumsum(y)*d;scf=sc(nt)z=trapz(y)*dqd=quad(fun,0,3*pi,1)qd2=quadl(fun,0,3*pi,1)EV=0.9008407878;err(1)=abs(scf-EV);er
24、r(2)=abs(z-EV)*10;err(3)=abs(qd-EV)*10;err(4)=abs(qd2-EV)*10;bar(err),title(不同积分方法比较),colormap(summer)text(0.7,7e-4,矩形法,fontsize,20)text(1.5,5e-4,梯形法,fontsize,20)text(2.4,3e-4,Simpson法,fontsize,20)text(3.4,1e-4,Lobatto法,fontsize,20)%-function y=fun(t)y=exp(-0.5*t).*sin(t+pi/6),广义积分,1. 奇点积分,2. 无穷积分,可
25、先选取一个有限的积分区间，如0,100计算；在选择一个较大的积分区间，如0 200计算，如两次计算结果的差满足一定的精度要求，则可认为此值即为无穷积分的值,fun=inline(1./(sqrt(x).*(exp(x)+1);quadl(fun,0,1)quadl(fun,eps,1),多重数值积分,二重积分函数：SS=dblquad(fun,xmin,xmax,ymin,ymax,tol,method)三重积分函数：SSS=triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol,method),说明:Matlab只能处理积分限为常数的多重积分，对内积
26、分上下限为外积分变量的积分问题,需自行编写函数求解。,样条函数在数值积分与微分中的应用,已知x=(0:0.1:1)*2*pi;y=sin(x) ，求其积分和导数,x=(0:0.1:1)*2*pi;y=sin(x); pp=spline(x,y); int_pp=fnint(pp); der_pp=fnder(pp); % 在基础区间上，计算三个样条函数与理论值的最大误差xx=(0:0.01:1)*2*pi;err_yy=max(abs(ppval(pp,xx)-sin(xx)err_int=max(abs(ppval(int_pp,xx)-(1-cos(xx)err_der=max(abs(p
27、pval(der_pp,xx)-cos(xx) % 计算y(x)在区间1,2上的定积分DefiniteIntegral.bySpline=ppval(int_pp,1,2)*-1;1; DefiniteIntegral.byTheory=(1-cos(2)-(1-cos(1);% 计算dy(3)/dxDerivative.bySpline=fnval(der_pp,3);Derivative.byTheory=cos(3);Derivative.byDiference=(sin(3.01)-sin(3)/0.01; %前向差分近似DefiniteIntegral,Derivative fnpl
28、t(pp,b-);hold onfnplt(int_pp,m:),fnplt(der_pp,r-);hold offlegend(y(x),S(x),dy/dx),ppval(int_pp,1,2)给出一个行向量，其中第一个元素是原函数在01的定积分，第2个元素是原函数在02的定积分，所以在12的定积分应是第2个元素减第1个元素的值,反应器停留时间分布的混合特性,在t0的时刻，在一容器入口处突然向流进容器的流体脉冲注入一定量的示踪剂，同时在容器出口处测量流出物料中示踪剂浓度随时间的变化，实验数据如下表：,试计算流体在容器中的平均停留时间以及扩散准数。,数学模型：,平均停留时间,方差：,扩散特征
29、数为：,function Cha4demo7t = 0:20:200;C = 0, 0, 0, 0, 0.4, 5.5, 16.2, 11.1, 1.7, 0.1, 0;plot(t,C,o)sp1= spline(t,C);hold onfnplt(sp1,b-); xlabel(Time (s)，ylabel(C (kmol/m3)103)hold off% By spaps():figure，plot(t,C,o)sp2= spaps(t,C,1); hold onfnplt(sp2,b-); xlabel(Time (s)，ylabel(C (kmol/m3)103)hold offs
30、p=input(Which function should be used to integrate?);% Integrationt0 = 0;tf = t(end);IC = quadl(Func,t0,tf,sp); ICt = quadl(Func1,t0,tf,sp); ICt2 = quadl(Func2,t0,tf,sp);tm1 = ICt/ICss1 = ICt2/IC - tm12DL2uL1 = ss1/tm12/2% -function y = Func(x,sp) % f= Cy = fnval(sp,x); % -function y = Func1(x,sp) %
31、 f= Cty = fnval(sp,x).*x;% -function y = Func2(x,sp) % f= Ct2y = fnval(sp,x).*x.2;,微分法进行动力学数据分析,反应物A在一等温间歇反应器中发生反应为： ，测量得到反应器中不同时间下反应物A的浓度CA如下表所示。试根据表中数据确定其反应速率方程。,数学模型：,function Cha4demo8 % 动力学数据t = 0 20 40 60 120 180 300;CA = 10 8 6 5 3 2 1;% 用最小二乘样条拟合法计算微分dCA/dtknots = 3;K = 3; % 三次B样条sp = spap2(
32、knots,K,t,CA);sp = spap2(newknt(sp),K,t,CA);pp = fnder(sp); % 计算B样条函数的导函数dCAdt = 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)figurefnplt(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),
