OpenCASCADE Root-Finding Algorithm
eryar@163.com
Abstract. A root-finding algorithm is a numerical method, or algorithm, for finding a value x such that f(x)=0, for a given function f. Such an x is called a root of the function f. In OpenCASCADE math package, implemente Newton-Raphson method to find the root for a function. The algorithm can be used for finding the extrema value for curve or surface, .i.e Point Inversion, find the parameter for a point on the curve or surface. The paper focus on the usage of OpenCASCADE method and its application.
Key Words. OpenCASCADE, Extrema, Newton-Raphson, Root-Finding, FunctionRoot
1. Introduction
代數(shù)方程求根問題是一個(gè)古老的數(shù)學(xué)問題，早在16世紀(jì)就找到了三次、四次方程的求根公式。但直到19世紀(jì)才證明n>=5次的一般代數(shù)方程式不能用代數(shù)公式求解。在工程和科學(xué)技術(shù)中，許多問題常常歸結(jié)為求解非線性方程的問題，因此，需要研究用數(shù)值方法求得滿足一定精度的代數(shù)方程式的近似解。
我國(guó)古代宋朝數(shù)學(xué)家秦九韶在他1247年所著的《數(shù)書九章》中，給出一個(gè)求代數(shù)方程根的近似方法，這個(gè)方法一般書上都稱為和納Horner方法（英國(guó)數(shù)學(xué)家W.G.Horner）。實(shí)際上Horner在1819年才提出這個(gè)方法，比秦九韶晚五百多年。每當(dāng)看到教科書中這樣的介紹不知是該驕傲，還是該嗤之以鼻。古人發(fā)明創(chuàng)造的東西比外國(guó)人早，而現(xiàn)在國(guó)內(nèi)用于CAD、CAM的軟件大都是國(guó)外進(jìn)口的，像CATIA，AutoCAD，Pro/E，UG NX，SolidWorks，AVEVA Plant/Marine，Intergraph，ACIS，Parasolid……等等不勝枚舉，很少看到中國(guó)軟件的身影。而這些軟件廣泛應(yīng)用于航空、造船、機(jī)械設(shè)計(jì)制造、工廠設(shè)計(jì)等各個(gè)行業(yè)，每年的軟件授權(quán)費(fèi)用不知幾何？衷心希望當(dāng)代國(guó)人奮發(fā)作為，為世界增添色彩。
閑話少說，本文主要關(guān)注非線性方程的數(shù)值解法，重點(diǎn)介紹了Newton-Rophson法及在OpenCASCADE中應(yīng)用，即求點(diǎn)到曲線曲面的極值，也就是曲線曲面點(diǎn)的反求參數(shù)問題。對(duì)數(shù)值算法感興趣的讀者，可以參考《數(shù)值分析》、《計(jì)算方法》之類的書籍以獲取更詳細(xì)信息。
2.Numerical Methods
方程求根的方法有很多，在《數(shù)學(xué)手冊(cè)》中列舉了如下一些方法：
v 秦九韶法；
v 二分法；
v 迭代法；
v 牛頓法Newton’s Method；
v 弦截法；
v 拋物線法；
v 林士諤－趙訪熊法；
其中二分法是求實(shí)根的近似計(jì)算中行之有效的最簡(jiǎn)單的方法，它只要求函數(shù)是連續(xù)的，因此它的使用范圍很廣，并便于在計(jì)算機(jī)上實(shí)現(xiàn)，但是它不能求重根也不能求虛根，且收斂較慢。
Newton法在單根鄰近收斂快，具有二階收斂速度，但Newton法對(duì)初值要求比較苛刻，即要求初值選取充分靠近方程的根，否則Newton法可能不收斂。擴(kuò)大初值的選取范圍，可采用Newton下山法。
Newton’s Method的實(shí)現(xiàn)原理的演示動(dòng)畫如下圖所示：
http://upload.wikimedia.org/wikipedia/commons/e/e0/NewtonIteration_Ani.gif
Figure 2.1 Newton’s Method(Newton-Raphson)
由上面的動(dòng)畫可以清楚理解Newton法的原理。用數(shù)學(xué)的文字描述如下：設(shè)f(x)二次連續(xù)可導(dǎo)，xk是f(x)=0的第k次近似解。我們用曲線y=f(x)過點(diǎn)(xk,yk)的切線Lk：
來近似曲線f(x)。取Lk與X軸的交點(diǎn)為f(x)=0的第k+1次近似解為：
Figure 3.2 Newton-Raphson Method
其中：
稱為Newton公式。
Newton法對(duì)初值x0要求苛刻，在實(shí)際應(yīng)用中往往難以滿足。Newton下山法是一種降低對(duì)初值要求的修正Newton法。
關(guān)于Newton方法的公開課的視頻我找到網(wǎng)易上有節(jié)課，介紹了Newton方法的原理及用法，網(wǎng)址為：http://v.163.com/movie/2006/8/T/V/M6GLI5A07_M6GLLGSTV.html，在后半部分。老師用實(shí)際例子來講解還挺有意思的，感興趣的讀者也可以完整地看看，也可復(fù)習(xí)下微積分的知識(shí)點(diǎn)。
3.OpenCASCADE Function Root
OpenCASCADE的math包中實(shí)現(xiàn)了方程求根的算法，相關(guān)的類有math_FunctionRoot，math_FunctionRoots，math_NewtonFunctionRoot等。在《Fundation Classes User’s Guide》中有對(duì)通用數(shù)學(xué)算法的介紹，即OpenCASCADE中實(shí)現(xiàn)了常見的數(shù)學(xué)算法：
v 求解線性代數(shù)方程的根的算法；
v 查找方程極小值的算法；
v 查找非線性方程（組）的根；
v 查找對(duì)角矩陣的特征值及特征向量的算法；
所有的數(shù)學(xué)算法以相同的方式來實(shí)現(xiàn)，即：在構(gòu)造函數(shù)中來做大部分的計(jì)算，從而給出適當(dāng)?shù)膮?shù)。所有相關(guān)數(shù)據(jù)都保存到結(jié)果對(duì)象中，因此所有的計(jì)算盡量以最高效的方式來求解。函數(shù)IsDone()表明計(jì)算成功。如下所示分別為采用不同的算法來計(jì)算如下方程在[0，2]區(qū)間上的根：
實(shí)現(xiàn)程序代碼如下所示：
/*
* Copyright (c) 2014 eryar All Rights Reserved.
*
* File : Main.cpp
* Author : eryar@163.com
* Date : 2014-10-20 18:52
* Version : 1.0v
*
* Description : Test OpenCASCADE function root algorithm.
*
* Key words : OpenCASCADE, Newton-Raphson, Root-Finding Algorithm, FunctionRoot
*/
#define WNT
#include <Precision.hxx>
#include <math_FunctionWithDerivative.hxx>
#include <math_BissecNewton.hxx>
#include <math_FunctionRoot.hxx>
#include <math_NewtonFunctionRoot.hxx>
#pragma comment(lib, "TKernel.lib")
#pragma comment(lib, "TKMath.lib")
class TestFunction : public math_FunctionWithDerivative
{
public:
virtual Standard_Boolean Value(const Standard_Real X, Standard_Real& F)
{
F = pow(X, 6) - X - 1;
return Standard_True;
}
virtual Standard_Boolean Derivative(const Standard_Real X, Standard_Real& D)
{
D = 6 * pow(X, 5) - 1;
return Standard_True;
}
virtual Standard_Boolean Values(const Standard_Real X, Standard_Real& F, Standard_Real& D)
{
Value(X, F);
Derivative(X, D);
return Standard_True;
}
};
void TestFunctionRoot(void)
{
TestFunction aFunction;
math_FunctionRoot aSolver1(aFunction, 1.5, 0.0, 0.0, 2.0);
math_BissecNewton aSolver2(aFunction, 0.0, 2.0, 0.0);
math_NewtonFunctionRoot aSolver3(aFunction, 1.5, Precision::Confusion(), Precision::Confusion());
std::cout << aSolver1 << std::endl;
std::cout << aSolver2 << std::endl;
std::cout << aSolver3 << std::endl;
}
int main(int argc, char* argv[])
{
TestFunctionRoot();
return 0;
}
由上述代碼可知，要想使用求根算法，必須從math_FunctionWithDerivative派生且重載其三個(gè)純虛函數(shù)Value(), Derivative(), Values()，在這三個(gè)純虛函數(shù)中計(jì)算相關(guān)的值及導(dǎo)數(shù)值即可。所以實(shí)際使用時(shí)，正確重載這三個(gè)函數(shù)是正確使用求根算法的關(guān)鍵。
求根用了三個(gè)不同的類，即三種方法來實(shí)現(xiàn)：
v math_FunctionRoot：即Newton-Raphson法；
v math_BissecNewton：是Newton-Raphson和二分法的組合算法；
v math_NewtonFunctionRoot：Newton Method；
計(jì)算結(jié)果如下圖所示：
Figure 3.1 Root-Finding result of OpenCASCADE
由計(jì)算結(jié)果可知，三種方法計(jì)算的結(jié)果相同，都是1.13472，與書中結(jié)果吻合。但是math_NewtonFunctionRoot的迭代次比math_FunctionRoot多一次，且計(jì)算精度要低很多。
使用math_BissecNewton求根不用設(shè)置初始值，比較方便，精度與math_FunctionRoot一致。
4.Application
在工程和科學(xué)技術(shù)中，許多問題常常歸結(jié)為求解非線性方程的問題。在OpenCASCADE中的應(yīng)用更多了，從下面一張類圖可見一斑：
Figure 4.1 math_FunctionWithDerivative class diagram
由圖可知，從類math_FunctionWithDerivative派生出了很多可導(dǎo)函數(shù)的類，這些函數(shù)都可用于求根的類中，從而計(jì)算出對(duì)應(yīng)方程的根。
下面給出一個(gè)實(shí)際應(yīng)用，即曲線曲面上點(diǎn)的反求問題，來說明如何應(yīng)用上述求根算法來解決實(shí)際的問題。由于曲線曲面的參數(shù)表示法，通過參數(shù)u或(u,v)可以方便地求出曲線上的點(diǎn)或曲面上的點(diǎn)。若給定一個(gè)點(diǎn)P(x,y,z)，假設(shè)它在p次NURBS曲線C(u)上，求對(duì)應(yīng)的參數(shù)u’使得C(u’)=P，這個(gè)問題稱為點(diǎn)的反求（point inverse）。在OpenCASCADE中提供了簡(jiǎn)單的函數(shù)接口來實(shí)現(xiàn)點(diǎn)的反求，使用類為GeomLib_Tool：
如何將點(diǎn)的反求問題歸結(jié)為方程求根的問題，就要根據(jù)條件來建立方程了。一種簡(jiǎn)單并完全可以解決這一問題的方法是：利用Newton迭代法來最小化點(diǎn)P和C(u)的距離，如下圖所示。如果最小距離小于一個(gè)事先指定的精度值，則認(rèn)為點(diǎn)P在曲線上。這種方法有效地解決了更一般的“點(diǎn)在曲線上的投影”的問題。
因?yàn)榉匠糖蟾腘ewton方法需要指定初值u0，所以可按如下方法得到一個(gè)用于Newton法的初值u0：
v 如果已知點(diǎn)P在給定精度內(nèi)位于曲線上，則用強(qiáng)凸包性確定候選的段，對(duì)于一般的點(diǎn)到曲線的投影問題，則選擇所要的段作為候選段；
v 在每個(gè)候選段上，計(jì)算n個(gè)按參數(shù)等間隔分布的點(diǎn)。計(jì)算出所有這些點(diǎn)和點(diǎn)P的距離，選擇其中距點(diǎn)P最近的點(diǎn)的參數(shù)作為u0。點(diǎn)數(shù)n一般利用某種啟發(fā)的方法來選擇。
Figure 4.2 Point projection and Point inversion
需要強(qiáng)調(diào)的是使用Newton方法，一個(gè)好的初值對(duì)于迭代的收斂性及收斂速度是非常重要的。現(xiàn)在假設(shè)已經(jīng)確定了初值u0，利用數(shù)量積定義函數(shù)：
不管P點(diǎn)是否位于曲線上，當(dāng)f(u)=0時(shí)，點(diǎn)P到C(u)的距離達(dá)到最小。對(duì)f(u)求導(dǎo)得：
代入Newton迭代公式得：
在OpenCASCADE中曲線點(diǎn)的反求主要是使用了派生自math_FunctionWithDerivative的類Extrema_PCFOfEPCOfExtPC，類圖如下所示：
Figure 4.3 class diagram for point inverstion
所以需要實(shí)現(xiàn)三個(gè)純虛函數(shù)Value(), Derivative(), Values()，將其實(shí)現(xiàn)代碼列出如下所示：
Standard_Boolean Extrema_FuncExtPC::Value (const Standard_Real U, Standard_Real& F)
{
if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
myU = U;
Vec D1c;
Tool::D1(*((Curve*)myC),myU,myPc,D1c);
Standard_Real Ndu = D1c.Magnitude();
if (Ndu <= Tol) { // Cas Singulier (PMN 22/04/1998)
Pnt P1, P2;
P2 = Tool::Value(*((Curve*)myC),myU + delta);
P1 = Tool::Value(*((Curve*)myC),myU - delta);
Vec V(P1,P2);
D1c = V;
Ndu = D1c.Magnitude();
if (Ndu <= Tol) {
Vec aD2;
Tool::D2(*((Curve*)myC),myU,myPc,D1c,aD2);
Ndu = aD2.Magnitude();
if(Ndu <= Tol)
return Standard_False;
D1c = aD2;
}
}
Vec PPc (myP,myPc);
F = PPc.Dot(D1c)/Ndu;
return Standard_True;
}
//=============================================================================
Standard_Boolean Extrema_FuncExtPC::Derivative (const Standard_Real U, Standard_Real& D1f)
{
if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
Standard_Real F;
return Values(U,F,D1f); /* on fait appel a Values pour simplifier la
sauvegarde de l'etat. */
}
//=============================================================================
Standard_Boolean Extrema_FuncExtPC::Values (const Standard_Real U, Standard_Real& F, Standard_Real& D1f)
{
if (!myPinit || !myCinit) Standard_TypeMismatch::Raise();
myU = U;
Vec D1c,D2c;
Tool::D2(*((Curve*)myC),myU,myPc,D1c,D2c);
Standard_Real Ndu = D1c.Magnitude();
if (Ndu <= Tol) {// Cas Singulier (PMN 22/04/1998)
Pnt P1, P2;
Vec V1;
Tool::D1(*((Curve*)myC),myU+delta, P2, V1);
Tool::D1(*((Curve*)myC),myU-delta, P1, D2c);
Vec V(P1,P2);
D1c = V;
D2c -= V1;
Ndu = D1c.Magnitude();
if (Ndu <= Tol) {
myD1Init = Standard_False;
return Standard_False;
}
}
Vec PPc (myP,myPc);
F = PPc.Dot(D1c)/Ndu;
D1f = Ndu + (PPc.Dot(D2c)/Ndu) - F*(D1c.Dot(D2c))/(Ndu*Ndu);
myD1f = D1f;
myD1Init = Standard_True;
return Standard_True;
}
根據(jù)代碼可知，實(shí)現(xiàn)原理與上述一致。下面給出一個(gè)簡(jiǎn)單的例子，來說明及方便調(diào)試點(diǎn)的反求算法。示例程序代碼如下所示：
/*
* Copyright (c) 2014 eryar All Rights Reserved.
*
* File : Main.cpp
* Author : eryar@163.com
* Date : 2014-10-20 18:52
* Version : 1.0v
*
* Description : Test OpenCASCADE function root algorithm.
*
* Key words : OpenCascade, Extrema, Newton's Method
*/
#define WNT
#include <math_FunctionRoots.hxx>
#include <math_NewtonFunctionRoot.hxx>
#include <Extrema_PCFOfEPCOfExtPC.hxx>
#include <GC_MakeCircle.hxx>
#include <GeomAdaptor_Curve.hxx>
#pragma comment(lib, "TKernel.lib")
#pragma comment(lib, "TKMath.lib")
#pragma comment(lib, "TKG3d.lib")
#pragma comment(lib, "TKGeomBase.lib")
void TestExtrema(void)
{
Handle_Geom_Curve aCircle = GC_MakeCircle(gp::XOY(), 2.0);
GeomAdaptor_Curve aCurve(aCircle);
Extrema_PCFOfEPCOfExtPC aFunction(aCircle->Value(0.123456789), aCurve);
math_FunctionRoots aSolver1(aFunction, -2.0, 2.0, 10);
math_NewtonFunctionRoot aSolver2(aFunction, 0.0, 0.0, 0.0);
aSolver1.Dump(std::cout);
std::cout << "========================================" << std::endl;
aSolver2.Dump(std::cout);
}
int main(int argc, char* argv[])
{
TestExtrema();
return 0;
}
根據(jù)圓上一點(diǎn)，求出對(duì)應(yīng)的參數(shù)值，計(jì)算結(jié)果如下所示：
5.Conclusion
Newton法可以選作對(duì)導(dǎo)數(shù)能有效求值，且導(dǎo)數(shù)在根的鄰域中連續(xù)的任何函數(shù)方程的求根方法。Newton法在單根鄰近收斂快，精度高，具有二階收斂速度，但Newton法對(duì)初值要求比較高，即要求初值選取充分靠近方程的根，否則Newton法可能不收斂。
OpenCASCADE的math包中提供了求根的幾種實(shí)現(xiàn)算法，雖然代碼有些亂，但是這種抽象的思想還是相當(dāng)不錯(cuò)的，便于擴(kuò)展應(yīng)用。理解了math_FunctionRoot的算法，進(jìn)而可以理解從math_FunctionWithDerivative派生的類的原理了。
通過曲線上點(diǎn)的反求問題引出使用求根算法的具體實(shí)例，從中可以看出關(guān)鍵還是要將實(shí)際問題抽象成一個(gè)方程。有了方程，根據(jù)Newton迭代公式，求出相應(yīng)的值和導(dǎo)數(shù)值，就可以得到方程的高精度的根了。
對(duì)數(shù)值算法感興趣的讀者，可以參考《計(jì)算方法》、《數(shù)值分析》等相關(guān)書籍，從而可以在理解OpenCASCADE的代碼的基礎(chǔ)上，可以自己來實(shí)現(xiàn)相關(guān)算法。
6. References
1. 數(shù)學(xué)手冊(cè)編寫組. 數(shù)學(xué)手冊(cè). 高等教育出版社. 1979
2. 趙罡，穆國(guó)旺，王拉柱譯. 非均勻有理B樣條. 清華大學(xué)出版社. 2010
3. Les Piegl, Wayne Tiller. The NURBS Book. Springer-Verlag. 1997
4. 易大義，沈云寶，李有法編. 計(jì)算方法. 浙江大學(xué)出版社. 2002
5. 易大義，陳道琦編. 數(shù)值分析引論. 浙江大學(xué)出版社. 1998
6. 李慶楊，王能超，易大義.數(shù)值分析.華中理工大學(xué)出版社. 1986
7. 同濟(jì)大學(xué)數(shù)學(xué)教研室. 高等數(shù)學(xué)(第四版). 高等教育出版社. 1996
8. Newton's Method video:
http://v.163.com/movie/2006/8/T/V/M6GLI5A07_M6GLLGSTV.html
9. http://en.wikipedia.org/wiki/Root-finding_algorithm
10. http://mathworld.wolfram.com/Root-FindingAlgorithm.html
11. http://mathworld.wolfram.com/NewtonsMethod.html