牛顿插值法的C语言学习知识实现

1、.*牛顿插值法的C 语言实现摘要：拉格朗日插值法具有明显的对称性，公式中的每一项与所有的插值节点有关。因此， 如果需要增加一个插值节点，则拉格朗日插值公式中的每一项都要改变， 在有的应用中就显得不太方便。 因此，可以利用另外一种差值方法来弥补这种缺陷，就牛顿插值法。 本文通过对牛顿插值法的数学分析，主要给出其C 语言实现方法。关键字：差商差分C 语言算法1 差商及其牛顿插值公式1.1 差商及其主要性质定义若已知函数f (x) 在点 xi (i0,1,2, , n) 处的函数值f ( xi ) 。则称：f x0 f (x0 ) 为函数 f (x) 在点 x0 的 0 阶差商；f x0 , x1

2、f x0 f x1 为函数 f ( x) 过点 x0 , x1 的 1阶差商；x0x1f x0 , x1 , x2 f x0 , x1f x0 , x2 为函数 f (x) 过点 x0, x1 , x2 的 2 阶差商；x1x2以此类推，一般地称f x0 , x1 , x2 , xk1, xk f x0 , x1 , xk 2 , xk 1 f x0 , x1 , xk 2 , xk 为函数 f (x) 过xk 1xk点 x0 , x1, xk 的 k 阶差商。性质 1阶差商表示为函数值f ( x0 ), f ( x1 ), , f ( xk ) 的线性组合。即kf (x j )f x0 ,

3、x1 , x2 , , xk =( x j x j 1)( x jx j 1)( xj xk )j 0 ( x jx0 )k f ( xj )j 0k( x jxi )i 0i j.*性质2 若函数f ( x) 在包含节点的区间 a, b 上存在 k 阶导数，则 k 阶差商与导数的关系为f x0 , x1 , , xk f (k ) ( ) , a, bk !1.2 牛顿插值公式通过 n1个互异点上的次数不超过n 的插值多项式Pn ( x) 可以表示如下形式：Pn ( x) f x0 f x0 , x1 ( x0x1)f x0 , x1 , x2 ( x x0 )( x x1)f x0 , x

4、1, , xn ( xx0 )( x xn 1 )这种形式的插值多项式称为牛顿插值多项式，一般记为N n(x)f x0 f x0 , x1 ( x0x1 )f x0 , x1, x2 ( xx0)( xx1 )f x0 , x1 , xn ( xx0 )( xxn 1 )由牛顿插值多项式可以看出， 当增加一个插值点时， 当前已有的各项不变， 只需要在后面增加一项即可。并且， 在牛顿插值公式中， 每一项的系数就是各阶差商， 比拉格朗日插值公式计算量小，且便于程序设计。根据差商性质2，即f x, x , x f (n 1) ( ) , a,b01n( n 1)!就可以将拉格朗日插值公式的余项转化成

5、牛顿插值公式的余项，即Rn ( x)f ( n 1) ()(n1)!f x0 , x1, , xn 牛顿插值公式余项更具有一般性，它对于列表函数或者f (x) 导数不存在的情形都适用。2 差分与等距结点插值公式2.1 差分及其主要性质定义 设函数 yfx 在等距结点xkx0kh, k0,1, n上的函数值为fkfxk, k0,1, n其中，为常熟，h 称为步长。则.*f (xk )f ( xk 1 )f ( xk ) 称为函数 f (x) 在处以 h 为步长的一阶向前差分， 并简记为f kfk 1fk ；f ( xk )f ( xk )f (xk 1 ) 称为函数 f ( x) 在处以 h 为

6、步长的一阶向后差分，并简记为f kfkf k 1f ( xk )f ( xk 1/2 )f ( xk 1/2 ) 称为函数 f ( x) 在处以 h 为步长的一阶向中心差分，并简记为fkf k 1/2f k1/2性质 1各阶差分可以表示成函数值的线性组合。即(1) nn1)k Cnk f j n kf i(k 0 n(2)n f j(1)kCnkf jkk0n1)kCnk f(3)f j(nkk 0j2性质 2 差商与差分具有如下关系：nf0n fnf x0, x1, , xk n!hnn!hn2.2 等距结点插值公式2.2.1 牛顿前插公式N n ( x)f 0t f0 +t (t 1)2f

7、0t (t1) (tn1) nf0！n!2nt(t1)(tk1)k f 0k 0k !其余项公式为Rn ( x)Rn ( x0th)f (n1)( )(tn +1 x0, xn (n！ t(t 1)n) h，1)2.2.2 牛顿后插公式N n ( x) f nt(t1)2f n( 1)n t (t 1)(t n 1) nfnt fn +！n!2n1)k t(t1) (tk1) k f n(k 0k !.*其余项公式为Rnf ( n 1) ()n 1(t n)hn +1, xn (x) Rn (xn th)t (t 1)( 1)， x0(n 1)！在用牛顿前插和后插公式计算时候，要涉及到各阶前插

8、和后插计算，下面是各阶向前和向后差分的计算格式，如下图所示。表 1 各阶向前差分和向前差分的计算公式xkfk1 阶差分2 阶差分3 阶差分4 阶差分用于前插公式x0f0f02 f03 f04 f0gggx1f1f12 f3 f1gggx2f2f22 f2gggx3f3f3gggx4f4gggggggggxn4f n4gggxn3f n3fn 3gggxn 2f n 2fn 22 fn 2gggxn 1f n 1fn 12 f n 13 fn 1gggxfnfn2 fn3 fn4 fn用于后插公式n3 牛顿差值公式的C 程序设计和应用实例3.1 牛顿差值法的应用步骤步骤 首先我们按照表1，求得各

9、点的差商. 然后利用牛顿前差或后差公式，把数值带入 .即可以求得n 次多项式。它在计算机上的应用步骤如下：步骤 1 输入所要求的牛顿多项式的次数n ，并依次输入n1 个节点 ( xi , yi ) .for : i=0ton+1.*scanf("%f",&xi)scanf("%f",&y0i);步骤 2 计算各阶差商for : i=1ton+1for : j=iton+1if(i>1)yij=(yi-1j-yi-1j-1)/(xj-xj-i);elseyij=(yi-1j-yi-1j-1)/(xj-xj-1);步骤 3 代入牛顿插值

10、公式, 可计算得出结果。for： i=1 ton+1temp=temp*(xx-xi-1);牛顿 =牛顿 +yii*temp;3.2 利用牛顿插值法程序的实例为了更方便的应用牛顿插值法，我们进行了与计算机的结合，下面我们将展示几个例子。例 2.1 已知 ysin x 的一组数据为x06432sinx01231222（1）构造牛顿插值函数并作图分析。（2）并分别利用程序估计 sin， sin14 的估计值。1215分析首先我们可以通过程序求出差商表：xkf (xk )一阶差商二阶差商三阶差商四阶差商.*0010.9549292686220.7910896312.0860764230.607024

11、4240.351538650.1364890932210.255872630.447100350.091254700.028797106带入定义1.5 中可求得如下牛顿插值多项式如下N 2( x)0.95492926x0.2086076x( x)(2.1)6N 3( x) 0.95492926x0.2086076x(x) 0.13648909x( x)( x)(2.2)664N 4( x)0.95492926 x 0.2086076 x( x)0.13648909 x(x)( x)6640.028797106 x( x)( x)( x)643第二步 利用 C+程序计算 sin和 sin 14的

12、值。步骤如下：利用C 程序：首先输入121514 的所要求的牛顿多项式的次数n，然后输入n+1 个节点的值 . 即可以得出 sin和 sin值为 0.2586 和 0.0804 ；1215例 2.2 设 f ( x) ln(1 x) 的函数表如下：x0.250.300.360.390.45f(x)0.2231440.2623640.3074850.3293040.371564试计算 ln(1.275) ， f ( x)ln 2.3分析：同上题步骤我们先求差商表，进而代入公式可得N4 ( x)0.22314 0.7844(x0.25)0.294393939(x0.25)( x0.30)0.141

13、1736357(x0.25)( x0.30)( x0.36)0.057720032( x 0.25)( x 0.30)( x0.36)( x0.39)N3 ( x)0.22314 0.7844(x0.25)0.294393939(x0.25)( x0.30)0.1411736357(x0.25)( x0.30)( x0.36).*N2( x) 0.22314 0.7844(x 0.25)0.294393939 (x0.25)( x 0.30)利用 C 程序我们可以得到计算结果：N4(0.275)0.242946，N3 (0.275) 0.242945，N2(0.275)0.242938N 4

14、(1.3)0.928823，N 3 (1.3) 0.877002，N2(0.275)0.737653从上述两个例子我们可以看出, 多项式在区间0, 周围与原函数逼近的较好。离这个4区间越远与原函数的误差越大在x1.5 处时，该图像就已经开始背离图像了. 所以该多项式只能在一个小的区间里可以逼近原函数, 不适合作为原函数的逼近函数. 也可以看出多项式在 0, 区间的周围逼近的较好，但是x2处时，该图像就离原图像误差较大. 多项式在区3间0,2.5 都逼近的挺好，从图中我们看出在远离这个区间的图像误差相对较大，但是在这三个多项式中是逼近最好的。于是可以得出节点越多，函数逼近的相对较好. 在节点附近

15、逼近的越好，越远离节点误差越大， 所以公式适用于计算节点附近的值于是为了减小误差， 在下一节的等距节点下的插值公式根据所求的点的函数值的不同分别采取了前插和后插公式。3.4 等距节点下的牛顿插值算法与程序设计前面我们讲述了一般节点下的牛顿插值公式，为了计算方便于是有了对等距节点下的牛顿多项式的研究，本节将对等距节点下的插值多项式进行总结讨论.3.4.1 等距节点下的牛顿插值法的程序算法步骤步骤先求差分. 然后利用牛顿前插公式或牛顿后插公式并把数值带入. 即可以求得n 次多项式。它在计算机上的应用步骤如下：步骤1 输入所要求的牛顿多项式的次数n ，步长h ，并依次输入n1个节点( xi , yi

16、 )。for : i=0ton+1scanf("%f",&xi)scanf("%f",&y0i);步骤2 求得各界差分for : i=1ton+1 for : j=iton+1yij=(yi-1j+1-yi-1j-1)/(xj-xj);.* / 求向前差分for : i=1ton+1 for : j=iton+1yij=(yi-1j-yi-1j-1)/(xj-xj-i);/求向后差分步骤 3 代入牛顿插值公式, 可计算得出结果for(i=1;i<n+1;i+)temp=temp*(t-i+1)/i);牛顿 =牛顿 +yii*temp

17、;printf(" 求得的结果为:N(%.4f)=%9fn",xx,牛顿 );/ 求得运用前插公式的值for(i=n;i<0;i-)temp=temp*(t-i+1)/i);牛顿 =牛顿 +yn-in-i*temp;printf(" 求得的结果为 :N(%.4f)=%9fn",xx, 牛顿 );/求得运用后插公式的值3.4.2等距节点下牛顿插值的实例例 已知 tan x 的值列表如下：x1. 31.311.321.33tan x3.60213.74713.90334.0723近似计算tan1.325tan1.305,。采用牛顿向后插公式 . 为此，

18、做差分表xifif i2 fi3 fi1.33.60210.14500.01120.00161.313.74710.15620.01281.323.90330.16901.334.0723从而，有p3 (x)4.0723 0.1690 t0.01280.00162).t (t1)t (t 1)(t1.3251.332!3!将 t0.5 代入上式，得0.01tan1.325p3 (1.325)3.9869。.*将 t1.305 1.332.5带入后插多项式中可以得到0.01tan1.305p3 (1.305)3.0797现在我们利用C 程序首先输入所求插值的次数3 和步长0.01. 然后输入各个

19、节点，并输入所要求的点1.325 既可以求出该点的函数值。即f (1.325) 3.9869 。 f (1.26)p3 (1.26)3.0797例设函数在各节点的取值如下x00.20.40.60.81.0f ( x)1.00.8187310.6703200.5488120.4493290.367879试利用插值公式求f (0.3) 的值。采用牛顿向前插公式，同上题我们先做差分表，然后相应带入到差分公式. 中求得后插公式。利用 C 语言程序步骤如下:首先输入所求插值的次数5 和步长0.2。然后输入各个节点，并输入所要求的点0.3 既可以求出该点的函数值。即f (0.3)0.7880 。由以上例子

20、我们看到例1 用了牛顿后插公式，例 2 用了牛顿前插公式， 我们该怎样选取。这个经过验证得出，如果所要求的点较靠近节点x1 ，则采用前插公式；如果靠近xn ，则采用牛顿后插公式。4 结束语用牛顿差值方法处理测量数据，具有使差值多项式通过选定测量值的特点，所以在数据处理中有一定的应用场合。当测量值较多、较密时，为了减轻计算工作量及提高准确性，应选取合适的测量值作为差商计算的依据求得N n ( x) 。参考文献 1李庆扬等编 . 数值分析 . 华中科技大学出版社 , 1989 2储钟武等编译 . 数值分析 . 黑龙江科学技术出版社 , 1987 3徐士良编 . 数值分析与算法 . 机械工业出版社

21、, 2007附录 A ：牛顿插值法的程序#include<stdio.h>Newtonvoid main()float x11,y1111,xx,temp,牛顿 ;.*int i,j,n;printf(" 牛顿插值 :n 请输入要运算的值:x=");scanf("%f",&xx);printf(" 请输入插值的次数(n<11):n=");scanf("%d",&n);printf(" 请输入 %d 组值 :n",n+1);for(i=0;i<n+1;i+)

22、printf("x%d=",i); scanf("%f",&xi);printf("y%d=",i);scanf("%f",&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;牛顿 =y00;for(i=1;i<n+1;i+)temp=te

23、mp*(xx-xi-1);牛顿 =牛顿 +yii*temp;printf("求得的结果为:N(%.4f)=%9fn",xx,牛顿 );附录 B：等距节点下的前插公式的C 语言程序#include<stdio.h>void main()float x11,y1111,xx,temp,牛顿 ,t,h;.*int i,j,n;printf(" 牛顿插值 :n 请输入要运算的值:x=");scanf("%f",&xx);printf(" 请输入插值的次数(n<11):n=");scanf(&quo

24、t;%d",&n);printf(" 步长 :n 请输入要运算的值:h=");scanf("%f",&h);printf(" 请输入 %d 组值 :n",n+1);for(i=0;i<n+1;i+)printf("x%d=",i);scanf("%f",&xi);printf("y%d=",i);scanf("%f",&y0i);t=(xx-x0)/h;for(i=1;i<n+1;i+)for(j=i;j<n+1;j+)yij=(yi-1j+1-yi-1j);printf("%fn",yii);temp=1;牛顿 =y00;for(i=1;i<n+1;i+)temp=temp*(t-i+1)/i);牛顿 =牛顿 +yii*temp;printf(" 求得的结果为 :N(%.4f)=%9fn",xx,牛顿 );附录 C：等距节点下的后插公式的C 语言程序#include<stdio.h>void main()float x11