OpenCASCADE Interpolation - Lagrange
eryar@163.com
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，原來他也是OpenCASCADE的發源地的人：法國人。Joseph-Louis Lagrange，法國著名數學家、物理學家。1736年1月25日生于意大利都靈，1813年4月10日卒于法國巴黎。他在數學、力學和天文學三個學科領域都有歷史性的貢獻，其中尤以數學方面的成就最為突出。以下圖片來自網易公開課《數學傳奇》：從笛卡爾到龐加萊—法國數學的人文傳統，公開課網址：http://open.163.com/special/cuvocw/shuxuechuanqi.html
為什么優雅浪漫的法國人的數學大師層出不窮呢？因為他們最優秀的人在學習數學。數學已經成為法國人傳統文化中最優秀的一部分了。
數學也是中國傳統文化的一部分，中國古代數學成就也很多：《周髀算經》； 《九章算術》（三國時劉徽著）； 祖沖之； 算盤。天文學：天象觀察記錄， 發明觀測儀器：圭表；渾儀；簡儀；高表；仰儀，制定歷法(農歷)。想想后來為什么沒有大的發展，原因可能是科舉制度造成的，考試的內容偏文。
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};
Standard_Real aResult = 0.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
由上圖可知，經過Lagrange拋物插值得到的結果與直接計算得到的結果吻合。為了更好的說明OpenCASCADE中PLib包的Lagrange插值的用意，下面將《計算方法》[6]中Lagrange插值部分中的例題進行計算來結計算結果進行比較。已知
分別用線性插值和拋物插值求
的近似值。相關計算代碼如下所示：
// 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
http://m.shnenglu.com/eryar/archive/2013/05/08/200118.html
6. 易大義,沈云寶,李有法. 計算方法. 浙江大學出版社. 2002
7. 蔣爾雄,趙風光,蘇仰鋒. 數值逼近. 復旦大學出版社. 2012
8. 王仁宏,李崇君,朱春鋼. 計算幾何教程. 科學出版社. 2008
9. 蔡天新. 數學傳奇. http://open.163.com/special/cuvocw/shuxuechuanqi.html
10. 科學出版社名詞室. 新漢英數學詞匯. 科學出版社. 2004