Ch5数值逼近

化工应用数学AppliedMathematicsforChemicalEngineering第5章数值逼近化学化工学院汤吉海2数值逼近的主要方法数值插值曲线拟合数值积分数值微分35.1数值插值在生产和科学实验中，自变量x与因变量y的函数y=f(x)的关系式有时不能直接写出表达式，而只能得到函数在若干个点的函数值或导数值（离散点）。当要求知道观测点之外的函数值时，需要估计函数值在该点的值。——插值内插：插值点在观测点范围之内外推：插值点在观测点范围之外数据点与插值点关系示意图4插值的步骤Step1：根据观测点的值，构造一个比较简单的函数y=φ(x)，使函数在观测点的值等于已知的数值或导数值。Step2：用简单函数y=φ(x)在点x处的值来估计未知函数y=f(x)在x点的值。关键问题：寻找函数φ(x)（办法是很多的）φ(x)可以是一个代数多项式，或是三角多项式，也可以是有理分式,如拉格朗日插值。φ(x)可以是任意光滑（任意阶导数连续）的函数或是分段函数，如分段低次插值、三次样条插值。函数类的不同，自然地有不同的逼近效果。5MATLAB的插值方法与命令一维数据插值：interp1二维数据内插值：interp2三维数据插值：interp3用快速Fourier算法作一维插值：interpft数据格点：griddata三次样条数据插值：splinen维数据插值：interpn生成用于画三维图形的矩阵数据：meshgrid生成用于多维函数计算或多维插值用的阵列：ndgrid6一维数据插值：interp1该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。yi=interp1(x,Y,xi)返回插值向量yi，每一元素对应于参量xi，同时由向量x与Y的内插值决定。参量x指定数据Y的点。若Y为一矩阵，则按Y的每列计算。yi是阶数为length(xi)*size(Y,2)的输出矩阵。yi=interp1(Y,xi)假定x=1:N，其中N为向量Y的长度，或者为矩阵Y的行数。yi=interp1(x,Y,xi,method)用指定的算法计算插值：‘nearest’：最近邻点插值，直接完成计算；‘linear’：线性插值（缺省方式），直接完成计算；‘spline’：三次样条函数插值。对于该方法，命令interp1调用函数spline、ppval、mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函数。命令spline用它们执行三次样条函数插值；‘pchip’：分段三次Hermite插值。对于该方法，命令interp1调用函数pchip，用于对向量x与y执行分段三次内插值。该方法保留单调性与数据的外形；‘cubic’：与‘pchip’操作相同；‘v5cubic’：在MATLAB5.0中的三次插值。7一维数据插值：interp1对于超出x范围的xi的分量，使用方法nearest、linear、v5cubic的插值算法，相应地将返回NaN。对其他的方法，interp1将对超出的分量执行外插值算法。yi=interp1(x,Y,xi,method,'extrap')对于超出x范围的xi中的分量将执行特殊的外插值法extrap。yi=interp1(x,Y,xi,method,extrapval)确定超出x范围的xi中的分量的外插值extrapval，其值通常取NaN或0。8一维数据插值：示例例1：一元函数插值例2：离散数据的一维插值x=0:10;y=x.*sin(x);xx=0:.25:10;yy=interp1(x,y,xx);plot(x,y,'kd',xx,yy)year=1900:10:2010;product=[75.99591.972105.711123.203131.669150.697179.323203.212226.505249.633256.344267.893];p1995=interp1(year,product,1995)x=1900:1:2010;y=interp1(year,product,x,'pchip');plot(year,product,'o',x,y)p1995=252.99012345678910-6-4-2024681900192019401960198020002020501001502002503009一维数据插值：示例例3：不同插值方法的比较做一个函数sin(x2)在区间[0，2]上40个函数值的表：用interp1来计算中间点的函数值，而不用sin(x)求：用样条插值来计算中间点的函数值：用sin计算得到的正确结果x=linspace(0,2*pi,40);y=sin(x.^2);values=00.604980.35591values=interp1(x,y,[0pi/23])better=interp1(x,y,[0pi/23],'spline')correct=sin([0pi/23].^2)better=00.624140.40982correct=00.624270.41212样条插值可得更高精度的结果，在表中用更多的数据点也可以使精度更好。10二维数据内插值：interp2ZI=interp2(X,Y,Z,XI,YI)返回矩阵ZI，其元素包含对应于参量XI与YI（可以是向量、或同型矩阵）的元素，即Zi(i,j)←[Xi(i,j),yi(i,j)]。用户可以输入行向量和列向量Xi与Yi，此时，输出向量Zi与矩阵meshgrid(xi,yi)是同型的。同时取决于由输入矩阵X、Y与Z确定的二维函数Z=f(X,Y)。参量X与Y必须是单调的，且相同的划分格式，就像由命令meshgrid生成的一样。若Xi与Yi中有在X与Y范围之外的点，则相应地返回nan（NotaNumber）。ZI=interp2(Z,XI,YI)缺省地，X=1:n、Y=1:m，其中[m,n]=size(Z)。再按第一种情形进行计算。ZI=interp2(Z,n)作n次递归计算，在Z的每两个元素之间插入它们的二维插值，这样，Z的阶数将不断增加。interp2(Z)等价于interp2(z,1)。ZI=interp2(X,Y,Z,XI,YI,method)用指定的算法method计算二维插值：’linear’：双线性插值算法（缺省算法）；’nearest’：最临近插值；’spline’：三次样条插值；’cubic’：双三次插值。11三次样条数据插值：splineyy=spline(x,y,xx)对于给定的离散的测量数据x,y（称为断点），要寻找一个三项多项式，以逼近每对数据(x,y)点间的曲线。该命令用三次样条插值计算出由向量x与y确定的一元函数y=f(x)在点xx处的值。若参量y是一矩阵，则以y的每一列和x配对，再分别计算由它们确定的函数在点xx处的值。则yy是一阶数为length(xx)*size(y,2)的矩阵。pp=spline(x,y)返回由向量x与y确定的分段样条多项式的系数矩阵pp，它可用于命令ppval、unmkpp的计算。12三次样条数据插值：示例例：对离散地分布在y=exp(x)sin(x)函数曲线上的数据点进行样条插值计算x=[024581212.817.219.920];y=exp(x).*sin(x);xx=0:.25:20;yy=spline(x,y,xx);plot(x,y,'o',xx,yy)02468101214161820-1012345x108135.2曲线拟合在科学技术及生产实践中，常常需要寻找某些参量之间的定量关系式，即有已知数据确定经验或半经验的数学模型，以便分析预测。当这些参量之间的数学关系式不能从理论上导出或者理论公式过于复杂时，常用的方法是将观测到的离散数据标记在平面图上，这只是对一个变量的情况而言，描成一条光滑的曲线（也包括直线，或对数坐标下的直线等）。为了进一步分析，在许多情况下，常常希望将曲线用一个简单的数学表达式加以描述，这就是曲线拟合，或是说经验建模。化工技术和生产中所采用的大量经验式和半经验式是通过最小二乘意义下的曲线拟合来定的。传递工程中的准数关联式；化工热力学中的比热容(或焓)与温度间的关系式；饱和蒸汽压方程、状态方程、活度系数方程；化工动力学中的动力学方程等半经验式；14曲线拟合经验的或半经验的数学关系式按其各待定系数与函数的关系可分为三类：线性关系，可用最小二乘法对试验数据直接进行拟合。非线性关系，但可将其线性化变为次生的线性关系，则可按照次生的线性关系式对实验数据进行最小二乘法拟合，然后再将最小二乘法用于原非线性关系式，最终求定待定参数。不能线性化的非线性关系，一般用优化方法提供初值，然后用最小二乘法求定参数最佳值。将最小二乘法用于线性关系式拟合，称作线性最小二乘法；直线函数拟合和代数多项式拟合是线性最小二乘法的两个特例。将最小二乘法用于非线性关系式拟合称作非线性最小二乘法。155.2.1多项式曲线拟合多项式曲线拟合有最小二乘法、牛顿法、拉格朗日法、牛顿一格雷高里法等多种方法，其中最小二乘法是最普遍的多项式拟合方法。在Matlab中，最小二乘拟合用函数polyfit()对一组数据进行定阶数的多项式拟合，其基本用法为：p＝polyfit(x,y,n)用最小二乘法对输入的数据x和y用n阶多项式进行逼近，函数返回多项式的系数，为一个长度为n+1的向量，包含多项式的系数。[p,s]＝polyfit(x,y,n)函数不仅返回多项式的系数向量p，还返回用函数polyval()获得的误差分析报告，保存在结构体变量S中。在多项式拟合中，如果n的值为1，就相当于用最小二乘法进行直线拟合。16多项式曲线拟合：示例例1某实验中测得一组数据，其值如下：xl2345y1.31.82.22.93.5已知x和y成线性关系，即y=kx+b，求系数k和bx=[12345];y=[1.31.82.22.93.5];[p,s]=polyfit(x,y,1)y1=polyval(p,x);plot(x,y1)holdon;plot(x,y,'b*');p=0.550.69s=R:[2x2double]df:3normr:0.16432可见，k为0.55，b为0.6911.522.533.544.5511.522.533.517多项式曲线拟合：示例例2：分别对向量x和y拟合3，4，5不同阶的多项式图形。x=[-3-1025.57];y=[3.34.52.01.52.5-1.2];p3=polyfit(x,y,3);%用向量x和y中元素拟合不同p4=polyfit(x,y,4);%次数的多项式p5=polyfit(x,y,5);xcurve=-3.5:0.1:7.2;%生成x值p3curve=polyval(p3,xcurve);%计算在这些x点的多项式值p4curve=polyval(p4,xcurve);p5curve=polyval(p5,xcurve);plot(xcurve,p3curve,'--',xcurve,p4curve,'-.',xcurve,p5curve,'-',x,y,'*');lx=[-11.5];ly=[00];holdon;plot(lx,ly,'--',lx,ly-1.3,'-.',lx,ly-2.6,'-');text(2,0,'degree3');text(2,-1.3,'egree4');text(2,-2.6,'degree5');holdoff;18多项式曲线拟合：示例可以看出，次数越高的多项式精度越好。5次的多项式经过所有的6个点。对于这个特殊的数据集合来说，这个多项式可称为内部插值多项式。-4-202468-3-2-101234567degree3egree4degree519多项式拟合的阶数选择问题例：计算函数在