<a href="mailto:[email protected]">[email protected]</a>
Abstract. Power basis polynomial is the most simple polynomial function. It also be called power series. OpenCASCADE provides basic computation functions for polynomial functions, such as evaluate the result for a given polynomial, Lagrange interpolation, Hermite interpolation, .etc. The package named PLib, means Polynomial functions Library. The paper focus on the Lagrange interpolation usage of PLib.
Key Words. OpenCASCADE, PLib, Interpolation, Lagrange, 插值
1.Introduction
無窮級數是高等數學的一個重要組成部分，它是表示函數、研究函數性質及進行數值計算的一種工具。由高等數學中的無窮級數的概念可知，函數項級數中簡單而常見的一類級數就是各項都是幂函數的函數項級數即所謂的幂級數（Power Series）。因為幂級數的形式簡單，易于了解，且可以高效計算曲線上的點及各階導數，是以在幾何造型中經常用幂級數來近似表示曲線曲面。由于将函數展開成幂級數是有條件的，所有并不是所有的曲線曲面都可以用幂級數的多項式來逼近。
對無窮級數概念比較陌生的讀者可以把《高等數學》的書找出來翻翻看，重溫一下大學的時光。一打開做了筆記有點泛黃的舊書，就會回想起青澀的校園時光。
當時學習《高等數學》的時候感覺很抽象難了解，因為無法将理論與實踐聯系起來，看不到直覺效果，有時也會冒出“學數學有什麼用？”這種問題。當你遇到相關的問題再去看國内的教材時，覺得國内的教材寫得還是很細緻用心的。現在擷取資訊已經很便利了，如OpenCASCADE這個開源的庫，其TKMath工具箱可以看成是數值計算理論聯系實踐的一個具體執行個體。結合OpenCASCADE的源碼來對相關理論的學習，效率會事半功倍。
本文對幂級數的概念做簡單介紹，并結合源程式詳細說明OpenCASCADE中的PLib包中關于幂級數多項式的計算及Lagrange插值的用法。
2.Polynomial Evaluation
高等數學的書中把幂級數的重點落在如何将其他的函數展開成幂級數，即用幂級數來逼近函數，而沒有介紹如何用數值的方法來對幂級數進行計算。在算法導論一書[2]中找到相關多項式的表示及計算的實作方法，給出了多項式在資料結構上的表示方式及求值算法。對于一個幂級數多項式：
多項式的表示有系數表示法和點值表示法。系數表示法（Coefficient Representation）就是将多項式的系數組成一個向量來表示這個多項式。對用系數表示法表示的多項式的求值計算可以采用Horner法則，也是是著名的秦九韶算法。此算法的實作代碼在《The NURBS Book》[3]一書中給出了，此處略去，隻給出OpenCASCADE中對幂級數多項式的計算的函數的使用。
void testPolynomialEvaluation(void)
{
// evaluate 1 dimension polynmoial.
Standard_Real aCoeff[3] = {2.0, 2.0, 3.0};
Standard_Real aResult = 0.0;
for (int i = 0; i < 3; ++i)
{
PLib::EvalPolynomial(i, 0, 2, 1, aCoeff[0], aResult);
std::cout << "x=" << i << ", (2.0 + 2.0*x + 3.0*x^2): " << aResult << std::endl;
}
}
從上述代碼可以看出，OpenCASCADE的PLib包中對多項式的表示方法是采用的系數表示法。其系數分别為：2.0，2.0，3.0，即表示了幂級數：
當x=0, 1,2時計算結果如下圖所示：
Figure 2.1 Polynomial Evaluation
3.Polynomial Interpolation
多項式的插值問題是多項式求值的逆問題。多項式求值問題幾何意義就是已知曲線的表達式計算曲線上的點；而插值問題是已經曲線上的一些點來求通過這些點的曲線表達式。多項式插值中最常見最基本的問題是求一次數不超過n的代數多項式：
使
滿足插值條件的多項式稱為函數f(x)在節點xi處的n次插值多項式。由插值條件可知，插值多項式的系數滿足線性方程組：
由線性代數可知，其系數行列式是n+1階Vandermonde行列式，且：
因為插值的點是不同的點，是以行列式V不為0，即線性方程組有唯一解。這也是算法導論一書中這樣說“從某種意義上說，多項式的系數表示法與點值表示法是等價的，即用點值形式表示的多項式都對應唯一一個系數形式的多項式”的理論依據。即插值多項式的唯一性。這裡不僅指出了插值多項式存在的唯性，而且也提供了一種解法，即通過解線性方程組來确定系數。根據這個思路，給出對上面已知系數求值的多項式進行插值，代碼如下所示：
void testPolynomialInterpolation(void)
// given three points: (0, 2), (1, 7), (2, 18) to interpolate a polynomial
math_Matrix A(1, 3, 1, 3, 0.0);
math_Vector B(1, 3);
math_Vector X(1, 3);
A(1, 1) = 1.0; A(1, 2) = 0.0; A(1, 3) = 0.0; B(1) = 2.0;
A(2, 1) = 1.0; A(2, 2) = 1.0; A(2, 3) = 1.0; B(2) = 7.0;
A(3, 1) = 1.0; A(3, 2) = 2.0; A(3, 3) = 4.0; B(3) = 18.0;
// solve functions: Ax = B
math_Gauss aSolver(A);
aSolver.Solve(B, X);
if (aSolver.IsDone())
std::cout << X << std::endl;
已知多項式通過三個點（0, 2），（1, 7），（2, 18)，求通過這三個點的多項式表達式。根據插值條件列出線性方程組如下：
将系數a0,a1,a2看成線性方程組的待求變量，使用類math_Gauss來對線性方程組進行求解，計算結果如下所示：
Figure 3.1 Polynomial Interpolation Result
由上圖可知，對線性方程組的求解結果與上節點的系數對應。
4.Lagrange Interpolation
為什麼優雅浪漫的法國人的數學大師層出不窮呢？因為他們最優秀的人在學習數學。數學已經成為法國人傳統文化中最優秀的一部分了。
數學也是中國傳統文化的一部分，中國古代數學成就也很多：《周髀算經》； 《九章算術》（三國時劉徽著）； 祖沖之； 算盤。天文學：天象觀察記錄， 發明觀測儀器：圭表；渾儀；簡儀；高表；仰儀，制定曆法(農曆)。想想後來為什麼沒有大的發展，原因可能是科舉制度造成的，考試的内容偏文。
Lagrange父親是法國陸軍騎兵裡的一名軍官，後由于經商破産，家道中落。據Lagrange本人回憶，如果幼年時家境富裕，他也就不會作數學研究了。經曆了挫折之後沒被打倒的人後期成就會更大，像《紅樓夢》的作者曹雪芹。嚴重跑題了，回到Lagrange插值問題上來。從Lagrange中值定理的證明及Lagrange乘數法求極值的方式中可以看出Lagrange有個特點，那就是喜歡引入輔助函數來解決問題，可以看出Lagrange是非常精明的。對多項式插值也不例外，通過構造了一個Lagrange插值基函數來簡化多項式的插值，如下公式為Lagrange插值基函數：
則Lagrange插值基函數可以表示為簡潔的形式：
則n次多項式
滿足插值條件。OpenCASCADE的PLib包中也提供了Lagrange插值的函數來進行多項式插值計算。其用法的代碼如下所示：
void testLagrangeInterpolation(void)
// given three points: (0,2), (1,7), (2,18) to interpolate a polynomial
Standard_Real aValues[3] = {2.0, 7.0, 18.0};
Standard_Real aParameters[3] = {0.0, 1.0, 2.0};
// this do not output the coeff of the interpolate polynomial
PLib::EvalLagrange(1.5, 0, 2, 1, aValues[0], aParameters[0], aResult);
std::cout << "Result: " << aResult << std::endl;
下面對PLib::EvalLagrange()函數的7個參數進行說明：
Parameter
待根據Lagrange插值多項式求值的參數
DerivativeRequest
插值多項式導數次數
Degree
插值多項式的次數
Dimension
插值多項式的維次
Values
插值多項式的值
Parameters
插值多項式的參數
Results
參數Parameter在Lagrange插值多項式中的值
上述代碼計算的是這樣一個問題，已知f(0)=2; f(1)=7, f(2)=18，求f(1.5)的近似值。計算結果如下圖所示：
Figure 4.1 Lagrange Interpolation Result
// sqrt(100)=10, sqrt(121)=11, sqrt(144)=12, evaluate sqrt(115) value.
Standard_Real aSqrtValues[3] = {10.0, 11.0, 12.0};
Standard_Real aSqrtParameters[3] = {100.0, 121.0, 144.0};
// linear interpolation
PLib::EvalLagrange(115.0, 0, 1, 1, aSqrtValues[0], aSqrtParameters[0], aResult);
std::cout << "Linear Interpolate Result: " << aResult << std::endl;
// Parabolic Interpolation
PLib::EvalLagrange(115.0, 0, 2, 1, aSqrtValues[0], aSqrtParameters[0], aResult);
std::cout << "Parabolic Interpolate Result: " << aResult << std::endl;
計算結果如下所示：
Figure 4.2 Linear Interpolate and Parabolic Interpolate Result
将上圖4.2的結果與書中的計算結果進行對比發現，下面的結果為《計算方法》書中的結果：
PLib中Lagrange插值結果準确，精度較高。
從上面的結果來看，Lagrange插值方法比直接解線性方程組的方法要簡單，且Lagrange插值法比解線性方程組的實作要簡單很多，隻用一個函數即可。
5.Conclusion
通過将最簡單的多項式進行求值，及用求解線性方程組的方法插值和Lagrange插值多項式，來學習曲線拟合中要求較高“插值”，因為其要求曲線嚴格通過插值點。曲線拟合中的逼近就沒有這個要求，隻是要求曲線與插值點之間的容差盡量小。
通過應用OpenCASCADE的PLib包中的函數，可以發現對幂級數的多項式的表示法一般會用系數表示法，且都會使用高效的Horner法則，也是秦九韶算法。
直接根據定義來插值幂次多項式時，可以使用Gauss消元法來求解線性方程組。這種方式計算工作量大，而Lagrange插值法結構緊湊，便于程式設計實作，且代碼相對簡單。通過對《計算方法》書中例題的計算，來驗證PLib::EvalLagrange()函數的用法及計算結果。
6.References
1. 同濟大學數學教研室. 高等數學. 高等教育出版社. 1996
2. 同濟大學應用數學系. 線性代數. 高等教育出版社. 2003
3. Thomas H. Cormen. Introduction to Algorithms. The MIT Press. 2001
4. Les Piegl, Wayne Tiller. The NURBS Book. Springer-Verlag. 1995
5. Shing Liu. Polynomial Library in OpenCASCADE. 2013
<a href="http://www.cppblog.com/eryar/archive/2013/05/08/200118.html">http://www.cppblog.com/eryar/archive/2013/05/08/200118.html</a>
6. 易大義,沈雲寶,李有法. 計算方法. 浙江大學出版社. 2002
7. 蔣爾雄,趙風光,蘇仰鋒. 數值逼近. 複旦大學出版社. 2012
8. 王仁宏,李崇君,朱春鋼. 計算幾何教程. 科學出版社. 2008
10. 科學出版社名詞室. 新漢英數學詞彙. 科學出版社. 2004