posts - 393, comments - 566, trackbacks - 0, articles - 0

OpenCASCADE解非线性方程组

Posted on 2018-09-06 21:20 eryar 阅读(284) 评论(0)  编辑 收藏 引用 所属分类: 2.OpenCASCADE

OpenCASCADE解非线性方程组

eryar@163.com

 

Abstract. 在科学技术领域里常常提出求解非线性方程组的问题,例如,用非线性函数拟合实验数据问题、非线性网络问题、几何上的曲线曲面求交问题等。OpenCASCADE中有关于非线性方程组定义的类及其求解类,本文主要介绍如何在OpenCASCADE中定义非线性方程组,及对其进行求解。

Key Words. Function Set, Function Set Root, Newton Raphson Algorithm

1.Introduction

在科学技术领域里常常提出求解非线性方程组的问题,例如,用非线性函数拟合实验数据问题、非线性网络问题等。在几何造型中很多问题也可以利用非线性方程组来解决。如曲线的光顺,曲线求交、曲面求交、Blend造型问题等。

OpenCASCADE提供了非线性方程组的类math_FunctionSet,可以先从类图上来看看有哪些算法使用了这个类:

  

图1 曲线光顺包FaireCurve

 

图2 Blending Surface between two surfaces

感兴趣的同学可以自己打开OpenCASCADE的类索引文件查看。可以看到很多算法涉及到方程组的求解问题。本文主要介绍如何定义非线性方程组及对其进行求解。理解这些套路后,对math_FunctionSet相关的派生类及其用用途就会有个清晰的认识,便于对源码的理解。

2.Function Set Definition

设有非线性方程组

 

为实变量的非线方程函数。引入向量形式表示,引进记号:

 

于是非线性方程组可以简单记作:F(x)=0。我们的问题是寻求X使F(X)=0,这个X就是非线性方程组的解。

OpenCASCADE中使用类math_FunctionSet来表示方程组,这是个抽象类,定义了如下纯虚函数:

l NbVariables():变量的个数,即末知量的个数;

l NbEquations():方程的个数,即方程组中有几个方程;

l Value(const math_Vector&X, math_Vector& F):方程组的值,即代入变量每个方程的值;

3.Function Set Root Algorithm

解非线性方程组的牛顿法和解方程式的思路一样,要求方程有一阶导数。而非线性方程组即是要求有偏导数。由fi(x)偏导数作成的矩阵记为J(x)F’(x),称为F(x)Jacobi矩阵:

 

求解非线性方程组的牛顿法为:

 

其中xk为方程线的近似解向量。

 

OpenCASCADE中也提供了非线性方程组的求解类,如:math_FunctionSetRootmath_NewtonFunctionSetRoot。而使用这些类的输入都是要求具有一阶偏导数的线性方程组的定义math_FunctionSetWithDerivaties。这个类定义了具有一阶偏导数的非线性方程组,其纯虚函数除了前面说明的几个以外,还增加了如下两个:

l Derivatives(const math_Vector& X, math_Matrix& D):一阶偏导数值,即计算Jacobi矩阵;

l Values(const math_Vector& X, math_Vector& F, math_Matrix& D):计算方程的值及一阶偏导数矩阵Jacobi矩阵。

4.Code Example

下面给出一个具体的例子来说明这些类的用法。设有非线性方程组:

 

从几何上看其解就是圆心在原点,半径为2的圆与曲线的交点:

 

圆与曲线求交

下面我们使用OpenCASCADE来对上述问题进行求解。首先定义这个非线性方程组: 

#include <math_FunctionSet.hxx>
#include 
<math_FunctionSetWithDerivatives.hxx>
#include 
<math_FunctionSetRoot.hxx>
#pragma comment(lib, 
"TKernel.lib")
#pragma comment(lib, 
"TKMath.lib")
class MyFunctionSet : public math_FunctionSetWithDerivatives
{
public:
    
virtual Standard_Integer NbVariables() const
    {
        
return 2;
    }
    
virtual Standard_Integer NbEquations() const
    {
        
return 2;
    }
    
virtual Standard_Boolean Value(const math_Vector& X, math_Vector& F)
    {
        F(
1= X(1* X(1+ X(2* X(2- 4.0;
        F(
2= Pow(M_E, X(1)) + X(2- 1.0;
        
return Standard_True;
    }
    
virtual Standard_Boolean Derivatives(const math_Vector& X, math_Matrix& D)
    {
        
// matrix D is Jacobi matrix.
        D(11= 2.0 * X(1);
        D(
12= 2.0 * X(2);
        D(
21= Pow(M_E, X(1));
        D(
22= 1.0;
        
return Standard_True;
    }
    
virtual Standard_Boolean Values(const math_Vector& X, math_Vector& F, math_Matrix& D)
    {
        Value(X, F);
        Derivatives(X, D);
        
return Standard_True;
    }
private:
};
void test()
{
    MyFunctionSet aFunctionSet;
    math_FunctionSetRoot aSolver(aFunctionSet);
    math_Vector aStartingPoint(
12);
    
// 1. (1.0, 1.0)
    aStartingPoint(1= 1.0;
    aStartingPoint(
2= 1.0;
    aSolver.Perform(aFunctionSet, aStartingPoint);
    
if (aSolver.IsDone())
    {
        aSolver.Dump(std::cout);
    }
    
// 2. (1.0, -1.0)
    aStartingPoint(1= 1.0;
    aStartingPoint(
2= -1.0;
    aSolver.Perform(aFunctionSet, aStartingPoint);
    
if (aSolver.IsDone())
    {
        aSolver.Dump(std::cout);
    }
}
int main(int argc, char* argv[])
{
    test();
    
return 0;
}

上述代码先定义了带有一阶偏导数的非线性方程组类:MyFunctionSet,因为有两个变量及两个方程,再分别实现计算方程值及偏导数的虚函数。

然后使用类math_FunctionSetRoot来对方程组进行求解,求解的结果如下图所示:

 

非线性方程组求解结果

由图3可知,两个曲线相交有两个交点,但是使用类math_FunctionSetRoot一次只能计算一个解。从图4的计算结果还可以看出,初值的选择对解的影响很大,既影响计算结果,也影响迭代次数。

 5.Conclusion

综上所述,OpenCASCADEmath工具箱中提供了方程组的定义、求解功能。其中对非线性方程组求解使用的是Newton迭代法,所以要求方程组必须实现计算一阶偏导数的虚函数,即计算Jacobi矩阵。

OpenCASCADE类图中可以看出,方程组定义类用在了很多地方,所以理解上述对方程组的定义及解的用法,对其他使用这个派生类的地方更容易其源码。

 6.References

  1. 同济大学数学教研室高等数学 第四版高等教育出版社. 2004
  2. 易大义沈云宝李有法计算方法浙江大学出版社. 2002

  


为了方便大家在移动端也能看到我的博文和讨论交流,现已注册微信公众号,欢迎大家扫描下方二维码关注。
Shing Liu(eryar@163.com)

只有注册用户登录后才能发表评论。
【推荐】超50万行VC++源码: 大型组态工控、电力仿真CAD与GIS源码库
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理