1、數值實驗用Newton商差公式進行插值姓名：陳輝學號：13349006院系：數據科學與計算機學院專業：計算機科學與技術班級：計科一班日期：2015-10-11指導老師：紀慶革- 2 -目錄一、實驗目的3二、實驗題目3三、實驗原理與基礎理論3四、實驗內容6五、實驗結果10六、心得體會14七、參考資料14八、附錄（源代碼）14- 19 -一、 實驗目的編寫一個程序，用牛頓差商公式進行插值。二、 實驗題目編寫一個程序，用牛頓差商公式進行插值。三、 實驗原理與基礎理論牛頓插值公式為：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個節點xi, fi,i=0,1,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)，直觀上來看，表中的某格的差商值可由其左邊和左上邊的值求得，從而牛頓插值多項式

3、變為N(x) = t00 + t11(x-x0) + + tnn(x-x0)(x-xn-1)做到這一步，我們已經可以通過上面的數據計算對于給出的x，f(x)的近似值。為了更規范地表示，下面我把N(x)變換成標準的降冪多項式N(x) = anxn + an-1x(n-1) + + a2x2 + a1x + a0。為了便于運算，不妨先把x-xi中的減號換成加號（在最后令yi=-xi，在令xi=yi，仍可以得到原本的結果），那么有：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)中的每一項展開然后x次數相同的系數相加就可以得到數組a。于是可以列出下表：x0x1xkxn-1xnt010000t11x11000t22x12x1000tkn-kxn-kn-kxn-k-110tn-1n-1xn-1n-1xn-2n-1xn-k-110tnnxnnxn-1nxn-knx11表中xi列的和就是N(x)中ai的值，下面就來求這個和，記cgh=0i0ihg-1xi0xih-1=gxhcgh的意義為g個數中所有h個數乘積之和，可以由g-1個數中所有h-1個數乘積之和，和

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

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

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

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

9、ALLAB GHOSH 著，徐士良 譯，清華大學出版社八、 附錄（源代碼）/ 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 = "