第四部分 插值、拟合与数值微分和积分天津科大海洋学院天津科大海洋学院2009-10 本章知识要点•插值方法(interp,spline)•拟合方法(polyfit,csaps)•数值微分(polyder, fnder)•数值积分(quad, quadl, fnint) 插值、拟合、数值微分、数值积分在化工计算中的作用•表格式物性数据的内插•离散实验数据点的处理•状态方程计算流体的焓和熵•微分法反应动力学方程拟合•等温活塞流反应器的设计计算•微观离析反应器的计算 插值简介插值的数学问题可以描述为：已知n+1个数对{xi, f(xi)}，其中i＝0,1…n，（xi互不相同，称之为节点），求取函数g(xi)=f(xi)当{xi, f(xi)}有相当的精确度，但它们的函数关系难以确定或难以计算时，则可利用这些数据点来构造一个较简单的函数来近似表达原函数关系根据逼近函数的不同，常见的插值方法：根据逼近函数的不同，常见的插值方法： •Lagrange多项式插值（线性插值）多项式插值（线性插值）•分段插值分段插值•三次样条插值三次样条插值•三角插值三角插值•有理式插值有理式插值 插值方法的选择已知熔盐在423～473K的密度和粘度如下表所示，估计450K时的密度和粘度。
温度(K)密度(kg/m3)粘度(Pa·S)4231976177.584331967146.514431959122.794531951104.6463194390.26473193478.79 Matlab的插值（Interpolation）函数插插值值方法方法Matlab函数函数插插值值方法方法Matlab函数函数一维插值interp1使用FFT方法的一维插值interpft快速一维插值interpq分段三次Hermite插值pchip二维插值interp2三次样条插值spline三维插值interp3N维插值interpn 一维插值一维插值interp1调用格式：yi＝interp1(T,P,Ti， spline) 已知数据向量（x,y），计算并返回在插值向量xi处的函数值yi=interp1(x,y,xi, ‘method’)yi=interp1(x,y,xi, ‘method’, ‘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.5440];(y＝sinx)，比较一维线性、线性最近、立方和三次样条插值所得xi＝0,0.15,0.30,0.45,…,10处的值yi如果初始数据点为x＝0,2,4,…,10，y＝sinx，以上方法插值效果 一维插值一维插值方法比较方法比较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,'线性插值\rightarrow')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])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与pchipSpline()的调用格式为：yi=spline(x,y,xi) 此函数等同于yi=interp1(x,y,xi, ‘spline’)pp=spline(x,y) 返回三次样条插值的分段多项式形式的向量spline函数可以保证插值函数的三阶导数连续pchip()的调用格式为：yi=pchip(x,y,xi) 此函数等同于yi=interp1(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)])grid on, axis equallegend('Orignal data','pchip','spline') 其它样条插值函数－csapecsape()的调用格式为：pp=csape(x,y)pp=csape(x,y,conds)conds可为以下字符串：‘complete’指定端点处一阶导数‘second’指定端点处二阶导数‘periodic’左右端点处一、二阶导数相等‘not-a-knot’第二个和倒数第二个节点处三阶导数连续‘variational’自然边界条件，即端点二阶导数为零除csape外，Matlab的样条工具箱还提供了其它样条插值函数，如csapi，spapi等 csape的使用设f(x)为区间[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’) 几个有用的函数1.[breaks coff k m n]=unmkpp(pp)命令可以看到样条函数的具体信息其中coff是一个矩阵，其第i行是第i个三次多项式的系数，多项式形式如下：2.fnval或ppval可以计算样条函数在指定点的函数值3.fnder和fnint可分别计算样条函数的导数和积分4.fnplt可用于绘制样条函数的图形 二维插值：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,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) 拟合简介•拟合和插值的区别在于：1.拟合时，所得函数不需要过所有数值点2.插值函数不宜外推，拟合函数在某些情况则可以•拟合方法中最常用的是最小二乘曲线拟合•最小二乘法的基本思路是使拟合因变量y在给定点xi上使残差平方和最小•用于拟合的函数可以是多样的，本讲只介绍多项式函数拟合和样条函数拟合 最小二乘多项式拟合：最小二乘多项式拟合：polyfit调用格式如下：p=polyfit(x,y,n)[p,s]=polyfit(x,y,n)说明：n输入参数：(x,y)为已知数据向量，n为多项式阶数；输出参数p为拟合生成的多项式的系数向量（长度为n+1），s为结构参数，供函数polyval()调用以获得误差估计值n函数polyval()常常与polyfit()联合使用，其调用格式为：y=polyval(p,x)n [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.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,'k-'),hold on,text(1.931,9.532,'七次多项式拟合\leftarrow'),hold offx0.51.01.52.02.53.0y1.752.453.814.808.008.60已知下表数据，用polyfit进行多项式拟合 最小二乘法拟合生成样条曲线函数名曲线类型拟合准则是否平滑处理csaps()三次样条曲线最小二乘法是spap2()B样条曲线最小二乘法否spaps()B样条曲线最小二乘法是•拟合得到曲线函数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时，相当于“自然的”三次样条拟合 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 拟合得到的样条函数 功能：平滑生成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')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;endx0.51.01.52.02.53.0y1.752.453.814.808.008.60对于微小的数据扰动，多项式拟合通常比样条函数更为敏感 数值微分•对于列表型函数往往需要用数值方法计算函数的微分•数值微分的基本方法1.差分2.利用插值（拟合）多项式求微分3.利用三次样条插值（拟合）函数求微分•数值微分可以放大误差，应谨慎使用 Matlab数值微分实现方法有限差分法：用差分函数diff()近似计算导数多项式拟合方法：离散数据 多项式拟合函数 导函数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=[1 3 8; 2 4 6]diff(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') 函数gradientFX = gradient(F)[FX,FY] = gradient(F)[Fx,Fy,Fz] = gradient(F)[...] = gradient(F,h)[...] = gradient(F,h1,h2,...)其中1.Fx=dF/dx，FY=dF/dy，Fz=dF/dz；当F为向量时，dF=gradient(F)是一维梯度2.h1, h2，……是步长，缺省值为13.如果h1，h2，……为向量,其长度必须匹配F的维数 例题某一液相反应浓度随时间变化的实验数据如下表：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];plot(x,C,'bo'),hold onpp=csaps(x,C);fnplt(pp)dC=fnder(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)respectively\n',dC1,dC2) 数值积分•数值积分在数值计算中有着重要作用,许多数值计算问题可以转化为数值积分问题,如常微分方程初值问题等•通常可用逼近多项式Pn(x)来代替被积函数f(x)，计算积分•构造数值积分的方法很多，主要有Newton-Cotes系列数值积分法、Gauss积分法和Romberg积分法等 数值积分数值积分MATLAB函数函数公式公式quad自适应Simpson求积公式（低阶）quadl自适应Lobatto求积公式；精度高，最常用trapz梯形求积公式；速度快，精度差cumtrapz梯形法求一个区间上的积分曲线cumsum等宽距法求一个区间上的积分曲线，精度很差fnint利用样条函数求不定积分；与spline，ppval配合使用，主要应用于表格“函数”积分 梯形法数值积分：梯形法数值积分：trapz()调用格式：z=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,b]tol 绝对误差限，默认值为1.e-6p1,p2…… 直接传递给函数fun的已知参数输出参数：q 积分结果自适应Lobatto法数值积分：quadl() 调用格式同quad自适应Simpson法数值积分：quad() 不同积分函数的比较不同积分函数的比较求积分：比较cumsum，trapz，quad，quadl的积分精度，该积分的精确解为0.9008407878function 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);err(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. 无穷积分可先选取一个有限的积分区间，如[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只能处理积分限为常数的多重积分，对内积分上下限为只能处理积分限为常数的多重积分，对内积分上下限为外积分变量的积分问题外积分变量的积分问题,需自行编写函数求解。
需自行编写函数求解 样条函数在数值积分与微分中的应用已知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(ppval(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 fnplt(pp,'b-');hold onfnplt(int_pp,'m:'),fnplt(der_pp,'r--');hold offlegend('y(x)','S(x)','dy/dx') ppval(int_pp,[1,2])给出一个行向量，其中第一个元素是原函数在0～1的定积分，第2个元素是原函数在0～2的定积分，所以在1～2的定积分应是第2个元素减第1个元素的值 反应器停留时间分布的混合特性在t＝0的时刻，在一容器入口处突然向流进容器的流体脉冲注入一定量的示踪剂，同时在容器出口处测量流出物料中示踪剂浓度随时间的变化，实验数据如下表：t/s020406080100120140160180200C(kmol/m3)×10300000.45.516.211.11.70.10试计算流体在容器中的平均停留时间以及扩散准数。
数学模型：平均停留时间方差： 扩散特征数为： 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/m^3)¡Á10^3')hold off% By spaps():figure，plot(t,C,'o')sp2= spaps(t,C,1); hold onfnplt(sp2,'b-'); xlabel(‘Time (s)’)，ylabel('C (kmol/m^3)¡Á10^3')hold offsp=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 - tm1^2DL2uL1 = ss1/tm1^2/2% ------------------------------------------------------------------function 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= Ct^2y = fnval(sp,x).*x.^2; 微分法进行动力学数据分析反应物A在一等温间歇反应器中发生反应为： ，测量得到反应器中不同时间下反应物A的浓度CA如下表所示。
试根据表中数据确定其反应速率方程数学模型：t/s0204060120180300CA(mol/L)10865321 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(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) 。