数学建模插值方法

.,插值与拟合,.,前言,函数是多种多样的，在科研与工程实际中有的函数表达式过于复杂而不便于计算，但又需要计算多点的函数值；有的函数甚至给不出数学式子，只能通过实验和测量得到一些离散数据（如某些点的函数值和导数值）。面对这种情况，很自然的一个想法就是构造某个简单的函数作为要考察的函数的近似。如果要求近似函数满足给定的离散数据，则称之为插值函数。实用上，我们常取结构相对比较简单的代数多项式作为插值函数,这就是所谓的代数插值。,.,插值部分,.,定理1设为给定的彼此互异的个插值节点，则存在唯一的次数不超过的多项式，满足条件,.,二、存在性与唯一性,.,证明:设,其中为待定系数.利用插值条件,我们得到一个线性代数方程组,其中观察发现矩阵A是范德蒙矩阵,那么,由几代知识知道矩阵A的行列式为,由定理中条件,插值结点为彼此互异的,那么行列式不为零.故由Cramer法则知线性代数方程组存在唯一解.,.,三、Lagrange插值法,（1）Lagrange插值多项式可以表示为,.,引入记号,易证,从而Lagrange插值多项式可表示为,.,（2）插值误差估计定理2设在上连续，在内存在,节点，是拉格朗日插值多项式，则对任意,插值余项其中且依赖于.,.,例2.求过点(2,0)(4,3)(6,5)(8,4)(10,1)的拉格朗日型插值多项式。解：用4次插值多项式对5个点插值,.,.,于是有,.,functionyi=lagrcz(x,y,xi)n=length(x);m=length(xi);fors=1:myi(s)=0;fori=1:nw(i)=1;dw(i)=1;forj=1:nif(j=i)w(i)=(xi(s)-x(j)*w(i);dw(i)=(x(i)-x(j)*dw(i);endendyi(s)=y(i)*w(i)/dw(i)+yi(s);endend,.,.,缺点:当增加或减少插值节点时,基函数需要重新构造,不便于实际的计算使用,.,.,(2)Newton插值公式由差商定义把以上各式由后向前代入,可得,.,差商表,.,例2:已知求满足以上插值条件的牛顿型插值多项式。解:,.,由上述差商表对角线上取得的值则牛顿三次插值多项式为,.,functionyi=newtcz(x,y,xi)n=length(x);m=length(xi);nt=zeros(n,n);nt(:,1)=y;fori=2:nforj=i:nnt(j,i)=(nt(j,i-1)-nt(j-1,i-1)/(x(j)-x(j-(i-1);endEndfori=1:nnt(i,i)Endfori=1:myi(i)=nt(1,1);forj=2:nt=1;fors=1:j-1t=t*(xi(i)-x(s);endyi(i)=yi(i)+t*nt(j,j);endend,.,五、Hermite插值多项式,给定的是节点上的函数值和导数值,问题：已知,求3次多项式，使得,.,.,*多项式插值的问题,前面介绍了构造插值公式的方法，并分析了它们的余项。在实际应用插值函数作近似计算时，总希望插值公式余项的绝对值小一些，即使得逼近的精度好。从表达式看，似乎提高插值多项式的次数便可达到目的，但实际上并非如此。,.,例如给定函数,取其等距节点，构造的Lagrange插值多项式为当时，只能在内收敛，而在这个区间以外是发散的。这种畸形现象通常叫做Runge现象。如下图所示。,.,.,六、分段插值,所谓分段插值，就是将被插值函数逐段多项式化。在每个子段上构造插值多项式，然后把它们装配在一，作为整个区间上的插值函数，即称为分段多项式。如果函数在每个子段上都是次式，则称为次式。,一般（低次：k=1,2,3）,.,（1）分段线性插值的构造（k=1)易知在每个子区间上是一次插值多项式分段线性插值的余项其中,.,（2）分段抛物线插值（K=2）,（3）分段三次Hermite插值(K=3),.,（4）三次样条插值,在分段插值中，分段线性插值在节点上仅连续而不可导，分段三次埃尔米特插值有连续的一阶导数，如此光滑程度常不能满足物理问题的需要，而引入的样条函数则可以同时解决这两个问题,使插值函数既是低阶分段函数,又是光滑的函数。,.,则称是在该区间上关于该划分的一个三次样条函数。,.,其中四个待定系数为,子区间共有n个所以要确定S(x)需要4n个待定系数。另一方面,要求分段三次多项式S(x)及其导数和在整个插值区间a,b上连续,则要求它们在各个子区间的连接点上连续，即满足条件,由样条函数的定义可知,三次样条插值函数S(x)是一个分段三次多项式,要求出S(x),在每个小区间xi,xi+1上要确定4个待定参数,若用Si(x)表示它在第i个子区间xi,xi+1上的表达式，则,.,（1）插值条件（2）连接条件式共给出了4n-2个条件,而待定系数有4n个,因此还需要2个条件才能确定S(x),通常在区间端点上各加一个条件,称为边界条件,常用边界条件有三种类型。,.,第一种类型：给定两端点的一阶导数值：第二种类型：给定两端点f(x)的二阶导数值：作为特例,称为自然边界条件。满足自然边界条件的三次样条插值函数称为自然样条插值函数。第三种类型：当是以为周期的函数时，则要求S(x)也是周期函数,这时边界条件应满足当时，,.,这样，由上给定的任一种边界条件加上插值条件和连接条件，就能得出4n个方程，可以惟一确定4n个系数。从而得到三次样条插值函数S(x)在各个子区间xi,xi+1上的表达式S(xi）(i=1,2,)。但是，这种做法当n较大时，计算工作很大，不便于实际应用。因此我们希望找到一种简单的构造方法。,.,三次样条插值函数的求法设S(x)在节点xi处的二阶导数为因为在子区间xi-1,xi上是三次多项式,所以在此小区间上是x的线性函数,且因为用线性插值,可知其表达式为,记，则有,.,其中,Ai,Bi为积分常数,可利用插值条件确定,即要求Ai,Bi满足并记，则得,连续两次积分得,.,由上讨论可知,只要确定这n+1个值,就可定出三样条插值函数S(x)。为了求出,利用一阶导数在子区间连接点上连续的条件，求导一次,得在区间xi-1,xi上的表达式为,.,也就是在右端点xi上有,在左端点xi-1上有,将上式中的i-1改为i,即得在子区间xi,xi+1上的表达式,并由此得,利用在内接点的连续性,即就可得到关于参数的一个方程,.,上式两边同乘以,即得方程,若记,.,则所得方程可简写成,即,这是一个含有n+1个未知数、n-1个方程的线性方程组.要完全确定的值还需要补充两个条件,这两个条件通常根据实际问题的需要，根据插值区间a,b的两个端点处的边界条件来补充。边界条件的种类很多，常见的有以下3种：,.,第一种边界条件：即已知插值区间两端的一阶导数值：则可得到包含Mi的两个线性方程,S(x)在子区间上的导数为,由条件得,即,同理,由条件得,.,即得确定的线性方程组,其中,.,第二种边界条件:即已知插值区间两端的二阶导数值:,由于在区间端点处二阶导数，所以方程中实际上只包含有n-1个未知数，从而得方程组,.,第三种边界条件:由与，可得和,其中,.,得关于的线性方程组。,利用线性代数知识,可以证明方程组的系数矩阵都是非奇异的，因此有惟一解。,.,用三次样条绘制的曲线不仅有很好的光滑度，而且当节点逐渐加密时，其函数值在整体上能很好地逼近被插函数，相应的导数值也收敛于被插函数的导数，不会发生龙格现象。因此三次样条在计算机辅助设计中有广泛的应用。,.,用MATLAB作插值计算,一维插值函数：,yi=interp1(x，y，xi，method),nearest：最邻近插值linear：线性插值；spline：三次样条插值；cubic：立方插值。缺省时：分段线性插值。,注意：所有的插值方法都要求x是单调的，并且xi不能够超过x的范围。,.,例：在1-12的11小时内，每隔1小时测量一次温度，测得的温度依次为：5，8，9，15，25，29，31，30，22，25，27，24。试估计每隔1/10小时的温度值。,x=1:12;y=589152529313022252724;xi=1:0.1:12;yi=interp1(x,y,xi,spline);plot(x,y,+,xi,yi,r),.,.,三次样条插值的Matlab实现如果三次样条插值没有边界条件，最常用的方法，就是采用非扭结（not-a-knot）条件。这个条件强迫第1个和第2个三次多项式的三阶导数相等。对最后一个和倒数第2个三次多项式也做同样地处理。Matlab中三次样条插值也有现成的函数：y=interp1(x0,y0,x,spline)；y=spline(x0,y0,x)；pp=csape(x0,y0,conds,valconds)，y=ppval(pp,x)。其中x0,y0是已知数据点，x是插值点，y是插值点的函数值。,.,对于三次样条插值，我们提倡使用函数csape，csape的返回值是pp形式，要求插值点的函数值，必须调用函数ppval。pp=csape(x0,y0)：使用默认的边界条件，即Lagrange边界条件。pp=csape(x0,y0,conds,valconds)中的conds指定插值的边界条件，其值可为：complete边界为一阶导数，一阶导数的值在valconds参数中给出，若忽略valconds参数，则按缺省情况处理。not-a-knot非扭结条件periodic周期条件second边界为二阶导数，二阶导数的值在valconds参数中给出，若忽略valconds参数，二阶导数的缺省值为0,0。variational设置边界的二阶导数值为0,0。,.,对于一些特殊的边界条件，可以通过conds的一个12矩阵来表示，conds元素的取值为0，1，2。conds(i)=j的含义是给定端点i的j阶导数，即conds的第一个元素表示左边界的条件，第二个元素表示右边界的条件conds=2,1表示左边界是二阶导数，右边界是一阶导数，对应的值由valconds给出。,.,最小二乘法拟合,已知一批离散数据(xi,yi),i=0,1,.,n，且x0x1xn,寻找一个函数f(x)，使达到最小.这个过程称为最小二乘拟合,f(x)称为拟合函数.,拟合部分,.,一、线性拟合,若设拟合函数f(x)=b+ax，则有令,.,即,这是一个关于a,b的2元线性方程组.求解即可得到f(x)的表达式.,.,二、多项式拟合有时所给数据点的分布并不一定近似地呈一条直线,这时仍用直线拟合显然是不合适的,可用多项式拟合。对于给定的一组数据寻求次数不超过m(mN)的多项式，,来拟合所给定的数据，与线性拟合类似，使偏差的平方和,为最小,.,由于可以看作是关于(j=0,1,2,m)的多元函数,故上述拟合多项式的构造问题可归结为多元函数的极值问题。令,得,即有,.,这是关于系数的线性方程组，通常称为正规方程组。可以证明，正规方程组有惟一解。,.,三、可化为线性拟合的非线性拟合有些非线性拟合曲线可以通过适当的变量替换转化为线性曲线，从而用线性拟合进行处理，对于一个实际的曲线拟合问题，一般先按观测值在直角坐标平面上描出散点图，看一看散点的分布同哪类曲线图形接近，然后选用相接近的曲线拟合方程。再通过适当的变量替换转化为线性拟合问题，按线性拟合解出后再还原为原变量所表示的曲线拟合方程。下表列举了几类经适当变换后化为线性拟合求解的曲线拟合方程及变换关系,.,曲线拟合方程变换关系变换后线性拟合方程,.,多项式曲线拟合函数：polyfit()调用格式：p=polyfit(x,y,n)p,s=polyfit(x,y,n)说明：x,y为数据点，n为多项式阶数，返回p为幂次从高到低的多项式系数向量p。矩阵s用于生成预测值的误差估计。,例2：由离散数据,拟合出多项式。,.,x=0:.1:1;y=.3.511.41.61.9.6.4.81.52n=3;p=polyfit(x,y,n)xi=linspace(0,1,100);z=polyval(p,xi);plot(x,y,o,xi,z,k:,x,y,b),.,.,二维插值,第一种（网格节点）：,.,注意：最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。,最邻近插值,二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值即为所求。,.,将四个插值点（矩形的四个顶点）处的函数值依次简记为：,分片线性插值,f(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,yj+1)=f4,.,双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：,其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。,双线性插值,.,(3)分片双三次样条插值,.,第二种（散乱节点）：,.,要求x0,y0单调；x，y可取为矩阵，或x取行向量，y取为列向量，x,y的值分别不能超出x0,y0的范围。,z=interp2(x0,y0,z0,x,y,method),用MATLAB作网格节点数据的插值,nearest最邻近插值linear双线性插值cubic双三次插值缺省时,双线性插值,.,例：测得平板表面3*5网格点处的温度分别为：828180828479636165818484828586试作出平板表面的温度分布曲面z=f(x,y)的图形。,输入以下命令：x=1:5;y=1:3;temps=8281808284;7963616581;8484828586;mesh(x,y,temps),1.先在三维坐