以曲线弧长作为曲线方程的参数,Gauss型求积公式是一种高精度的数值积分公式

OpenCASCADE Gauss Integration

eryar@163.com

Abstract. Numerical integration is
the approximate computation of an integral using numerical techniques.
The numerical computation of an integral is sometimes called quadrature.
The most straightforward numerical integration technique uses the
Newton-Cotes formulas(also called quadrature formulas), which
approximate a function tabulated sequence of regularly spaced intervals
by various degree polynomials. If the functions are known analytically
instead of being tabulated at equally spaced intervals, the best
numerical method of integrations is called Gauss Integration(Gaussian
quadrature). By picking the abscissas at which to evaluate the function,
Gaussian quadrature produces the most accurate approximations possible.
In OpenCASCADE math package it implement the Gauss-Legendre integration.
So I will focus on the usage of the class in OpenCASCADE. 

Key Words. OpenCASCADE, Gauss
Integration, Gauss-Legendre, Numerical Analysis

1. Introduction

在正确和工程测算难点中,日常要计算一些定积分或微分,它们的精确值不可以算出或统计量太大,只好用数值的措施给出具有指定误差限的近似值。最直观的数值积分方法有牛顿-Cotes,其将积分区间等分之,并取分点为积分节点。那种做法尽管简化了计算,但却下落了所得公式的代数精度。 

Gauss型求积公式是一种高精度的数值积分公式。在求积节点数相同的情事下,即计算工作量相近的事态下,利用Gauss型求积公式往往可以拿到精确程序较高的积分结果,只是它在不相同距的无理数上测算被积函数。 

OpenCASCADE的math包中贯彻了Gauss-Legendre积分算法。本文主要介绍其利用格局,进而对其选择举办了然。

2. The Gauss-Legendre Integration

Gauss型求积公式是数值稳定的,且对个别闭区间上的连日函数,Gauss求积的数值随节点数目的充实而化为乌有到标准积分值。 

常用的Gauss型求积公式有Gauss-Legendre求积公式,Gauss-Chebyshev求积公式,Gauss-Laguerre求积公式和Gauss-Hermite求积公式等。 

对于一般区间[a,
b]上的Gauss型求积公式,可经过变量变换,由Gauss-Legendre求积公式得到: 

图片 1其中: 

图片 2

OpenCASCADE中对应的类有math_GaussSingleIntegration,首要完毕的函数为Perform(),总结进程如下: 

v 查表求得Gauss点及求积周全;

//Recuperation des points de Gauss dans le fichier GaussPoints.
  math::GaussPoints(Order,GaussP);
  math::GaussWeights(Order,GaussW);

v 依据Gauss-Legendre求积公式统计; 

 

// Changement de variable pour la mise a l'echelle [Lower, Upper] :
  xm = 0.5*(Upper + Lower);
  xr = 0.5*(Upper - Lower);
  Val = 0.;

  Standard_Integer ind = Order/2, ind1 = (Order+1)/2;
  if(ind1 > ind) { // odder case
    Ok1 = F.Value(xm, Val);
    if (!Ok1) return;
    Val *= GaussW(ind1);
  }
// Sommation sur tous les points de Gauss: avec utilisation de la symetrie.
  for (j = 1; j <= ind; j++) {
    dx = xr*GaussP(j);
    Ok1 = F.Value(xm-dx, F1);
    if(!Ok1) return;
    Ok1 = F.Value(xm+dx, F2);
    if(!Ok1) return;
    // Multiplication par les poids de Gauss.
    Standard_Real FT = F1+F2;
    Val += GaussW(j)*FT;  
  }
  // Mise a l'echelle de l'intervalle [Lower, Upper]
  Val *= xr;

比较Gauss-Legendre求积公式来了然上述代码如故相比较清晰的。上边给出使用此类的一个有血有肉实例: 

/*
*    Copyright (c) 2014 eryar All Rights Reserved.
*
*        File    : Main.cpp
*        Author  : eryar@163.com
*        Date    : 2014-09-11 20:46
*        Version : 1.0v
*
*    Description : Demo for Gauss-Legendre Integration usage.
*
*      Key words : OpenCascade, Gauss-Legendre Integration
*/

#define WNT
#include <math_Function.hxx>
#include <math_GaussSingleIntegration.hxx>

#pragma comment(lib, "TKernel.lib")
#pragma comment(lib, "TKMath.lib")

class Test_GaussFunction : public math_Function
{
public:
    virtual Standard_Boolean Value(const Standard_Real x, Standard_Real &y)
    {
        y = x;

        return Standard_True;
    }

private:
};

void TestGaussIntegration(void)
{
    Test_GaussFunction aFunction;
    math_GaussSingleIntegration aSolver(aFunction, 1, 10, 10);

    std::cout << aSolver << std::endl;
}

int main(int argc, char* argv[])
{
    TestGaussIntegration();

    return 0;
}

要害是从math_Function派生一个类来在虚函数Value()中重定义求积函数即可。上述实例中统计的是之类积分: 

图片 3

计量结果如下图所示: 

图片 4

Figure 2.1 Gauss-Legendre Integtation Result 

3. Application

由高级数学知识可见,积分的施用关键用以计算图形面积,体量及曲线的弧长,功等。 

积分在OpenCASCADE中的紧要运用有总计曲线长度,曲面面积及实体的体积等。如下图所示: 

图片 5

Figure 3.1 Compute Area of a Surface 

以身作则代码如下所示:

TopoDS_Shape S = BRepBuilderAPI_MakeFace(BSS, Precision::Confusion()).Face();

GProp_GProps System;
BRepGProp::SurfaceProperties(S,System);
gp_Pnt G = System.CentreOfMass ();
Standard_Real Area = System.Mass();
gp_Mat I = System.MatrixOfInertia();

 

 

4. Conclusion

OpenCASCADE中贯彻的Gauss-Legendre求积算法,由于是查表求得Gauss点及求积周到,所以测算速度快。唯一不足是对高斯点数有限制。 

综述,可见数值总结在OpenCASCADE中最主要职能。一个TKMath库也就是落成了一本《数值分析》课本中的一大半情节。所以有趣味的朋友可整合《数值分析》或《统计方法》之类的书本,来对OpenCASCADE的数学库TKMath举办驳斥联系实际的力透纸背驾驭。

5. References

  1. Wolfram MathWorld, Numerical Integration,  

http://mathworld.wolfram.com/NumericalIntegration.html

  1. 易大义,沈云宝,李有法编. 统计方法. 新疆大学出版社. 2002 

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

  3. 李庆杨,王能超,易大义.数值分析.华中地质大学出版社. 1986 

  4. 同济大学数学教研室. 高等数学(第四版). 高等教育出版社. 1996

PDF Version: OpenCASCADE Gauss
Integration

OpenCASCADE Curve Length
Calculation

OpenCASCADE 参数曲面面积

eryar@163.com

Abstract.
本文介绍了参数曲面的率先中坚公式,并利用曲面的第一骨干公式,结合OpenCASCADE中总计多重积分的类,对私下参数曲面的面积举行统计。

Key Words. Parametric Curve, Parametric Surface, Gauss Integration,
Global Properties

eryar@163.com

1.Introduction

小编们通晓一元函数y=f(x)的图像是一条曲线,二元函数z=f(x,y)的图像是一张曲面。但是,把曲线曲面表示成参数方程则越来越便利于探究,那种代表方法首先是由欧洲瑞士联邦地数学家Euler引进的。例如,在空间中的一条曲线可以象征为多个一元函数:

X=x(t), Y=y(t), Z=z(t)

在向量的定义出现后,空间中的一条曲线可以自然地表示为一个一元向量函数:

r=r(t)=(x(t), y(t), z(t))

用向量函数来表示曲线和曲面后,使曲线曲面一些量的盘算方法相比统一。如曲线可以象征为一元向量函数,曲面可以代表为二元向量函数。

正文结合OpenCASCADE来介绍参数曲线曲面积分分别计算曲线弧长和曲面的面积。结合《微分几何》来更好地知道曲线曲面相关知识。

Abstract. The natural parametric equations of a curve are parametric
equations that represent the curve in terms of a coordinate-independent
parameter, generally arc length s, instead of an arbitray variable like
t or u. According to the natural equations, the curve length is the
integration of the curve parametric equation’s derivation. So the core
algorithm for curve length calculation is the numerical integration
method. OpenCASCADE use Gauss-Legendre to calculate the integration for
single variable and multiple variables. Because of curve in OpenCASCADE
is single variable function, so can use the Gauss-Legendre to calculate
the arc length for the curve. 

2.Curve Natural Parametric Equations

设曲线C的参数方程是r=r(t),命:

图片 6

则s是该曲线的一个不变量,即它与上空中的坐标系的抉择非亲非故,也与该曲线的参数变换无关。前者是因为在笛Carl直角坐标变换下,切向量的长短|r’(t)|是不变的,故s不变。关于后世可以经过积分的变量的交替来验证,设参数变换是:

图片 7

并且

图片 8

因此

图片 9

依据积分的变量替换公式有:

图片 10

不变量s的几何意义就是曲线段的弧长。那表达曲线参数t可以是随机的,但挑选不相同的参数拿到的参数方程会有两样,不过曲线段的弧长是不变的。以曲线弧长作为曲线方程的参数,那样的方程称为曲线的当然参数方程Natural
Parametric Equations。

由曲线的参数方程可见,曲线弧长的计算公式为:

图片 11

几何意义就是在种种微元处的切向量的长短求和。

Key Words. OpenCASCADE, The Gauss-Legendre Integration, Parametric
Curve,  

3.Gauss Integration for Arc Length

曲线弧长的测算就是一元函数的积分。OpenCASCADE中是如何总括任意曲线弧长的啊?直接找到相关的源码列举如下:(在类CPnts_AbscissaPoint中)

// auxiliary functions to compute the length of the derivative

static Standard_Real f3d(const Standard_Real X, const Standard_Address C)

{

  gp_Pnt P;

  gp_Vec V;

((Adaptor3d_Curve*)C)->D1(X,P,V);

return V.Magnitude();

}

static Standard_Real f2d(const Standard_Real X, const Standard_Address C)

{

  gp_Pnt2d P;

  gp_Vec2d V;

((Adaptor2d_Curve2d*)C)->D1(X,P,V);

return V.Magnitude();

}

//==================================================================

//function : Length

//purpose  : 3d with parameters

//==================================================================

Standard_Real CPnts_AbscissaPoint::Length(const Adaptor3d_Curve& C,

const Standard_Real U1,

const Standard_Real U2)

{

  CPnts_MyGaussFunction FG;

//POP pout WNT

  CPnts_RealFunction rf = f3d;

  FG.Init(rf,(Standard_Address)&C);

//  FG.Init(f3d,(Standard_Address)&C);

  math_GaussSingleIntegration TheLength(FG, U1, U2, order(C));

if (!TheLength.IsDone()) {

throw Standard_ConstructionError();

}

return Abs(TheLength.Value());

}

//==================================================================

//function : Length

//purpose  : 2d with parameters

//==================================================================

Standard_Real CPnts_AbscissaPoint::Length(const Adaptor2d_Curve2d& C,

const Standard_Real U1,

const Standard_Real U2)

{

  CPnts_MyGaussFunction FG;

//POP pout WNT

  CPnts_RealFunction rf = f2d;

  FG.Init(rf,(Standard_Address)&C);

//  FG.Init(f2d,(Standard_Address)&C);

  math_GaussSingleIntegration TheLength(FG, U1, U2, order(C));

if (!TheLength.IsDone()) {

throw Standard_ConstructionError();

}

return Abs(TheLength.Value());

}

 

上述代码的趣味是平素对曲线的一阶导数的长度求积分,即是弧长。OpenCASCADE的代码写得有点难懂,依照意思把对三维曲线求弧长的代码改写下,更有利领会:

//! Function for curve length evaluation.

class math_LengthFunction : public math_Function

{

public:

    math_LengthFunction(const Handle(Geom_Curve)& theCurve)

: myCurve(theCurve)

{

}

virtual Standard_Boolean Value(const Standard_Real X, Standard_Real& F)

{

        gp_Pnt aP;

        gp_Vec aV;

        myCurve->D1(X, aP, aV);

        F = aV.Magnitude();

return Standard_True;

}

private:

    Handle(Geom_Curve) myCurve;

};

 

          The Natural Parametric Equations,  Arc Length, 

4.First Fundamental Form of a Surface

曲面参数方程是个二元向量函数。依照《微分几何》中曲面的第一主干公式(First
Fundamental Form of a Surface)可见,曲面上曲线的表明式为:

r=r(u(t), v(t)) = (x(t), y(t), z(t))

若以s代表曲面上曲线的弧长,则由复合函数求导公式可得弧长微分公式:

图片 12

在古典微分几何中,上式称为曲面的第一大旨公式,E、F、G称为第一基本量。在曲面上,每一点的率先基本量与参数化无关。

选取曲面第一为君王式可以用于总括曲面的面积。参数曲面上与u,v参数空间的要素dudv对应的面积元为:

图片 13

由参数曲面法向的测算可见,曲面的面积成分即为u,v方向上的偏导数的乘积的模。

图片 14

其几何意义可以知晓为参数曲面的面积微元是由u,v方向的偏导数的向量围成的一个四边形的面积,则全体曲面的面积即是对面积成分求积分。由于参数曲面有七个参数,所以若要总结曲面的面积,只须求对面积元素计算二重积分即可。

  

5.Gauss Integration for Area

OpenCASCADE的math包中提供了多重积分的计算类math_GaussMultipleIntegration,由类名可见积分算法接纳了Gauss积分算法。上面通过切实代码来表明OpenCASCADE中计算曲面积分的经过。要统计积分,先要定义被积函数。因为参数曲面与参数曲线差距,参数曲线唯有一个参数,而参数曲面有多个参数,所以是一个多元函数。

//! 2D variable function for surface area evaluation.

class math_AreaFunction : public math_MultipleVarFunction

{

public:

    math_AreaFunction(const Handle(Geom_Surface)& theSurface)

: mySurface(theSurface)

{

}

virtual Standard_Integer NbVariables() const

{

return 2;

}

virtual Standard_Boolean Value(const math_Vector& X, Standard_Real& Y)

{

        gp_Pnt aP;

        gp_Vec aDu;

        gp_Vec aDv;

        mySurface->D1(X(1), X(2), aP, aDu, aDv);

        Y = aDu.Crossed(aDv).Magnitude();

return Standard_True;

}

private:

    Handle(Geom_Surface) mySurface;

};

 

鉴于参数曲面是多元函数,所以从类math_MultipleVarFunction派生,并在虚函数NbVariables()中验证有三个变量。在虚函数Value()中统计面积成分的值,即依照曲面第一宗旨公式中面积成分的概念,对参数曲面求一阶导数,总结五个偏导数向量的叉乘的模。

有了被积函数,只需求在定义域对其总计二重积分,相应代码如下所示:

void evalArea(const Handle(Geom_Surface)& theSurface, const math_Vector& theLower, const math_Vector& theUpper)

{

    math_IntegerVector aOrder(1, 2, math::GaussPointsMax());

    math_AreaFunction aFunction(theSurface);

    math_GaussMultipleIntegration anIntegral(aFunction, theLower, theUpper, aOrder);

if (anIntegral.IsDone())

{

        anIntegral.Dump(std::cout);

}

}

 

透过theLower和theUpper指定定义域,由于接纳了Gauss-Legendre算法统计二重积分,所以必要指定阶数,且阶数越高积分结果精度越高,那里运用了OpenCASCADE中最高的阶数。

上边通过对大旨曲面的面积统计来验证结果的不利,并将总括结果和OpenCASCADE中总括面积的类BRepGProp::SurfaceProperties()结果进行自查自纠。

1. Introduction

6.Elementary Surface Area Test

上面通过对OpenCASCADE中多少个初等曲面的面积进行测算,代码如下所示:

/*
Copyright(C) 2017 Shing Liu(eryar@163.com)

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files(the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and / or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions :

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/

#include <gp_Pnt.hxx>
#include <gp_Vec.hxx>

#include <math.hxx>
#include <math_Function.hxx>
#include <math_MultipleVarFunction.hxx>
#include <math_GaussMultipleIntegration.hxx>

#include <Geom_Plane.hxx>
#include <Geom_ConicalSurface.hxx>
#include <Geom_CylindricalSurface.hxx>
#include <Geom_SphericalSurface.hxx>
#include <Geom_ToroidalSurface.hxx>
#include <Geom_BSplineSurface.hxx>
#include <Geom_RectangularTrimmedSurface.hxx>

#include <GeomConvert.hxx>

#include <GProp_GProps.hxx>

#include <TopoDS_Face.hxx>

#include <BRepGProp.hxx>

#include <BRepBuilderAPI_MakeFace.hxx>

#pragma comment(lib, "TKernel.lib")
#pragma comment(lib, "TKMath.lib")

#pragma comment(lib, "TKG2d.lib")
#pragma comment(lib, "TKG3d.lib")
#pragma comment(lib, "TKGeomBase.lib")
#pragma comment(lib, "TKGeomAlgo.lib")

#pragma comment(lib, "TKBRep.lib")
#pragma comment(lib, "TKTopAlgo.lib")


//! 2D variable function for surface area evaluation.
class math_AreaFunction : public math_MultipleVarFunction
{
public:
    math_AreaFunction(const Handle(Geom_Surface)& theSurface)
        : mySurface(theSurface)
    {

    }

    virtual Standard_Integer NbVariables() const
    {
        return 2;
    }

    virtual Standard_Boolean Value(const math_Vector& X, Standard_Real& Y)
    {
        gp_Pnt aP;
        gp_Vec aDu;
        gp_Vec aDv;

        Standard_Real E = 0.0;
        Standard_Real F = 0.0;
        Standard_Real G = 0.0;

        mySurface->D1(X(1), X(2), aP, aDu, aDv);

        E = aDu.Dot(aDu);
        F = aDu.Dot(aDv);
        G = aDv.Dot(aDv);

        Y = Sqrt(E * G - F * F);

        //Y = aDu.Crossed(aDv).Magnitude();

        return Standard_True;
    }

private:
    Handle(Geom_Surface) mySurface;
};

void evalArea(const Handle(Geom_Surface)& theSurface, const math_Vector& theLower, const math_Vector& theUpper)
{
    math_IntegerVector aOrder(1, 2, math::GaussPointsMax());

    math_AreaFunction aFunction(theSurface);
    math_GaussMultipleIntegration anIntegral(aFunction, theLower, theUpper, aOrder);
    if (anIntegral.IsDone())
    {
        anIntegral.Dump(std::cout);
    }
}

void evalArea(const Handle(Geom_BoundedSurface)& theSurface)
{
    math_IntegerVector aOrder(1, 2, math::GaussPointsMax());

    math_Vector aLower(1, 2, 0.0);
    math_Vector aUpper(1, 2, 0.0);

    theSurface->Bounds(aLower(1), aUpper(1), aLower(2), aUpper(2));

    math_AreaFunction aFunction(theSurface);
    math_GaussMultipleIntegration anIntegral(aFunction, aLower, aUpper, aOrder);
    if (anIntegral.IsDone())
    {
        anIntegral.Dump(std::cout);
    }
}

void testFace(const TopoDS_Shape& theFace)
{
    GProp_GProps aSurfaceProps;

    BRepGProp::SurfaceProperties(theFace, aSurfaceProps);

    std::cout << "Face area: " << aSurfaceProps.Mass() << std::endl;
}

void testPlane()
{
    std::cout << "====== Test Plane Area =====" << std::endl;
    Handle(Geom_Plane) aPlaneSurface = new Geom_Plane(gp::XOY());

    math_Vector aLower(1, 2);
    math_Vector aUpper(1, 2);

    // Parameter U range.
    aLower(1) = 0.0;
    aUpper(1) = 2.0;

    // Parameter V range.
    aLower(2) = 0.0;
    aUpper(2) = 3.0;

    evalArea(aPlaneSurface, aLower, aUpper);

    // Convert to BSpline Surface.
    Handle(Geom_RectangularTrimmedSurface) aTrimmedSurface =
        new Geom_RectangularTrimmedSurface(aPlaneSurface, aLower(1), aUpper(1), aLower(2), aUpper(2));

    Handle(Geom_BSplineSurface) aBSplineSurface = GeomConvert::SurfaceToBSplineSurface(aTrimmedSurface);
    evalArea(aBSplineSurface);

    // Test Face.
    TopoDS_Face aFace = BRepBuilderAPI_MakeFace(aTrimmedSurface, Precision::Confusion()).Face();
    testFace(aFace);

    aFace = BRepBuilderAPI_MakeFace(aBSplineSurface, Precision::Confusion()).Face();
    testFace(aFace);
}

void testCylinder()
{
    std::cout << "====== Test Cylinder Area =====" << std::endl;
    Handle(Geom_CylindricalSurface) aCylindrialSurface = new Geom_CylindricalSurface(gp::XOY(), 1.0);

    math_Vector aLower(1, 2);
    math_Vector aUpper(1, 2);

    aLower(1) = 0.0;
    aUpper(1) = M_PI * 2.0;

    aLower(2) = 0.0;
    aUpper(2) = 3.0;

    evalArea(aCylindrialSurface, aLower, aUpper);

    // Convert to BSpline Surface.
    Handle(Geom_RectangularTrimmedSurface) aTrimmedSurface =
        new Geom_RectangularTrimmedSurface(aCylindrialSurface, aLower(1), aUpper(1), aLower(2), aUpper(2));

    Handle(Geom_BSplineSurface) aBSplineSurface = GeomConvert::SurfaceToBSplineSurface(aTrimmedSurface);
    evalArea(aBSplineSurface);

    // Test Face.
    TopoDS_Face aFace = BRepBuilderAPI_MakeFace(aTrimmedSurface, Precision::Confusion()).Face();
    testFace(aFace);

    aFace = BRepBuilderAPI_MakeFace(aBSplineSurface, Precision::Confusion()).Face();
    testFace(aFace);
}

void testSphere()
{
    std::cout << "====== Test Sphere Area =====" << std::endl;
    Handle(Geom_SphericalSurface) aSphericalSurface = new Geom_SphericalSurface(gp::XOY(), 1.0);

    math_Vector aLower(1, 2);
    math_Vector aUpper(1, 2);

    aSphericalSurface->Bounds(aLower(1), aUpper(1), aLower(2), aUpper(2));

    evalArea(aSphericalSurface, aLower, aUpper);

    // Convert to BSpline Surface.
    Handle(Geom_BSplineSurface) aBSplineSurface = GeomConvert::SurfaceToBSplineSurface(aSphericalSurface);
    evalArea(aBSplineSurface);

    // Test Face.
    TopoDS_Face aFace = BRepBuilderAPI_MakeFace(aSphericalSurface, Precision::Confusion()).Face();
    testFace(aFace);

    aFace = BRepBuilderAPI_MakeFace(aBSplineSurface, Precision::Confusion()).Face();
    testFace(aFace);
}

void test()
{
    testPlane();

    testSphere();

    testCylinder();
}

int main(int argc, char* argv[])
{
    test();

    return 0;
}

 

算算结果如下图所示:

图片 15

上述代码统计了曲面的面积,再将曲面转换成B样条曲面,再拔取算法计算面积。再将曲面和更换的B样条曲素不相识成拓朴面,利用OpenCASCADE中统计曲面面积效能举行对照。使用自定义函数math_AreaFunction利用多重积分类总括的结果与OpenCASCADE中统计曲面面积的值是一律的。当把曲面转换成B样条曲面后,OpenCASCADE总结的曲面面积偏大。

对于同样条曲线,选取的参数差别其表明格局也大相径庭。而曲线本身的弧长是曲线自己的不变量,它与坐标系的选料无关,由此取曲线的弧长作为参数来钻探曲线具有特别首要的意思。以曲线弧长作为曲线方程的参数,那样的方程称为曲线的当然参数方程。由曲线的当然参数方程可见曲线的弧长是曲线参数方程导数的积分。所以总括曲线弧长的骨干算法成了总计函数的积分了。 

7.Conclusion

在攻读《高等数学》的积分时,其重点的一个施用就是计算弧长、面积和容积等。学习高数抽象概念时,总会问学了高数有怎么着用?就从统计机图形方面来看,可以采纳数学工具对任意曲线求弧长,对随意曲面总计面积等,更具一般性。

透过自定义被积函数再使用积分算法来计算任意曲面的面积,将舌战与履行结合起来了。即将曲面的第一为国君式与现实的代码甚至可以动用OpenCASCADE生成对应的图片,这样抽象的辩解就直观了,更有益于掌握相应的定义。

OpenCASCADE中数学计算库TKMath中有一种高精度总计积分的算法Gauss-Legendre积分法,可对单变量和多变量的函数举行积分计算。由于OpenCASCADE中曲线是单变量的参数方程,所以可以用Gauss-Legendre来计量积分,进而赢得曲线的弧长。如下图所示: 

8.References

1.朱心雄. 自由曲线曲面造型技术. 科学出版社. 2000

2.陈维桓. 微分几何. 上海大学出版社. 2006

3.同济高等高校数学教研室. 高等数学. 高等教育出版社. 1996

图片 16

Figure 1.1 Length of Curve 

正文对曲线的自然参数方程的定义举行验证,并简要介绍Gauss-Legendre算法及其在OpenCASCADE中的应用,即求曲线的弧长。

2. The Natural Parametric
Equations

对于同一曲线若拔取的参数差距,则其表明式亦不相同,故用坐标系商量曲线时,具有人为的属性,而曲线本人的弧长则曲线的不变量,它与坐标系的取舍毫无干系。由此,大家取曲线本身的弧长作为参数来研商曲线的内在性质。无论是理论探索依旧实用商讨,弧长参数化都有很重点的含义。 

给定空间曲线r,在其下车取一点P0(x0, y0,
z0)作为计量弧长的源点。应用弧长积分公式,即可测算该曲线上任意点P(x,y,z)到P0之间的弧长。因而,曲线上点的岗位与该点处的弧长是逐一对应的,如下图所示: 

图片 17

Figure 2.1 Arc length parameterization 

即曲线上点的坐标是以弧长为参数的函数: 

图片 18

由其矢量方程可见,曲线是弧长为参数的矢函数,大家将其誉为弧长参数化(Arc
length parameterization)。弧长称为自然参数(Natural
Parameter),曲线方程称为自然参数方程(Natural Parametric
Equations)。现研商曲线的本来参数方程与参数方程的关联:已知曲线的矢量方程为: 

图片 19

则弧长的微分和积分公式分别为: 

图片 20

s(u)即为与参数u0和u1对应的两点P0和P之间的弧长。弧长参数化是一类重点的定义,可是只要用上式来总括弧长相比繁琐,可以经过累加弦长方法来近似计算。 

给定曲线: 

图片 21

图片 22是参数轴上的一个等距分划;又令图片 23为曲线上与参数ui对应的点列,则可用下式总括弦长: 

图片 24

图片 25

Figure 2.2 Parametric Curve 

其中校参数轴等分愈多,则求得的曲线的弧长越规范。 

3. The Gauss-Legendre
Integration

上述累加弦长的乘除方法应该是牛顿-Cotes积分计算法,牛顿-Cotes就是将积分区间等分,并取分点为求积节点。不难看到当积分区间较大时,直接使用Newton-Cotes公式所得积分近似值的精度是很难有限协助的。 

至于积分的数值算法有很多,其中Gauss-Legendre求积公式具有统计工作量小,所得近似值精确度高等优点,是一种高精度的求积公式。情势如下所示的求积公式: 

图片 26

代数精度达到了2n+1,则称它为高斯型求积公式,并称相应的求积节点x0, x1,
… xn为高斯点(Gauss Point)。Ak为求积周到。 

OpenCASCADE的主旨模块的数学库TKMath中有Gauss-Legendre求积算法,可用来对单变量和多变量的函数进行积分计算,对应的类分别为:math_GaussSingleIntegration和math_GaussMultipleIntegration。计算所需的高斯点及周详通过查表的法门得到,数据以数组的花样在文书math.cxx中罗列出,如下图所示: 

图片 27

Figure 3.1 Gauss Point for Gauss-Legendre Integration 

Gauss-Legendre积分用处之一就是基于曲线的本来参数方程总结曲线的弧长,代码完成如下所示: 

//=====================================================================
//function : Length
//purpose  : 3d with parameters
//=====================================================================

Standard_Real CPnts_AbscissaPoint::Length(const Adaptor3d_Curve& C,
                      const Standard_Real U1,
                      const Standard_Real U2) 
{
  CPnts_MyGaussFunction FG;
//POP pout WNT
  CPnts_RealFunction rf = f3d;
  FG.Init(rf,(Standard_Address)&C);
//  FG.Init(f3d,(Standard_Address)&C);
  math_GaussSingleIntegration TheLength(FG, U1, U2, order(C));
  if (!TheLength.IsDone()) {
    Standard_ConstructionError::Raise();
  }
  return Abs(TheLength.Value());
}

依照上述求曲线弧长的代码可见,只要求给定曲线及其求积区间,即可测算出此区间内曲线的弧长。因为使用直接求积计算,所得曲线的弧长值的精度照旧很高的。 

高斯型求积公式是一种高精度的求积公式。在求积节点数相同,即计算工作量相近的意况下,利用高斯型求积公式往往可以拿走纯粹程序较高的积分近似值。不过,它必须在差别距的节点上总计被积函数的值,而且当节罗列改变时,所用数据都要重复查表总结。 

4. Conclusion

曲线自个儿的弧长是曲线本人的不变量,它与坐标系的选用毫无干系,由此取曲线的弧长作为参数来切磋曲线具有非常重大的意义。以曲线弧长作为曲线方程的参数,那样的方程称为曲线的自然参数方程。由曲线的本来参数方程可见曲线的弧长是曲线参数方程导数的积分。所以计算曲线弧长的主导算法成了计算函数的积分了。 

在《计算办法》、《数值分析》等读本中都至于于求积分的算法,在OpenCASCADE中落成了Gauss-Legendre求积算法。Gauss-Legendre求积算法是一种高精度的求积方法。所以依照曲线的本来参数方程可见曲线弧长就是对曲线参数方程导数的求积分。所以拔取高斯求积法可以拿到曲线弧长较精确值。 

对数值积分的切实可行算法感兴趣的爱侣,可以参见《总括方法》、《数值分析》、《数值总结》等有关书籍。 

5. References

  1. 朱心雄. 自由曲线曲面造型技术. 科学出版社. 2008 

  2. 王仁宏, 李崇君, 朱春钢. 总计几何教程. 科学出版社. 2008 

  3. 易大义, 沈云宝, 李有法. 统计方法. 西藏学院出版社. 2002 

  4. 易大义, 陈道琦. 数值分析引论. 甘肃大学出版社. 1998 

  5. http://mathworld.wolfram.com/NaturalParametricEquations.html

  6. http://mathworld.wolfram.com/Legendre-GaussQuadrature.html

  7. http://en.wikipedia.org/wiki/Gaussian_quadrature

  8. http://pomax.github.io/bezierinfo/legendre-gauss.html

 

PDF Version: OpenCASCADE Curve Length
Calculation

相关文章