Newton插值法的程序设计与应用

1、引言在我们生活中，许多实际问题都都可以利用函数来表示某种内在规律的数量关 系，但是在很多应用领域，函数 f(x)有时不能直接写出表达式或者过于复杂不易计 算，而只能给出函数f(x)在若干个点上的函数值或导数值，通常也是造一张函数表， 当遇到要求表中未列出的变量的函数值时， 就必须做数值逼近.例如给定了函数f(x) 在几个特定点上的函数值，为了研究函数的变化规律，就要根据这个表，将其公式化 寻求某一函数y(x)去逼近f(x),在给定点上等于函数值，在其他点上约等于函数值， 这样y(x)既能反映f(x)的函数特性，又便于计算.称 y(x)的插值函数，f(x)为被插 信函数.我们可以求一个待定函数来

2、近似反映函数的特性，使得待定函数在给定点上等于 函数值，在其它点上的函数的值作为函数的近似值，这种方法称为插值法.利用差值函数很容易得到Lagrange插值多项式，公式结构紧凑，在理论分析中甚为方便.由于 公式中的每一项与所有的插值结点有关.因此，如果需要增加一个插值结点，则Lagrange插值公式中的每一项都需要改变计算量大，为了克服这一缺点，于是产生 了 Newton插值法.有的要求更高需要在某些点处的导数也相同于是前两种方法达不 到精度于是产生了 Hermite插值法.本题主要讨论牛顿插值多项式.主要是对牛顿插值法的理论总结，利用C+取序和Matlab语言程序编写，以及对其应用讨论.第一

3、章插值方法求得近似函数的方法通常有两种：一种是插值法，另一种是曲线拟合法.本章主 要讨论插值法用多项式来逼近列表函数问题，即具有唯一插值函数的多项式插值,对其中的多项式插值主要讨论n次多项式插值的方法，即给定 n+1各点处的函数值后， 怎样构造一个n次插值多项式的方法.§ 1.1 值问题的由来在许多实际问题及科学研究中，因素之间往往存在函数关系.然而这种关系经常 很难有明显的解析表达式，通常只是由观察与测试得到一些离散数值 .有时由于表达 式过于复杂，不仅使用不便，而且不易于进行计算与理论分析 .解决此类问题的方法 有两种：一种是插值法，另一种是拟合法.插值法是一种古老的数学方法，

4、他来自生产实践，早在一千多年前，我国科学家 在研究历法上就应用了线性插值与二次插值，但它的基本理论却是在微积分产生之后才逐渐完善的，具应用也逐步增多.特别是在计算机软件中，许三多库函数.等的计算 实际上归结于它的逼近函数的计算.在工程的实际问题中，由于航空、造船、精密机 械加工等实际问题的需要，插值方法在实践上和理论上显得更为重要，并得到了进一 步的发展和广泛的应用.因此我们希望能够得到一个“简单函数”逼近被计算函数的 函数值.于是就得到了一种方法.叫做插值逼近或者插值法.§ 1.2 个基本概念§ 1.2.1 值函数的概念及相关性质定义1.1 设连续函数y - f(x)在区

5、间a,b上有定义，已知在n+1个互异的点 Xo,X1，Xn上取值分别为y°, y1,yn (设a Mx。x 4 «b ).若在函数类中存在 以简单函数p(x),使得p(xi)= yi (i =0,1，n ),则称p(x)为f (x)的插值函数.称Xo,X1，Xn为插值节点,称a,b为插值区间.定义1.2设如下表XiXoX1 .XnV、VoV1 .yn求一个 p(x),使其满足：(1)p(x) w IPn(x) P(Xi) = f (Xi) = yi , i =0,1，n则称p(x)为f (x)的n次代数插值多项式，该问题称为n次代数插值问题.定理1.1n次代数插值问题的解存

6、在且唯§ 1.2.2 Lagrange插值函数的概念及相关性质定义1.3 若n次多项式1k(x) (k=0,1n)在n+1个互异节点x0,x1，xn上满 足%, j = k兀(七)=跖Jk,j=0,1;,n)P j # k则称1k(x) ( k =0,1，n )为节点x0,x1，xn上的n次插值基函数.n次差值多项式Ln(x)为nLn(x)=" yJk(x)k =0称之为拉格朗日多项式.定理1.2 设fwCn由a,b,则余项为:Rn(x)=fn1() (n 1)!Wn 1(x)其中 Wn i(X)=(x -Xo)(X -Xi) (x -xn).§ 1.2.3 商的

7、概念及相关性质定义1.4称函数f (x)于点Xi,Xj ( x #Xj)上的平均变化率为f(x)在Xi,Xj处的一阶差商，记作fx,Xj,即fX,Xjf(Xj)-f(x)Xj -XiXi 二 Xj .一阶差商的差商为f(x)在点xixXk处的二阶差商.记作flxixXk,即fXj,Xk 一 fXi,Xjfx,Xj,Xk:Xk -Xi般地，在求出f (x)的m -1阶差商后就可以构造f(x)的m阶差商.称f Xi0,Xi1, . ,XimfXi1,% - fXi0,xim . xi0为f (x)在为0?小xim处的m阶差商.特别地，规定零阶差商fxj = f (xi).fXo,Xi,X2二性质1

8、.1任意调换插值节点的次序其值不变即fXi,X0,X2 = fXo,X2,Xi性质1.2差商是一个数，可表示为某个函数的线性组合.性质1.3 当f(X)的m阶导数f(m)(X)在包含节点Xi0,Xi1，Xm的某个区间上存在时，在Xi0, 4,Xim之间至少存在一点之使一，f(m)()f Xi 0 , X 1, Xim =m!注1.1由性质2可以推出，n次多项式的一阶差商是n-1多项式.表1-1差商表Xkf(Xk)一阶差商二阶差商三阶差商四阶差商X0f(Xk)X1f(Xk)fx。, xX2f (Xk)fX1,X2f Xo,x,x2X3f(Xk)f X2,X3f X1, X2,X3fx,x,%,x

9、X4f (Xk)f X3,X4f X2, X3,X4f X1,X2, X3f%,x,x,x,x§ 1.2.3 Newton插值函数与差分 的概念及相关性质定义1.5我们定义牛顿插值公式为Nn(X)= f (X。)+ f Xo, x(x - X。)+ f Xo,X1,x2(x- Xo)(x- X1)(11)f X0,X1,Xn(X - X0)(X - X1) (X -Xni).定义1.6 设函数y = f (X)在等距节点Xk =x。+kh (k =0,1，n)上的函数值yk=f(x。(k=0,1，n),其中h为常数，称为步长.称函数在每个小区间上的增量y-yk为函数f(x)在Xk处以

10、h为步长的一阶向前差分，记作 ",即-7k =yk 1 -yk称一阶差分的差分Ayk书-Ayk为二阶差分，记作Hyk,即2yk =AykHtAyk = ykt2yk + yk由此可得m-1阶差分的差分myk=”'yk1 一即“为f(x)在Xk点处的m阶差分.称yk-yg为函数f(x)在Xk处以h为步长的一阶向后 差分，记作''yk.称Vmbk Vmyj为函数f(x)在Xk处以h为步长的m阶向后差分，记作Vmyk.即，r m, m 1I m' yk=yk-ykj表1-2向前差分表XiV、y、2y3yKvXoy1Xiy1y0X2y2K y。X3y3y22w

11、3y0X4VaK y2Nv。表1-3向后差分表X、y、y、V2ybybyXoy1X1y1VyoX2y2V2y0X3y3内2可矿y0X4Va内2五丁y0性质1.4各阶差分可以表示成函数值的线性组合.即n nyj = " (-1)kC：yj nA k =0 n(2) “j =" (7)kC：yj/k =0性质1.5差商与差分有如下关系fXo,Xi,Xnnfo n!hnfon!hn定义1.7 设n十1个等距插值节点的顺序为Xo,Xi，2二,4,则由差商表示的牛顿插值公式为Nn(X) = fXo fX),X1(X-Xo) fX),X1,X2(X-Xo)(X-X1) fXo,X1,

12、,Xn(X-Xo) (X- XnA)又设插值节点x = x0th (0 ： t ： n)由性质1.5差商与差分关系fX0,Xi； ,Xk:可得N n(X)= fo.2f. . t(t -1) (t-n 1) .nf!_ f o '- Ton!n=zk ZZ0t(t 1) (t k 1)(1.2)k!称为牛顿前插公式.其余项公式为fnZ )n1 _Rn(x)= Rn(Xo th)= -(-t(t -1) (t-n)h ,Xo,Xn(n 1)!定义1.8 设n+1个等距插值节点的顺序为4,4,X1,X0,则由差商表示的牛顿 插值公式为Nn(X)=fXn fXn,Xn(XXn)fXn,XnL

13、Xn/(XXn)(XXn)(XXn/),f Xn , Xn 1, Xo( X -Xn ) (X -X1 )又设插值节点x = Xn -th (o : t :二 n)由性质1.5差商与差分关系fXn,Xn4, Xn-k=fnk!hk可得、 t(t -1)、2n t(t -1) (t -n 1) nNn(x) = fn Yfn ' J fn -( -1)fn(1.3)2!n!= 1(7)5(tk1)：kfn小k!称为牛顿后插公式.其余项公式为fn1( )n1 n1 .R(X) =Rn(Xn-th)t(t-1) (t-n)(-1) h ,Xo,Xn(n 1)!第二章Newton插值法的程序设

14、计与应用利用插值基函数很容易得到 Lagrange插值多项式，公式结构紧凑，在理论分析 中也很方便，但当插值节点增减时全部插值基函数L(x) i =0,1，n都要随之变化，整个公式也将发生变化，这在实际计算中是很不方便的，于是我们导出插值多项式的 又一种表示形式一一Newton插值公式.远远减少了工作量.但由于通过相同结点的差 值多项式是惟一的，因此，拉格朗日差值公式与牛顿公式本质上是一个多项式，只不过形式不同而已.为了更方便的应用Newton插法，本章主要讨论Newton插值法的程序设计，并对其 应用进行讨论.§ 2.1 Newton插值法的算法与程序设计本节主要讨论n次多项式插值

15、的方法，即给定n+1各点处的函数值后，怎样构造 一个n次插值多项式的方法.为了更方便的计算并与计算机结合，我们便作出了相应 的程序.并利用Matlab进行绘图. 2.1.1 2.1.1 Newton插值法的应用步骤步骤 首先我们按照表1-1,求得各点的差商.然后利用Newton前差或后差公式 即(公式1.3)把数值带入.即可以求得n次多项式.它在计算机上的应用步骤如下：步骤1输入所要求的牛顿多项式的次数n ,并依次输入n+1个节点(Xi,yi). for : i=0 to n+1scanf("%f',&xi) scanf("%f",&y0i

16、);步骤2计算各阶差商 for : i=1 to n+1for : j=i to n+1 if(i>1)yij=(y1j-y11)/(xj-xi);elseyij=(yi-1j-yi-1j-1)/(xj-xj-1);步骤3代入牛顿插值公式，可计算得出结果.for : i=1 to n+1temp=temp*(xx-xi-1);newton=newton+yii*temp;由于篇幅所限，流程图、C语言程序和Matlab程序详见附录(附录A、B、C). 2.1.2 利用Newton插值法程序的实例为了更方便的应用牛顿插值法，我们进行了与计算机的结合，下面我们将展示几 个例子.例2.1已知y

17、= sin x的一组数据为0冗兀冗JIx6432sinx01V2V31222(1)构造牛顿插值函数并作图分析.(2)并分别利用程序估计sin-, sin啊的估计值.1215分析 首先我们可以通过程序求出差商表：xkf (xk)一阶差商二阶差商三阶差商四阶差商006120.954929268714三 20.791089631-2.0860762L3而 20.607024424-0.35153865-0.13648909匹 210.25587263-0.44710035-0.091254700.028797106带入定义1.5中可求得如下牛顿插值多项式如下冗N2(x) =0.95492926x -

18、0.2086076x(x -)(2.1)JIJIN3(x) =0.95492926x -0.2086076x(x - ) -0.13648909x(x -一)(x -一) 664(2.2)71JE31N4(x) = 0.95492926x -0.2086076x(x ) - 0.13648909x(x-)(x -)(2.3)664JIJJ0.028797106x(x -)(x -)(x -)643第二步 我们分别利用C+?序和Matlab分别计算sin工和sin粤的值.步骤如下： 1215利用C+E序：首先我们输入所要求的牛顿多项式的次数n,然后/&入n+1个节点的值.即可以得出sin

19、工和sin巴的值为0.2586和0.0804 ; 1215利用Matlab程序：首先我们在Matlab中新建一个文件，文件内容即为附录C中Matlab程序代码，然后输入以下步骤：x=0,pi/6,pi/4,pi/3,pi/2;y=0,0.5,sqrt(2)./2,sqrt(3)./2,1;x0=pi/12;Newton(x,y,x0ans =0.2586x=0,pi/6,pi/4,pi/3,pi/2;y=0,0.5,sqrt(2)./2,sqrt(3)./2,1;x0=14.*pi/15;Newton(x,y,x0) ans =0.0804 即结果为0.2586和0.0804 利用Matlab

20、分别对求得的牛顿差值多项式作图(代码详见附录D)图形如下.-1%IIM2X ifcpi/12.O.250e)闾i-如g)XI1-D511.52图2-1例2.2设f (x) =ln(1 +x)的函数表如下:x0.250.300.360.390.450.2231440.2623640.3074850.3293040.371564试计算 ln(1.275) , f(x) =ln2.3分析：同上题步骤我们先求差商表，进而代入公式可得N4(x) =0.22314 0.7844(x-0.25) 0.294393939(x - 0.25)(x - 0.30)+ 0.1411736357(x-0.25)(x-

21、0.30)(x-0.36)+0.057720032( 2.4 )(x -0.25)(x -0.30)(x-0.36)(x -0.39)N3(x) = 0.22314 + 0.7844(x 0.25) -0.294393939(x-0.25)(x-0.30)250.1411736357(x -0.25)(x -0.30)(x -0.36)M(x) =0.22314 +0.7844(x0.25) -0.294393939 (x - 0.25)(x - 0.30)(2.6)利用C+?序我们可以得到计算结果：W(0.275) =0.242946,电(0.275) = 0.242945,N2(0.275

22、) =0.242938 N4(1.3) = 0.928823 , N3(1.3) = 0.877002 , N2(0.275) = 0.737653我们利用Matlab进行绘图.1.J上上1一0 .cP :.遍NFSJ/i明出处100,511,5图2-2从上述两个例子我们可以看出，多项式(2.1 )在区间0,1周围与原函数逼近的 较好.离这个区间越远与原函数的误差越大在 x =1.5处时，该图像就已经开始背离图 像了.所以该多项式只能在一个小白区间里可以逼近原函数，不适合作为原函数的逼近函数.也可以看出多项式(2.2)在0,；区间的周围逼近的较好，但是x = 2处时， 该图像就离原图像误差较大

23、.多项式(2.3)在区间0,2.5都逼近的挺好，从图中我 们看出在远离这个区间的图像误差相对较大，但是在这三个多项式中是逼近最好的，同样在例2.2中也是如此.于是我们我们可以得出节点越多，函数逼近的相对较好.在节点附近逼近的越好， 越远离节点误差越大，所以公式适用于计算节点附近的值于是为了减小误差，在下一节的等距节点下的插值公式根据所求的点的函数值的不同分别采取了前插和后插公 式.§ 2.2 距节点下的Newton插值算法与程序设计前面我们讲述了一般节点下的牛顿插值公式，为了计算方便于是有了对等距节点下的牛顿多项式的研究，本节将对等距节点下的插值多项式进行总结讨论.§2.2

24、.1 等距节点下的Newton插值法的程序算法步骤步骤 我们按照表1-2或1-3 ,求差分.然后利用Newton前插公式(公式1.2)或 Newton后插公式(公式1.3)并把数值带入.即可以求得n次多项式.它在计算机上的应用步骤如下：步骤1输入所要求的牛顿多项式的次数n，步长h，并依次输入n+1个节点(Xi, yi).for : i=0 to n+1scanf("%f',&xi) scanf("%f',&y0i);步骤2求得各界差分for : i=1 to n+1 for : j=i to n+1yij=(yi-1j+1-yi-1j-1)/

25、(xj-xj);/求向前差分for : i=1 to n+1 for : j=i to n+1yij=(yi-1j-yi-1j-1)/(xj-xj-i)；/求向后差分步骤3代入牛顿插值公式，可计算得出结果for(i=1;i<n+1;i+) temp=temp*(t-i+1)/i);newton=newton+yii*temp;printf(" 求得的结果为：N(%.4f)=%9fn",xx,newton);/求得运用前插公式的值for(i=n;i<0;i-)temp=temp*(t-i+1)/i);newton=newton+yn-in-i*temp;print

26、f("求得的结果为：N(%.4f)=%9fn",xx,newton);/求得运用后插公式的值由于篇幅所限C语言程序(详见附录D E).§222等距节点下Newton插值的实例例2.4 已知tan x的值列表如下:x1.31.311.321.33tanx3.60213.74713.90334.0723近似计算 tan 1.325, tan1.305.采用Newton向后插公式.为此，做差分表xifiV2fiV3fi1.33.60210.14500.01120.00161.313.74710.15620.01281.323.90330.16901.334.0723从而

27、，有p3(x) =4.0723 0.1690t 00128t(t 1) 0.0016 t(t 1)(t 2).2!3!1.325 -1.33将t = -0.5代入上式,得0.01tan1.325 p3(1.325) =3.9869.将t = 1.305 -1.33 =_2.5带入后插多项式中可以得到0.01tan1.305 5(1.305) =3.0797现在我们利用C+E序首先输入所求插值的次数3和步长0.01.然后输入各个节点，并输入所要求的点1.325既可以求出该点的函数值.即f (1.325) =3.9869. f (1.26) ： p3(1.26) =3.0797卜面我们可得图图2-

28、4例2.4设函数在各节点的取值如下X00.20.40.60.81.0f凶1.00.8187310.6703200.5488120.4493290.367879试利用插值公式求f (0.3)的值.采用Newton向前插公式，同上题我们先做差分表，然后相应带入到差分公式 . 中求得后插公式.利用C语言程序步骤如下：首先输入所求插值的次数 5和步长0.2.然后输入各个节点，并输入所要求的点 0.3既可以求出该点的函数值.即f(0.3) =0.7880.由以上例子我们看到例1用了牛顿后插公式，例2用了牛顿前插公式，我们该怎 样选取.这个经过验证得出，如果所要求的点较靠近节点 %,则采用前插公式；如果

29、靠近Xn ,则采用牛顿后插公式.附录A:牛顿插值法的流程图Ni<nYj<n?i>1输入i<n+1输入n/分别输入 n+1/个输入各节点YNNYY/输出n的值yi-ij -yi-i,j-iy Xj - Xjyi-1j - yi-1,j-1 yj 二Xj - jt=t*(x- Xi J )N=N+y iii=1;t=1;N= y°o开始i=i+1i=i+1j=j+1j=ii=0附录B:牛顿插值法的程序#include<stdio.h>void main()float x11,y1111,xx,temp,newton;int i,j,n；printf(&

30、quot;Newton插值:n请输入要运算的值:x=");scanf("%f",&xx);printf("请输入插值的次数(n<11):n=");scanf("%d",&n);printf("请输入 %dffl 值:n",n+1);for(i=0;i<n+1;i+) printf("x%d=",i);scanf("%f",&xi);printf("y%d=",i);scanf("%f",&a

31、mp;y0i);for(i=1;i<n+1;i+)for(j=i;j<n+1;j+) if(i>1)yij=(yi-1j-yi-1j-1)/(xj-xj-i);elseyij=(yi-1j-yi-1j-1)/(xj-xj-1);printf("%fn",yii);temp=1;newton=y00;for(i=1;i<n+1;i+) temp=temp*(xx-xi-1);newton=newton+yii*temp;printf("求得的结果为：N(%.4f)=%9fn",xx,newton);附录C:牛顿插值法的Matlab程

32、序function f = Newton(x,y,x0)syms t;if(length(x) = length(y)n = length(x);c(1:n) = 0.0;elsedisp('x 和y的维数不相等！');return;endf = y(1);y1 = 0;l = 1;for(i=1:n-1)for(j=i+1:n)y1(j) = (y(j)-y(i)/(x(j)-x(i);endc(i) = y1(i+1);l = l*(t-x(i);f = f + c(i)*l;simplify(f);y = yi;if(i=n-1)if(nargin = 3)f = sub

33、s(f,'t',x0);elsef = collect(f);%将插值多项式展开f = vpa(f, 6);endendend附录D:等距节点下白前插公式的C语言程序#include<stdio.h> void main() float x11,y1111,xx,temp,newton,t,h;int i,j,n；printf("Newton插值:n请输入要运算的值:x=");scanf("%f",&xx);printf("请输入插值的次数(n<11):n=");scanf("%d",&n);printf("步长:n请输入要运算的值：h=");scanf("%f",&h);printf(" 请输入 %dffl 值:n&q