OpenCASCADE Root-Finding Algorithm

时间:2023-01-12 09:13:14

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

代数方程求根问题是一个古老的数学问题,早在16世纪就找到了三次、四次方程的求根公式。但直到19世纪才证明n>=5次的一般代数方程式不能用代数公式求解。在工程和科学技术中,许多问题常常归结为求解非线性方程的问题,因此,需要研究用数值方法求得满足一定精度的代数方程式的近似解。

我国古代宋朝数学家秦九韶在他1247年所著的《数书九章》中,给出一个求代数方程根的近似方法,这个方法一般书上都称为和纳Horner方法(英国数学家W.G.Horner)。实际上Horner在1819年才提出这个方法,比秦九韶晚五百多年。每当看到教科书中这样的介绍不知是该骄傲,还是该嗤之以鼻。古人发明创造的东西比外国人早,而现在国内用于CAD、CAM的软件大都是国外进口的,像CATIA,AutoCAD,Pro/E,UG NX,SolidWorks,AVEVA Plant/Marine,Intergraph,ACIS,Parasolid……等等不胜枚举,很少看到中国软件的身影。而这些软件广泛应用于航空、造船、机械设计制造、工厂设计等各个行业,每年的软件授权费用不知几何?衷心希望当代国人奋发作为,为世界增添色彩。

闲话少说,本文主要关注非线性方程的数值解法,重点介绍了Newton-Rophson法及在OpenCASCADE中应用,即求点到曲线曲面的极值,也就是曲线曲面点的反求参数问题。对数值算法感兴趣的读者,可以参考《数值分析》、《计算方法》之类的书籍以获取更详细信息。

2.Numerical Methods

方程求根的方法有很多,在《数学手册》中列举了如下一些方法:

v 秦九韶法;

v 二分法;

v 迭代法;

v 牛顿法Newton’s Method;

v 弦截法;

v 抛物线法;

v 林士谔-赵访熊法;

其中二分法是求实根的近似计算中行之有效的最简单的方法,它只要求函数是连续的,因此它的使用范围很广,并便于在计算机上实现,但是它不能求重根也不能求虚根,且收敛较慢。

Newton法在单根邻近收敛快,具有二阶收敛速度,但Newton法对初值要求比较苛刻,即要求初值选取充分靠近方程的根,否则Newton法可能不收敛。扩大初值的选取范围,可采用Newton下山法。

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

//bbsmax.ikafan.com/static/L3Byb3h5L2h0dHAvdXBsb2FkLndpa2ltZWRpYS5vcmcvd2lraXBlZGlhL2NvbW1vbnMvZS9lMC9OZXd0b25JdGVyYXRpb25fQW5pLmdpZg==.jpg

OpenCASCADE Root-Finding Algorithm

Figure 2.1 Newton’s Method(Newton-Raphson)

由上面的动画可以清楚理解Newton法的原理。用数学的文字描述如下:设f(x)二次连续可导,xk是f(x)=0的第k次近似解。我们用曲线y=f(x)过点(xk,yk)的切线Lk:

OpenCASCADE Root-Finding Algorithm

来近似曲线f(x)。取Lk与X轴的交点为f(x)=0的第k+1次近似解为:

OpenCASCADE Root-Finding Algorithm

OpenCASCADE Root-Finding Algorithm

Figure 3.2 Newton-Raphson Method

其中:

OpenCASCADE Root-Finding Algorithm

称为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 查找对角矩阵的特征值及特征向量的算法;

所有的数学算法以相同的方式来实现,即:在构造函数中来做大部分的计算,从而给出适当的参数。所有相关数据都保存到结果对象中,因此所有的计算尽量以最高效的方式来求解。函数IsDone()表明计算成功。如下所示分别为采用不同的算法来计算如下方程在[0,2]区间上的根:

OpenCASCADE Root-Finding Algorithm

实现程序代码如下所示:

/*
* 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, ) - X - ; return Standard_True;
} virtual Standard_Boolean Derivative(const Standard_Real X, Standard_Real& D)
{
D = * pow(X, ) - ; 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 ;
}

由上述代码可知,要想使用求根算法,必须从math_FunctionWithDerivative派生且重载其三个纯虚函数Value(), Derivative(), Values(),在这三个纯虚函数中计算相关的值及导数值即可。所以实际使用时,正确重载这三个函数是正确使用求根算法的关键。

求根用了三个不同的类,即三种方法来实现:

v math_FunctionRoot:即Newton-Raphson法;

v math_BissecNewton:是Newton-Raphson和二分法的组合算法;

v math_NewtonFunctionRoot:Newton Method;

计算结果如下图所示:

OpenCASCADE Root-Finding Algorithm

Figure 3.1 Root-Finding result of OpenCASCADE

由计算结果可知,三种方法计算的结果相同,都是1.13472,与书中结果吻合。但是math_NewtonFunctionRoot的迭代次比math_FunctionRoot多一次,且计算精度要低很多。

使用math_BissecNewton求根不用设置初始值,比较方便,精度与math_FunctionRoot一致。

4.Application

在工程和科学技术中,许多问题常常归结为求解非线性方程的问题。在OpenCASCADE中的应用更多了,从下面一张类图可见一斑:

OpenCASCADE Root-Finding Algorithm

Figure 4.1 math_FunctionWithDerivative class diagram

由图可知,从类math_FunctionWithDerivative派生出了很多可导函数的类,这些函数都可用于求根的类中,从而计算出对应方程的根。

下面给出一个实际应用,即曲线曲面上点的反求问题,来说明如何应用上述求根算法来解决实际的问题。由于曲线曲面的参数表示法,通过参数u或(u,v)可以方便地求出曲线上的点或曲面上的点。若给定一个点P(x,y,z),假设它在p次NURBS曲线C(u)上,求对应的参数u’使得C(u’)=P,这个问题称为点的反求(point inverse)。在OpenCASCADE中提供了简单的函数接口来实现点的反求,使用类为GeomLib_Tool:

OpenCASCADE Root-Finding Algorithm

如何将点的反求问题归结为方程求根的问题,就要根据条件来建立方程了。一种简单并完全可以解决这一问题的方法是:利用Newton迭代法来最小化点P和C(u)的距离,如下图所示。如果最小距离小于一个事先指定的精度值,则认为点P在曲线上。这种方法有效地解决了更一般的“点在曲线上的投影”的问题。

因为方程求根的Newton方法需要指定初值u0,所以可按如下方法得到一个用于Newton法的初值u0:

v 如果已知点P在给定精度内位于曲线上,则用强凸包性确定候选的段,对于一般的点到曲线的投影问题,则选择所要的段作为候选段;

v 在每个候选段上,计算n个按参数等间隔分布的点。计算出所有这些点和点P的距离,选择其中距点P最近的点的参数作为u0。点数n一般利用某种启发的方法来选择。

OpenCASCADE Root-Finding Algorithm

Figure 4.2 Point projection and Point inversion

需要强调的是使用Newton方法,一个好的初值对于迭代的收敛性及收敛速度是非常重要的。现在假设已经确定了初值u0,利用数量积定义函数:

OpenCASCADE Root-Finding Algorithm

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

OpenCASCADE Root-Finding Algorithm

代入Newton迭代公式得:

OpenCASCADE Root-Finding Algorithm

在OpenCASCADE中曲线点的反求主要是使用了派生自math_FunctionWithDerivative的类Extrema_PCFOfEPCOfExtPC,类图如下所示:

OpenCASCADE Root-Finding Algorithm

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

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

OpenCASCADE Root-Finding Algorithm

5.Conclusion

Newton法可以选作对导数能有效求值,且导数在根的邻域中连续的任何函数方程的求根方法。Newton法在单根邻近收敛快,精度高,具有二阶收敛速度,但Newton法对初值要求比较高,即要求初值选取充分靠近方程的根,否则Newton法可能不收敛。

OpenCASCADE的math包中提供了求根的几种实现算法,虽然代码有些乱,但是这种抽象的思想还是相当不错的,便于扩展应用。理解了math_FunctionRoot的算法,进而可以理解从math_FunctionWithDerivative派生的类的原理了。

通过曲线上点的反求问题引出使用求根算法的具体实例,从中可以看出关键还是要将实际问题抽象成一个方程。有了方程,根据Newton迭代公式,求出相应的值和导数值,就可以得到方程的高精度的根了。

对数值算法感兴趣的读者,可以参考《计算方法》、《数值分析》等相关书籍,从而可以在理解OpenCASCADE的代码的基础上,可以自己来实现相关算法。

6. References

1. 数学手册编写组. 数学手册. 高等教育出版社. 1979

2. 赵罡,穆国旺,王拉柱译. 非均匀有理B样条. 清华大学出版社. 2010

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

4. 易大义,沈云宝,李有法编. 计算方法. 浙江大学出版社. 2002

5. 易大义,陈道琦编. 数值分析引论. 浙江大学出版社. 1998

6. 李庆杨,王能超,易大义.数值分析.华中理工大学出版社. 1986

7. 同济大学数学教研室. 高等数学(第四版). 高等教育出版社. 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

OpenCASCADE Root-Finding Algorithm的更多相关文章

  1. HDU4889 Scary Path Finding Algorithm

    Fackyyj loves the challenge phase in TwosigmaCrap(TC). One day, he meet a task asking him to find sh ...

  2. HDU 4889 Scary Path Finding Algorithm

    其实这个题是抄的题解啦…… 题解给了一个图,按照那个图模拟一遍大概就能理解了. 题意: 有一段程序,给你一个C值(程序中某常量),让你构造一组数据,使程序输出"doge" 那段代码 ...

  3. Intersection between 2d conic in OpenCASCADE

    Intersection between 2d conic in OpenCASCADE eryar@163.com Abstract. OpenCASCADE provides the algori ...

  4. Intersection between a 2d line and a conic in OpenCASCADE

    Intersection between a 2d line and a conic in OpenCASCADE eryar@163.com Abstract. OpenCASCADE provid ...

  5. OpenCASCADE 3 Planes Intersection

    OpenCASCADE 3 Planes Intersection eryar@163.com Abstract. OpenCASCADE provides the algorithm to sear ...

  6. Class loading in JBoss AS 7--官方文档

    Class loading in AS7 is considerably different to previous versions of JBoss AS. Class loading is ba ...

  7. Matlab 中S-函数的使用 sfuntmpl

    function [sys,x0,str,ts,simStateCompliance] = sfuntmpl(t,x,u,flag) %SFUNTMPL General MATLAB S-Functi ...

  8. &lbrack;Swift&rsqb;LeetCode990&period; 等式方程的可满足性 &vert; Satisfiability of Equality Equations

    Given an array equations of strings that represent relationships between variables, each string equa ...

  9. 32 Profiling Go Programs 分析go语言项目

    Profiling Go Programs  分析go语言项目 24 June 2011 At Scala Days 2011, Robert Hundt presented a paper titl ...

随机推荐

  1. the request resource is not available

    form表单递交数据的问题 我的解决方法 将要访问的servlet地址写入form的action中 例如:访问地址为http://localhost:8080/Webprojectselfservic ...

  2. C&num;--GDI&plus;的LinearGradientBrush类

    命名空间:System.Drawing.Drawing2D LinearGradientBrush对象用颜色线性渐变填充图形.简言之,颜色渐变包含一种在两种指定的颜色之间渐变的颜色,渐变的方向是沿着指 ...

  3. Dictionary序列化和反序列化

    public class SerializeHelper { public static string XmlSerialize(List<CustomSearchEntity> obj) ...

  4. U盘安装Linux安装报错及解决方案

    导读 从网上看到了<Linux就该这么学>后,偏离软件行业多年的我下定决心回归!这篇文章是我这一个小白的亲身经历,希望能被采纳! 开始按照<Linux就该这么学>中所讲在自己的 ...

  5. Sharing

    To store English words, one method is to use linked lists and store a word letter by letter. To save ...

  6. 【转】vlc android 代码编译

    转自:http://blog.csdn.net/asircao/article/details/7734201 系统:ubuntu12.04代码:git://git.videolan.org/vlc- ...

  7. javascript第三课underfind和类型获取

    1.underfind一般发生于变量定义之后未赋值,因此变量的值就为underfind 2.var obj=new object(); 此时使用obj点,可以获取到obj对象的一些方法,使用alert ...

  8. LINQ To SQL 处理 DateTime&quest;

    LINQ To SQL 处理 DateTime? 类型 例子: 搜索栏含有最后扫描时间的日期(DateTime?)与多个其他条件(String) 现在需要写一个查询 : 查询符合最后扫描的日期的查询 ...

  9. JDK 8u131

    JDK 8u131 发布了.Java SE 8u131 包括重要的安全修复和bug修复.Oracle 强烈建议所有 JavaSE 8 用户升级到此版本.此次完整版本号为1.8.0_131-b11. J ...

  10. LAMP部署流水

    1.安装完成linux系统后,关闭防火墙: [root@localhost ~]# service iptables stop iptables: Setting chains to policy A ...