OpenCASCADE Root-Finding Algorithm

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

代数方程求根难点是3个古老的数学标题,早在16世纪就找到了3遍、4次方程的求根公式。但停止19世纪才说明n>=八次的一般代数方程式无法用代数公式求解。在工程和科技(science and technology)中,许多难点平日归结为求解非线性方程的难题,由此,供给研讨用数值方法求得知足一定精度的代数方程式的近似解。 

作者国西魏南宋科学家秦九韶在她1247年所著的《数书天问》中,给出3个求代数方程根的好像方法,这些主意一般书上都称之为和纳Horner方法(United Kingdom化学家W.G.Horner)。实际上Horner在1819年才建议那么些艺术,比秦九韶晚五百多年。每当看到教科书中如此的牵线不知是该骄傲,依然该置之不顾。古人发明成立的事物比美国人早,而现行反革命境内用于CAD、CAM的软件大都是国外进口的,像CATIA,AutoCAD,Pro/E,UG
NX,SolidWorks,AVEVA
Plant/马林e,Intergraph,ACIS,Parasolid……等等一种类,很少见到中华软件的身形。而那一个软件广泛应用于飞行、造船、机械设计创制、工厂设计等相继行业,每年的软件授权花费不知几何?衷心希望当代国人奋发作为,为世界扩充色彩。 

闲话少说,本文首要关注非线性方程的数值解法,重点介绍了Newton-Rophson法及在OpenCASCADE中应用,即求点到曲线曲面包车型大巴极值,也即是曲线曲面点的反求参数难点。对数值算法感兴趣的读者,能够参考《数值分析》、《总计方式》之类的书籍以博取更详细新闻。 

2.Numerical Methods

方程求根的方式有成都百货上千,在《数学手册》中历数了如下一些主意: 

v 秦九韶法; 

v 二分法; 

v 迭代法; 

v 牛顿法Newton’s Method; 

v 弦截法; 

v 抛物线法; 

v 林士谔-赵访熊法; 

当中二分法是求实根的近似总结中央银卓有成效的最简易的格局,它只须要函数是接连的,因而它的运用限制很广,并有利于在微型总括机上落实,然则它无法求重根也不能求虚根,且流失较慢。 

Newton法在单根邻近收敛快,具有二阶收敛速度,但牛顿法对初值供给相比较苛刻,即必要初值接纳丰富靠近方程的根,不然Newton法大概不流失。扩充初值的选拔范围,可应用Newton下山法。 

Newton’s Method的落到实处原理的演示动画如下图所示: 

http://upload.wikimedia.org/wikipedia/commons/e/e0/NewtonIteration_Ani.gif

科学技术 1

Figure 2.1 Newton’s Method(Newton-Raphson)  

由地点的卡通片能够知晓了然牛顿法的原理。用数学的文字描述如下:设f(x)3回接二连三可导,xk是f(x)=0的第k次近似解。大家用曲线y=f(x)过点(xk,yk)的切线Lk: 

科学技术 2

来就像是曲线f(x)。取Lk与X轴的交点为f(x)=0的第k+一次类似解为: 

科学技术 3

科学技术 4

Figure 3.2 Newton-Raphson Method 

其中: 

科学技术 5

称为Newton公式。 

Newton法对初值x0供给苛刻,在实际应用中屡屡难以满意。Newton下山法是一种下降对初值必要的勘误Newton法。 

关于Newton方法的公开课的摄像笔者找到微博上有节课,介绍了Newton方法的法则及用法,网址为:http://v.163.com/movie/2006/8/T/V/M6GLI5A07_M6GLLGSTV.html,在后半部分。老师用实际例子来讲学还挺有趣的,感兴趣的读者也得以完全地探访,也可复习下微积分的知识点。 

3.OpenCASCADE Function Root

OpenCASCADE的math包中贯彻了方程求根的算法,相关的类有math_FunctionRoot,math_FunctionRoots,math_NewtonFunctionRoot等。在《Fundation
Classes User’s
Guide》中有对通用数学算法的牵线,即OpenCASCADE中达成了周边的数学算法: 

v 求解线性代数方程的根的算法; 

v 查找方程非常小值的算法; 

v 查找非线性方程(组)的根; 

v 查找对角矩阵的特征值及特征向量的算法; 

装有的数学算法以相同的格局来兑现,即:在构造函数中来做当先59%的乘除,从而给出适当的参数。全部相关数据都保存到结果对象中,因而具有的总括尽量以最飞快的点子来求解。函数IsDone()评释计算成功。如下所示分别为使用分裂的算法来计算如下方程在[0,2]距离上的根: 

科学技术 6

兑现程序代码如下所示: 

/*
*    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派生且重载其多个纯虚函数Value(),
Derivative(),
Values(),在那四个纯虚函数中计算有关的值及导数值即可。所以实际上利用时,正确重载那四个函数是毋庸置疑利用求根算法的最首要。 

求根用了几个不等的类,即二种格局来达成: 

v math_FunctionRoot:即Newton-Raphson法; 

v math_BissecNewton:是Newton-Raphson和二分法的结缘算法; 

v math_NewtonFunctionRoot:Newton Method; 

总结结果如下图所示: 

科学技术 7

Figure 3.1 Root-Finding result of OpenCASCADE 

由总结结果可见,两种格局计算的结果一律,都以1.13472,与书中结果符合。不过math_牛顿FunctionRoot的迭代次比math_FunctionRoot多二遍,且计量精度要低很多。 

使用math_BissecNewton求根不用安装早先值,相比较便宜,精度与math_FunctionRoot一致。 

科学技术,4.Application

在工程和科学和技术中,许多题目平时总结为求解非线性方程的难题。在OpenCASCADE中的应用越来越多了,从上边一张类图知秋一叶: 

科学技术 8

Figure 4.1 math_FunctionWithDerivative class diagram 

由图可以,从类math_FunctionWithDerivative派生出了成都百货上千可导函数的类,那么些函数都可用于求根的类中,从而计算出相应方程的根。 

上边给出3个实际采纳,即曲线曲面上点的反求难题,来证实什么选择上述求根算法来消除实际的难点。由于曲线曲面包车型地铁参数表示法,通过参数u或(u,v)能够方便地求出曲线上的点或曲面上的点。若给定一个点P(x,y,z),假使它在p次NURBS曲线C(u)上,求对应的参数u’使得C(u’)=P,这些题材称为点的反求(point
inverse)。在OpenCASCADE中提供了简约的函数接口来兑现点的反求,使用类为吉优mLib_Tool: 

科学技术 9

如何将点的反求难题归纳为方程求根的标题,就要依据规则来树立方程了。一种简单并完全能够缓解这一题材的点子是:利用牛顿迭代法来最小化点P和C(u)的相距,如下图所示。借使最小距离小于一个事先钦赐的精度值,则认为点P在曲线上。那种格局有效地解决了更相像的“点在曲线上的黑影”的难题。 

因为方程求根的Newton方法供给钦定初值u0,所以可按如下方法获得一个用以Newton法的初值u0: 

v
假设已知点P在加以精度内放在曲线上,则用强凸包性鲜明候选的段,对于一般的点到曲线的黑影难点,则选用所要的段作为候选段; 

v
在每一种候选段上,计算n个按参数等间隔分布的点。总计出富有这几个点和点P的偏离,选拔在那之中距点P近期的点的参数作为u0。点数n一般选拔某种启发的章程来选拔。 

科学技术 10

Figure 4.2 Point projection and Point inversion 

需求强调的是利用Newton方法,3个好的初值对于迭代的收敛性及没有速度是特别重庆大学的。未来只要已经分明了初值u0,利用数据积定义函数: 

科学技术 11

任凭P点是不是位于曲线上,当f(u)=0时,点P到C(u)的相距达到最小。对f(u)求导得: 

科学技术 12

代入Newton迭代公式得: 

科学技术 13

在OpenCASCADE中曲线点的反求重倘若应用了派生自math_FunctionWithDerivative的类Extrema_PCFOfEPCOfExtPC,类图如下所示: 

科学技术 14

Figure 4.3 class diagram for point inverstion 

据此需求贯彻八个纯虚函数Value(), Derivative(),
Values(),将其落到实处代码列出如下所示: 

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;
}

听大人讲代码可知,实现原理与上述同样。下边给出三个简便的事例,来证实及便利调节和测试点的反求算法。示例程序代码如下所示: 

/*
*    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;
}

依据圆上一点,求出对应的参数值,总计结果如下所示: 

科学技术 15

5.Conclusion

牛顿法能够选作对导数能有效求值,且导数在根的邻域中年老年是的别的函数方程的求根方法。Newton法在单根邻近收敛快,精度高,具有二阶收敛速度,但Newton法对初值须要相比较高,即须求初值采取丰硕靠近方程的根,否则Newton法只怕不收敛。 

OpenCASCADE的math包中提供了求根的两种达成算法,就算代码有些乱,然则那种肤浅的商量照旧非凡不错的,便于扩充应用。精通了math_FunctionRoot的算法,进而能够理解从math_FunctionWithDerivative派生的类的法则了。 

透过曲线上点的反求难点引出使用求根算法的现实性实例,从中能够见到关键照旧要将实际难点抽象成1个方程。有了方程,依据Newton迭代公式,求出相应的值和导数值,就能够收获方程的高精度的根了。 

对数值算法感兴趣的读者,能够参见《总计方法》、《数值分析》等有关书籍,从而得以在知情OpenCASCADE的代码的底子上,能够团结来兑现相关算法。 

6. References

  1. 数学手册编写组. 数学手册. 高教出版社. 一九七六 

  2. 赵罡,穆国旺,王拉柱译. 非均匀有理B样条. 复旦东军政学院学出版社. 二零零六 

  3. Les Piegl, Wayne Tiller. The NURBS Book. Springer-Verlag. 1997 

  4. 易大义,沈云宝,李有法编. 计算方法. 吉林高校出版社. 二零零三 

  5. 易大义,陈道琦编. 数值分析引论. 海南大学出版社. 壹玖玖玖 

  6. 李庆杨,王能超,易大义.数值分析.华中理哲高校出版社. 1990 

  7. 同济数学教学讨论室. 高等数学(第六版). 高教出版社. 一九九六 

  8. Newton’s Method video:  

http://v.163.com/movie/2006/8/T/V/M6GLI5A07_M6GLLGSTV.html

  1. http://en.wikipedia.org/wiki/Root-finding_algorithm

  2. http://mathworld.wolfram.com/Root-FindingAlgorithm.html

  3. http://mathworld.wolfram.com/NewtonsMethod.html

 

PDF Version: OpenCASCADE Root-Finding
Algorithm

Leave a Comment.