c++實(shí)現(xiàn)牛頓插值法實(shí)驗(yàn)報(bào)告

1、數(shù)值實(shí)驗(yàn)用Newton商差公式進(jìn)行插值姓名：陳輝學(xué)號：13349006院系：數(shù)據(jù)科學(xué)與計(jì)算機(jī)學(xué)院專業(yè)：計(jì)算機(jī)科學(xué)與技術(shù)班級：計(jì)科一班日期：2015-10-11指導(dǎo)老師：紀(jì)慶革- 2 -目錄一、實(shí)驗(yàn)?zāi)康?二、實(shí)驗(yàn)題目3三、實(shí)驗(yàn)原理與基礎(chǔ)理論3四、實(shí)驗(yàn)內(nèi)容6五、實(shí)驗(yàn)結(jié)果10六、心得體會(huì)14七、參考資料14八、附錄（源代碼）14- 19 -一、 實(shí)驗(yàn)?zāi)康木帉懸粋€(gè)程序，用牛頓差商公式進(jìn)行插值。二、 實(shí)驗(yàn)題目編寫一個(gè)程序，用牛頓差商公式進(jìn)行插值。三、 實(shí)驗(yàn)原理與基礎(chǔ)理論牛頓插值公式為：Nnx=fx0+fx0,x1x-x0+fx0,xnx-x0(x-xn-1)其中，fx0,x1=fx1-f(x0)x1-x

2、0fx0,xk=fx1,xk-fx0,xkxk-x0我們將從鍵盤讀入n階牛頓插值的n+1個(gè)節(jié)點(diǎn)xi, fi,i=0,1,n，以此得出牛頓插值多項(xiàng)式。有了節(jié)點(diǎn)，我們只需要求fx0,xk即可。我們記fxm,xm-1,xm-k為tmk，則tmk在差商表表的位置為(m, k)：m k012n0f(x0)2f(x1)fx1,x03f(x2)fx2,x1fx2,x1,x0nf(xn)fxn,xn-1fxn,xn-1,xn-2fxn,x0由fx0,xk的公式可得tmk = (tmk-1 - tm-1k-1) / (xi xi-j)，直觀上來看，表中的某格的差商值可由其左邊和左上邊的值求得，從而牛頓插值多項(xiàng)式

3、變?yōu)镹(x) = t00 + t11(x-x0) + + tnn(x-x0)(x-xn-1)做到這一步，我們已經(jīng)可以通過上面的數(shù)據(jù)計(jì)算對于給出的x，f(x)的近似值。為了更規(guī)范地表示，下面我把N(x)變換成標(biāo)準(zhǔn)的降冪多項(xiàng)式N(x) = anxn + an-1x(n-1) + + a2x2 + a1x + a0。為了便于運(yùn)算，不妨先把x-xi中的減號換成加號（在最后令yi=-xi，在令xi=yi，仍可以得到原本的結(jié)果），那么有：x+x0x+x1x+xm-1=xm+xm-1i=0m-1xi+xm-20i0i1m-1xi0xi1+x00i0im-1m-1xi0xi1xim-1為了便于表示，把0i0i

4、km-1xi0xi1xik-1記為mxk那么x+x0x+x1x+xm-1=xm+xm-1mx1+x0mxm只要把N(x)中的每一項(xiàng)展開然后x次數(shù)相同的系數(shù)相加就可以得到數(shù)組a。于是可以列出下表：x0x1xkxn-1xnt010000t11x11000t22x12x1000tkn-kxn-kn-kxn-k-110tn-1n-1xn-1n-1xn-2n-1xn-k-110tnnxnnxn-1nxn-knx11表中xi列的和就是N(x)中ai的值，下面就來求這個(gè)和，記cgh=0i0ihg-1xi0xih-1=gxhcgh的意義為g個(gè)數(shù)中所有h個(gè)數(shù)乘積之和，可以由g-1個(gè)數(shù)中所有h-1個(gè)數(shù)乘積之和，和

5、g-1個(gè)數(shù)中所有h個(gè)數(shù)乘積之和求得，遞推公式為cgh=cg-1h-1xih+cg-1h。由cgh的意義可以得到遞推的邊界狀態(tài)為ci0=x0+xi，cii=x0xi。于是我們又可以得到一張數(shù)組c的表（初始狀態(tài)）：i j01kn0x01x0+x1x0x1kt=1kxtt=1kxtnt=1nxtt=1nxt表的左下角每個(gè)空格都可以由其上面的和左上邊的格中的值求出。這樣，所需要求的所有值都已經(jīng)得到，只需要通過簡單的循環(huán)結(jié)構(gòu)就可以算出數(shù)組a，從而得到一個(gè)降冪的多項(xiàng)式。四、 實(shí)驗(yàn)內(nèi)容實(shí)驗(yàn)環(huán)境：Mac OS X Yosemite，Sublime Text 2，Xcode，CMD實(shí)驗(yàn)步驟：我用一個(gè)類封裝牛頓插

6、值算法：NewtonInterpolation()為類的初始化：Interpolation()類似于類中的主函數(shù)，依此調(diào)用下面其他的函數(shù)：Input()從鍵盤讀入插值多項(xiàng)式的次數(shù)和節(jié)點(diǎn)，xi和fi為節(jié)點(diǎn)，yi和ti將會(huì)在后面解釋，這里對其進(jìn)行初始化：MakeTable()用差商表中tmk的遞推式求出數(shù)組t所有的元素，這里的t已經(jīng)在Input()中進(jìn)行過初始化：DivideDifference()為求值公式：PrintTable()輸出fx0,x1到fx0,xn的值：Calculate()計(jì)算第二部分中最后一張表c數(shù)組的值，其中y中的元素等于x中對應(yīng)元素的相反數(shù)，這步操作已在Input()中完成

7、：MakeFormula()運(yùn)用第二部分中的第二張表，每列累加計(jì)算出標(biāo)準(zhǔn)降冪多項(xiàng)式中的系數(shù)數(shù)組a，并輸出到屏幕上：Formula()顯示提示并從鍵盤讀入x，然后用整理后的插值多項(xiàng)式計(jì)算對應(yīng)的近似值f(x)：程序主函數(shù)很簡單，創(chuàng)建類并調(diào)用函數(shù)：由于我是在蘋果系統(tǒng)下完成這個(gè)實(shí)驗(yàn)，生成的可執(zhí)行文件只能在Mac系統(tǒng)下打開，所以最后我在win虛擬機(jī)中編譯了一個(gè)win下的可執(zhí)行文件。五、 實(shí)驗(yàn)結(jié)果第一組：我用練習(xí)第一題中的主句進(jìn)行檢驗(yàn)。提示輸入多項(xiàng)式次數(shù)：提示輸入節(jié)點(diǎn)：計(jì)算出fx0,x1到fx0,xn的值，和將此多項(xiàng)式的系數(shù)并顯示，這個(gè)結(jié)果和我手算得結(jié)果一致，然后提示輸入x：輸入x之后計(jì)算出f(x)的近似

8、值，然后程序結(jié)束：第二組：第一組是二階的牛頓插值多項(xiàng)式，第二組我用練習(xí)題中第五題來測試，這題有五個(gè)節(jié)點(diǎn)，是四階牛頓插值多項(xiàng)式。運(yùn)行結(jié)果如下，和我手算得結(jié)果一致：六、 心得體會(huì)在完成這個(gè)程序之前，我進(jìn)行了大量的計(jì)算，特別是把插值形式的式子整理成降冪的多項(xiàng)式那幾個(gè)步驟。最后編寫完成的時(shí)候發(fā)現(xiàn)結(jié)果怎么也不正確，然后我調(diào)試后發(fā)現(xiàn)是一個(gè)列表的行和列弄反了，改正這個(gè)錯(cuò)誤之后就正確了。在做這個(gè)實(shí)驗(yàn)的過程中，我對很多數(shù)學(xué)的方法有了新的領(lǐng)會(huì)和心得，對word中插入公式的操作也熟悉了很多。完成程序后我在百度和谷歌上搜索了一下，還沒有一個(gè)搜索結(jié)果是可以把多項(xiàng)式整理為降冪排列的。七、 參考資料數(shù)值方法（C+描述），P

9、ALLAB GHOSH 著，徐士良 譯，清華大學(xué)出版社八、 附錄（源代碼）/ main.cpp/ NewtonInterpolation/ Created by 陳輝 on 2015/10/11./ Copyright (c) 2015年 CHANFai. All rights reserved./#include <iostream>#include <cmath>using namespace std;class NewtonInterpolationprivate: int n; double x20, f20, a20, c2020, y20; double t

10、2020; / Divided Difference Tablepublic: NewtonInterpolation(); void Interpolation(); void Input(); void MakeTable(); double DividedDifference(double f1, double f0, double x1, double x0); void PrintTable(); void Calculate(); void MakeFormula(); void Formula();NewtonInterpolation:NewtonInterpolation()

11、 memset(x, 0, sizeof x); memset(f, 0, sizeof f); memset(a, 0, sizeof a); memset(c, 0, sizeof c); memset(t, 0, sizeof t);void NewtonInterpolation:Interpolation() Input(); MakeTable(); PrintTable(); Calculate(); MakeFormula(); Formula();void NewtonInterpolation:Input() cout << "Input the de

12、gree of Newton Interpolation: " << endl; cin >> n; cout << "Input " << n+1 << " points: x f(x)" << endl; for (int i = 0; i <= n; i +) cin >> xi >> fi; yi = -xi; ti0 = fi; cout << endl;void NewtonInterpolation:MakeTable

13、() for (int i = 1; i <= n; i +) for (int j = 1; j <= n; j +) if (j > i) continue; tij = DividedDifference(tij-1, ti-1j-1, xi, xi-j); double NewtonInterpolation:DividedDifference(double f1, double f0, double x1, double x0) double ans = (f1-f0) / (x1-x0); return ans;void NewtonInterpolation:P

14、rintTable() for (int i = 1; i <= n; i +) cout << "fx0" for (int j = 1; j <= i; j +) cout << ",x" << j; cout << " = " << tii << endl; cout << "N(x) = f(x0) + fx0,x1(x-x0) + . + fx0,.,xn(x-x0).(x-x.n-1)" << e

15、ndl;void NewtonInterpolation:Calculate() c00 = y0; for (int i = 1; i <= n; i +) ci0 = ci-10 + yi; for (int i = 1; i <= n; i +) cii = ci-1i-1 * yi; for (int j = 1; j <= n-1; j +) for (int i = j+1; i <= n; i +) cij = ci-1j-1*yi + ci-1j;void NewtonInterpolation:MakeFormula() for (int j = 0;

16、 j <= n; j +) for (int i = j; i <= n; i +) if (i = j) aj = aj + tii; else aj = aj + tii*ci-1i-j-1; cout << "N(x) = " for (int i = n; i >= 0; i -) cout << ai; if (i != 0) cout << "x" << i << " + " cout << endl; cout << endl;void NewtonInterpolation:Formula() double X, ans = 0; cout << "Input x = "