c++实现牛顿插值法实验报告

1、用newton商差公式进行插值姓名：陈辉学号：13349006院系：数据科学与计算机学院专业:计算机科学与技术班级：计科-班日期：2015-10-11指导老师：纪庆革目录一、实验目的3二、实验题目3三、实验原理与基础理论3四、实验内容6五、实验结果10六、心得体会14七、参考资料14八、附录（源代码）14一、实验冃的编写一个程序，用牛顿差商公式进行插值。二、实验题目编写一个程序，用牛顿差商公式进行插值。三、实验原理与基础理论牛顿插值公式为:nnm = fm + /xoik% 一 %0） + + fxof .,xn（x 一 %0）（x - xn_!）其中，丫_（衍）-心）/ 00*兀1一 衍-%

2、0耐/!巾耳一门兀0，尢订fxqf.9xk =xk 尤0我们将从键盘读入n阶牛顿插值的n+1个节点、= 0,1, 以此得 出牛顿插值多项式。有了节点，我们只需要求亢勺，,冷即可。我们记fxmf., s订为tm k,则tmk在差商表表的位置为（m, k）:mk012 n0/（%o）2心）fxlfx03心）f込衍fx2txlfx0 nfn> 兀n-lf 氏pn-2 f 兀？v > x°由f ro，,g的公式可得 tmk = tmk-l - tm-lk-l) / (xi直观上来看，表中的某格的差商值口j由其左边和左上边的值求得，从而牛顿插值 多项式变为 n(x) = t00 +

3、 tll(x-xo) + + tnn(x-xo).(x-xn-l)做到这一步，我们已经可以通过上面的数据计算对于给出的x, f(x)的近似 值。为了更规范地表示，下面我把n(x)变换成标准的降幕多项式n(x)二anxan + an-lxa(n-l) + . + a2xa2 + alx + a0o为了便于运算,不妨先把xxi中的减号换成加号(在最后令yi=-xi,在 令xi=yi,仍可以得到原本的结果),那么有:(% + xo)(x + xl) . (% + xm 1)m-1=xm + xm_1xii=0+ xm2 xioxil + + x°xioxil .xim 1为了便于表示，把x

4、toxil.% ik 1记为m那么(% + xo)(x + xl) . (% + xm 1)只要把n(x)中的每一项展开然后x次数相同的系数相加就可以得到数组ao于是可以列出下表：x0xl xk xn-lxnt010000tl11000t222000 tkxn k n_kxn k 1n-k10 tn-lx n 1n-1xn 2n-1 xn - k 1n-1 10tnxnnxn 1n xn kn »1n1表中xi列的和就是n(x)中ai的值，下面就来求这个和，记cgh= 工 xi0 .xih-l = xhcgh的意义为g个数屮所冇h个数乘积z和，可以由g-1个数中所冇h-1 个数乘积z

5、和，和g-1个数中所有h个数乘积z和求得，递推公式为 cgh=cg-lh-lxih+cg-lh o由cg h的意义可以得到递推的边界状态 为 ci0=x0+.+xi, cii=x0.xio于是我们又可以得到一张数组c的表（初始状态）:ij01 kn0x01x0+xlx0xl kk冰t=lkt=lnnt=lnt=l表的左下角每个空格都可以由其上血的和左上边的格中的值求;ii。这样，所需要求的所有值都已经得到，只需要通过简单的循环结构就可以算 出数组a,从而得到一个降幕的多项式。四、实验内容实验环境:mac os x yosemite, sublime text 2, xcode, cmd实验步骤

6、：我用一个类封装牛顿插值算法:678901234567890111122222222223private:int n;double x20r f20r a20, c2020f y20;double t2020; / divided difference table public:neonlnterpolation();void interpolation。；void input();void maketableo;double divideddifference(double flr double f0r double xlr double x0); void printtableo;void

7、 calculate!);void makeformula(); void formulao;；newtoninterpolation 0 为类的初始化:void nexrtonlnterpolation: :interpolation()interpolation()类似于类中的主函数，依此调用下面其他的函数:34me(nset(xfsizeofx);35metnsekf,0rsizeoff)；36mefnset(at0rsizeofa)；37me(nset(cf0rsizeofc);38me<nset(tfsizeoft);39414243444546474849inputt);ma

8、ketablet); printtable()j calculate();makeformula(); fonmula():inputo从键盘读入插值多项式的次数和节点，xi和fi为节点，yi和ti 将会在后面解释，这里对其进行初始化：51525354555657585960616263void newtoninterpolation: input()cout « "input the degree of newton interpolation: 11 « endl; cin » n;cout « hinput 11 « rv+1

9、« h points: x f(x)11 « endl; for (int 1 = i0； i <= nj i 卄)cin » xi » fi】； yi = -xi;ti 0 - fil;>cout « endl;maketableo用斧商表中tmk的递推式求出数组t所有的元素，这里的t 已经在input。中进行过初始化：65void newtoninterpolation:maketableo6667for (int i = 1; i <= n; i +)68for (int j = 1; j <= n; j +)69

10、70if (j > i) continue;71tij = divideddifference(tij-1r ti-lj-lr xilt xi-j);7273dividedifference ()为求值公式:double nextonlnterpolation:divideddifference(double flr double f0t double xlf double x0) 76double ans = (fl-f0) / (xl-x0);78 return ans;79 printtableq输出 fxo,xl到 fxoxn的值:void newtoninterpolation

11、:printtableofor (int i = 1; i <= n; i +)cout « ufx0h； for (int j = 1； j <= i； j +) cout «« j;cout « h = 11 « t 1 i « endl;cout « hn(x) = f(x0) + fx0rxl(x-x0) + . fx0t.fxn(x-x0).(x-x.n-1)11 « endl; 93949596979899100101102103104calculateo计算第二部分中最后一张表c数组的值，

12、其中y屮的元素等于x 中对应元素的相反数，这步操作已在input。中完成：void newtoninterpolation:calculate() c00 y0j;for (int i = 1; i <= n; i +) ci0 = ci-l0 + yi;for (int i=l; i <= n; i +) cii = ci-li-1 * yi;ci j = ci-l j-l*yi + ci-l j;for (int j = 1; j for (int i =asmakeformulao运用第二部分中的第二张表，每列累加计算出标准降幕多项 式中的系数数组a,并输出到屏幕上：3060

13、7080910111213141516171819 allllllllllllllvoid nexrtonlnterpolation: makefonnulat)for (int j = 0; j <= n; j +) for (int i = j; i <= n; i +)if (i = j) aj = aj elseaj = aj+ ti i*ci-l i-j-1;cout « "n(x) = 11; for (int i = n; i >= 0;0 1222222211111112 3 4 5 6-cout « ai】； if (i !=

14、0)cout « "x11 « cout « endl; cout « endl;formulaq显示捉示并从键盘读入x,然后用整理后的插值多项式计算对应的 近似值f(x)：12b void newtoninterpolation:fonnula()129 130 double xr ans = q;131 cout « "input x ="132 cin » x;133 for (int i = 0; i n; i +)134 ans = ans + ai*pow(xr i):135 cout

15、71; "nc « x « pl) = n « ans « endl;136 1"程序主函数很简单，创建类并调用函数:138 int main(int argct const char * argv)139 140 netrtonlnterpolation newton;141 netrton.interpolation();142 >由于我是在苹果系统卜完成这个实验，生成的可执行文件只能在mac系统 下打开，所以最后我在win虚拟机中编译了一个win下的可执行文件。五、实验结果第一组：我用练习第一题中的主句进行检验。提示输入多

16、项式次数：r1 c)仝 chanfai newtoninterpolation newtonlnterpolat 80x24last login: mon oct 12 00:02:09 on ttys000/users/chanfai/desktop/study/2015-2/nmc/lab/newtoninterpolation/build/debug/newt onlntepolation ; exit;chanfais-macbook-air: chanfai$ /users/chanfai/desktop/study/2015-2/nmc/lab/newt onlnterpolati

17、on/build/debug/newtontnterpolation ; exit;input the degree of newton interpolation:提示输入节点：全 chanfai newtonlnterpolation newtonlast login: mon oct 12 00:02:09 on ttys000 /users/chanfai/desktop/study/2015-2/nmc/lab/newtoninterpolation/build/debug/newt onlnterpolation ; exit;chanfais-macbook-air:<v

18、chanfai$ /users/chanfai/desktop/study/2015-2/nmc/lab/newt onlnterpolation/build/debug/newtonlnterpolation ; exit;input the degree of newton interpolation:2input 3 points: x f(x)计算出fxo,xl到fxo,xn的值，和将此多项式的系数并显示，这个结果和 我手算得结杲一致，然后提示输入x：r1 g)佥 chanfai newtoninterpolation newtonlnterpolat 80x24last login:

19、 mon oct 12 00:02:09 on ttys000 /users/chanfai/desktop/study/2015-2/nmc/lab/newtoninterpolation/build/debug/newt onlnterpolation ; exit;chanfals-macbook-air: chanfai$ /users/chanfai/desktop/study/2015-2/nmc/lab/newt onlnterpolation/build/debug/newtonlnterpolation ; exit;input the degree of newton in

20、terpolation:2input 3 points: x f(x)6 1ii 1；1 2 fx0fxl 1fx0fxl,x2 = 0.5)；n(x) = f(x0) + fx0fxl(x-x0) + . + fx0f.fxn(x-x0).(x-x.n-1) n(x) m 0.5x2 + 0.5xl + 1input x =输入xz后计算出f(x)的近似值，然后程序结束： (金 chanfai newtoninterpolation 80x24last login: mon oct 12 00:02:09 on ttys000 /users/chanfai/desktop/study/201

21、5-2/nmc/lab/newtoninterpolation/build/debug/newt onlnterpolation ; exit;chanfais-macbook-air:<v chanfai$ /users/chanfai/desktop/study/2015-2/nmc/lab/newt onlnterpolation/build/debug/newtonlnterpolation ; exit;input the degree of newton interpolation:input 3 points: x f(x)1 22 4fx0,xl = 1n(x) * f(

22、x0) + fx0fxl(x-x0) + . + fx0f.fxn(x-x0).(x-x.n-1)n(x)工 0.5x2 + 0.5x1 + 1input x = 3n(3) = 7logout程序完成】i第二组：笫一组是二阶的牛顿插值多项式，笫二组我用练习题中第五题来测试，这题 冇五个节点，是四阶牛顿插值多项式。运行结杲如下，和我手算得结杲一致：金 chanfai newtoninterpolation 80x26last login: mon oct 12 00:15:13 on ttys000chanfais-macbook-air: chanfai$ /users/chanfai/de

23、sktop/study/2015-2/nmc/lab/newt onlnterpolation/build/debug/newtonlnterpolation ; exit;input the degree of newton interpolation:4newtonlnteroolationinput 5 points: x f(x)1.615 2.414501.634 2.464591.702 2.652711.828 3.030351.921 3.34066fx0.x1 = 2.63632fx0fxltx2 = 1.49603fx0rxlfx2fx3 m -1-44131fx09xl9

24、x2fx39x4 8.82423n(x) = f(xo) + fx0fxl(x-x0) + . + fx0f.fxn(x-x0).(x-x.n-1)n(x) = 8.82423x4 + -61.2608x3 + 160.578x2 + -185.398x1 + 81.0281input x = 1.682n(l.682) = 2.59612logout程序完成】i金 chanfai newtoninterpolation 80x26last login: mon oct 12 00:18:26 on ttys000chanfais-macbook-air: chanfai$ /users/ch

25、anfai/desktop/study/2015-2/nmc/lab/newt onlnterpolation/build/debug/newtonlnterpolation ; exit;input the degree of newton interpolation:4newtonlnteroolationinput 5 points: x f(x)1.615 2.414501.634 2.464591.702 2.652711.828 3.030351.921 3.34066fx0.x1 = 2.63632fx0fxltx2 = 1.49603fx0rxlfx2fx3 m -1-4413

26、1fx09xl9x2fx39x4 8.82423n(x) = f(xo) + fx0fxl(x-x0) + . + fx0f.fxn(x-x0).(x-x.n-1)n(x) = 8.82423x4 + -61.2608x3 + 160.578x2 + -185.398x1 + 81.0281input x = 1.813n(l.813) = 2.98332logout程序完成】i六、心得体会在完成这个程序z前，我进行了大量的计算，特别是把插值形式的式子整理 成降幕的多项式那几个步骤。最后编写完成的时候发现结杲怎么也不止确，然后 我调试后发现是一个列表的行和列弄反了，改正这个错误之后就正确了。在

27、做这 个实验的过程屮，我对很多数学的方法有了新的领会和心得，对word中插入公 式的操作也熟悉了很多。完成程序后我在百度和谷歌上搜索了一下，还没有一个 搜索结果是可以把多项式整理为降幕排列的。七、参考资料数值方法（c+描述）,pal lab ghosh著,徐士良译,清华大学出版社八、附录（源代码）/ main.cpp/ newtoninterpolation/ created by 陈辉 on 2015/10/11 / copyright （c） 2015 年 chanfai. all rights reserved./# include<iostream># include<

28、;cmath>using namespace std;classnewtonlnterpolationprivate:int n;double x20, f20, a20, c2020, y20;double t2020; / divided difference tablepublic:newtonlnterpolationq；void interpolation();void inputf)；voidmaketablef);doubledivideddifferencefdouble fl, double fo, double xl, double xo); voidprinttab

29、leq;void calculatef);voidmakeformulaf);void formula。；;newtoninterpolation:newtoninterpolationomemsetfx, 0, sizeof x);memsetft 0, sizeof f;memsetfa, 0, sizeof a);memset(c, 0, sizeof c);memsettt, 0, sizeof t);voidnewtoninterpolation:lnterpolation()inputq;maketable()；printtableo；calculateo；makeformulao

30、；formulao；voidnewtoninterpolation:inputocout« "input the degree of newton interpolation:11 «endl;cin» n;cout« "input" « n+1 « " points: x f(x)" «endl;for (inti = 0; i<= n; i +)cin» xi » fi;yi=xi;tio=fi；cout«endl;voidnewton

31、interpolation:maketableofor (inti = 1; i<= n; i +)for (int j = 1; j <= n; j +)if (j >i) continue;ti j = divideddifference(ti j-1, ti-l j-1, xi, xi-jj;doublenewtoninterpolation:divideddifference(double fl, double fo, double xl, double xo)doubleans = (fl-fo) / (xl-xo); returnans;voidnewtoninterpolation:printtable()for (