计算方法插值法11223510 李晓东在许多实际问题及科学研究中，因素之间往往存在着函数关系，然而，这种关系经常很难有明显的解析表达， 通常只是一些离散数值。 有时即使给出了解析表达式， 却由于表达式过于复杂， 使用不便， 且不易于计算与分析。 解决这类问题我们往往使用插值法： 用一个“简单函数”)(x逼近被计算函数)(xf， 然后用)(x的函数值近似替代)(xf的函数值。插值法要求给出)(xf的一个函数表，然后选定一种简单的函数形式，比如多项式、分段线性函数及三角多项式等，通过已知的函数表来确定)(x作为)(xf的近似，概括地说，就是用简单函数为离散数组建立连续模型。一、 理论与算法（一）拉格朗日插值法在求满足插值条件n次插值多项式)(xPn之前，先考虑一个简单的插值问题：对节点),1,0(nixi中任一点)0(nkxk，作一 n 次多项式)(xlk，使它在该点上取值为1，而在其余点), 1, 1,1,0(nkkixi上取值为零，即kikixlik01)(（1.1 ）上式表明 n个点nkkxxxxx,1110都是 n次多项式)(xlk的零点，故可设)()()()(1110nkkkkxxxxxxxxxxAxl其中，kA为待定系数。由条件1)(kkxl立即可得)()()(1110nkkkkkkkxxxxxxxxA（1.2 ）故)()()()()()()(110110nkkkkkknkk kxxxxxxxxxxxxxxxxxl（1.3 ）由上式可以写出1n个 n 次插值多项式)(,),(),(10xlxlxln。我们称它们为在1n个节点nxxx,10上的 n次基本插值多项式或n次插值基函数。利用插值基函数立即可以写出满足插值条件的n次插值多项式)()()(1100xlyxlyxlynn（1.4 ）根据条件（ 1.1 ） ，容易验证上面多项式在节点ix处的值为),1,0(niyi，因此，它就是待求的 n次插值多项式)(xPn。形如式（ 1.4 ）的插值多项式就是拉格朗日插值多项式，记为)(xLn，即)()()()()()()()()()(11011001100nkkkkkknkknkknnnxxxxxxxxxxxxxxxxyxlyxlyxlyxL （1.5 ）为了便于上机，常将拉格朗日插值多项式（1.5 ）改写成：nknkjjjkj knxxxxyxL00)(（1.6 ）下图给出了利用式（ 1.6 ）计算x处函数值)(xLn的程序流程图：根据流程图用 matlab 语言编写的函数为：y输入x,y(i=1,2.,n)及yy+P*输出x,y Pj=1,2,.,n是 j=k? 否function y=lagrange(x0,y0,x) n=length(x0);m=length(x); for i=1:m s=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(x(i)-x0(j)/(x0(k)-x0(j); endends=p*y0(k)+s; endy(i)=s; end（二）分段线性插值法由于高次的插值多项式有不收敛现象, 有时会出现较大的误差, 不能有效的逼近被插函数, 所以人们提出了用分段的低次多项式逼近被插函数, 这就是分段插值方法。 常见的有分段线性插值和分段二次插值，本文使用的是分段线性插值法。当给定了 n 个点nxxx21上的函数值nyyy,21后，若要计算点ixx处函数值)(xf的近似值。可先选取两个节点ix与1ix，使1,iixxx，然后在小区间1,iixx上座线性插值，即得：iii i iii ixxxxyxxxxyxPxf11 11 1)(（2.1 ）相应的程序流程图如下：输入x ,y(i=1,2,.,n)输出按公式（2.1） 计算P1(x) x是否根据流程图用 matlab 语言编写的函数为：function y=seg_linear(x0,y0,x)n = length(x0);m=length(x);for i=1:mfor j=1:n-1if x(i)=x0(j)if x(i)=x0(j+1)y(i)=y0(j)*(x(i)-x0(j+1)/(x0(j)-x0(j+1)+y0(j+1)*(x(i)-x0(j)/(x0(j+1)-x0(j);endendendend（三）牛顿插值法设有函数xf，210,xxx, 为一系列互不相等的点， 称jixxxfxfxxfjiji ji,为xf关于点jixx ,的一阶差商。一般的称kkk kxxxxxfxxxfxxxf010110 10,，（3.1 ）为xf关于点kxxx,10的 k 阶差商。则xRxNxfnn其中nnnxxxfxxxxxxxxxfxxxxxxfxxxfxN,1011021010000（3.2 ）nnnxxxfxxxxxxxR,1110（3.3 ）显然，xNn是满足插值条件的至多n次的多项式。可得nixNxfini, 1 , 0。因而它是xf的n次的插值多项式。我们称xNn为牛顿插值多项式。在实际计算中，经常利用由式（3.1 ）构造出的差商表3.1 计算差商，即由表3.1中加下划横线的差商值直接构造插值多项式（3.2 ）kkxfx一阶差商二阶差商三阶差商四阶差商程序流程图为：yy+F(1,1输出F,y 输 入x,y(i=1,2,.,n)及j=2,3,.,ni=j,j+1,.,nF(i,j)F(i,j-1)-F(i-1,j-1)/(x-xt1t(x-x)*t j=2,3,.,yy+F(j,j)*t 表 3.1function y=newton(x0,y0,x)x0=x0(:);y0=y0(:);n=length(x0);m=length(x);F=zeros(n,n); %均差表二维矩阵，对角线上元素对应y=0*x; %插值点对应的近似值F(:,1)=y0; % 均差表第一列为已知节点函数值for j=2:nfor i=j:n F(i,j)=(F(i,j-1)-F(i-1,j-1)/(x0(i)-x0(i-j+1);endendfor j=1:m for k=2:n t=1; for i=1:k-1 t=(x(j)-x0(i)*t; endy(j)=y(j)+F(k,k)*t; endy(j)=y(j)+F(1,1); enddisp( 差商表： ) F=x0,F; disp(F)二、 例题分析给出)(xf的函数表x0.4 0.55 0.65 0.80 0.90 1.05 )(xf0.41075 0.57815 0.69675 0.88811 1.02652 1.25382 1、用 Lagrange 插值法计算)1.0()753(0.(0.596)42(0.ffff、的近似值；2、用分段线性插值法计算)1.0()753(0.(0.596)42(0.ffff、的近似值；3、构造差商表，用牛顿插值法计算)1.0()753(0.(0.596)42(0.ffff、的近似值。先将函数表存放到matlab 的 workspace 中以便调用：1、Lagrange 插值2、分段线性插值3、牛顿插值法拉格朗日插值的优点： 它的形式是对称的， 这样很容易编程上机实现。 它在理论上十分重要。牛顿插值的优点：在计算插值多项式及求解函数近似值都比较方便且计算量相对较小。从公式中可以看出：每增加一个节点，插值多项式只增加一项，因此便于递推运算，所以其具有灵活增加节点的优点。两者的不同点： 牛顿插值的计算量比拉格朗日要省，更利于程序设计。 由插值多项式的唯一性可知， n 次牛顿插值多项式与n 次拉格朗日插值多项式是等价的，它们只是表示形式不同。因此，牛顿余项和拉格朗日余项也是等价的，故matlab 计算结果二者一致。绘制所求点及已知节点在坐标中的位置，在matlab 中运行：plot(x,y3,rv,x,y2,g,x0,y0,-.y,x0,y0,ob) 得到坐标图：放大其中两组点可以看出， 牛顿插值法和拉格朗日插值法有较高的精度，误差很小，而分段线性插值的误差略大。但是牛顿仅对节点处的函数作了约束，如果插值条件再增加节点处对导数的限制的话， 解决这个问题的方法就是要利用 Hermite 插值多项式。但是由于高次插值存在不稳定性，会出现龙格（Runge ）现象，一般实际计算很少使用高次插值，更多的使用分段低次插值。