# Apply Newton Method to Find Extrema in OPEN CASCADE

### Apply Newton Method to Find Extrema in OPEN CASCADE

eryar@163.com

Abstract. In calculus, Newton’s method is used for finding the roots of a function. In optimization, Newton’s method is applied to find the roots of the derivative. OPEN CASCADE implement Newton method to find the extrema for a multiple variables function, such as find the extrema point for a curve and a surface.

Key Words. Nonlinear Programming, Newton Method, Extrema, OPEN CASCADE

1. Introduction

Figure 1.1 A Curve and A Surface

2.Construct Function

```//======================================================================
//function : value
//purpose  :
//======================================================================
void Extrema_GlobOptFuncCS::value(Standard_Real cu,
Standard_Real su,
Standard_Real sv,
Standard_Real &F)
{
F = myC->Value(cu).SquareDistance(myS->Value(su, sv));
}```

```//======================================================================
//purpose  :
//======================================================================
Standard_Real su,
Standard_Real sv,
math_Vector &G)
{
gp_Pnt CD0, SD0;
gp_Vec CD1, SD1U, SD1V;

myC->D1(cu, CD0, CD1);
myS->D1(su, sv, SD0, SD1U, SD1V);

G(1) = + (CD0.X() - SD0.X()) * CD1.X()
+ (CD0.Y() - SD0.Y()) * CD1.Y()
+ (CD0.Z() - SD0.Z()) * CD1.Z();
G(2) = - (CD0.X() - SD0.X()) * SD1U.X()
- (CD0.Y() - SD0.Y()) * SD1U.Y()
- (CD0.Z() - SD0.Z()) * SD1U.Z();
G(3) = - (CD0.X() - SD0.X()) * SD1V.X()
- (CD0.Y() - SD0.Y()) * SD1V.Y()
- (CD0.Z() - SD0.Z()) * SD1V.Z();
}```

```//======================================================================
//function : hessian
//purpose  :
//======================================================================
void Extrema_GlobOptFuncCS::hessian(Standard_Real cu,
Standard_Real su,
Standard_Real sv,
math_Matrix &H)
{
gp_Pnt CD0, SD0;
gp_Vec CD1, SD1U, SD1V, CD2, SD2UU, SD2UV, SD2VV;

myC->D2(cu, CD0, CD1, CD2);
myS->D2(su, sv, SD0, SD1U, SD1V, SD2UU, SD2VV, SD2UV);

H(1,1) = + CD1.X() * CD1.X()
+ CD1.Y() * CD1.Y()
+ CD1.Z() * CD1.Z()
+ (CD0.X() - SD0.X()) * CD2.X()
+ (CD0.Y() - SD0.Y()) * CD2.Y()
+ (CD0.Z() - SD0.Z()) * CD2.Z();

H(1,2) = - CD1.X() * SD1U.X()
- CD1.Y() * SD1U.Y()
- CD1.Z() * SD1U.Z();

H(1,3) = - CD1.X() * SD1V.X()
- CD1.Y() * SD1V.Y()
- CD1.Z() * SD1V.Z();

H(2,1) = H(1,2);

H(2,2) = + SD1U.X() * SD1U.X()
+ SD1U.Y() * SD1U.Y()
+ SD1U.Z() * SD1U.Z()
- (CD0.X() - SD0.X()) * SD2UU.X()
- (CD0.Y() - SD0.Y()) * SD2UU.Y()
- (CD0.Z() - SD0.Z()) * SD2UU.Z();

H(2,3) = + SD1U.X() * SD1V.X()
+ SD1U.Y() * SD1V.Y()
+ SD1U.Z() * SD1V.Z()
- (CD0.X() - SD0.X()) * SD2UV.X()
- (CD0.Y() - SD0.Y()) * SD2UV.Y()
- (CD0.Z() - SD0.Z()) * SD2UV.Z();

H(3,1) = H(1,3);

H(3,2) = H(2,3);

H(3,3) = + SD1V.X() * SD1V.X()
+ SD1V.Y() * SD1V.Y()
+ SD1V.Z() * SD1V.Z()
- (CD0.X() - SD0.X()) * SD2VV.X()
- (CD0.Y() - SD0.Y()) * SD2VV.Y()
- (CD0.Z() - SD0.Z()) * SD2VV.Z();
}```

H(2,1)=H(1,2)

H(3,1)=H(1,3)

H(3,2)=H(2,3)

3.Newton’s Method

A. 给定初始点，及精度；

B. 计算函数f(xk)的一阶导数（梯度），二阶导数（Hessian Matrix）：若|梯度|<精度，则停止迭代，输出近似极小点；否则转C；

C. 根据Newton迭代公式，计算x(k+1)；

```/*
*
*           File : main.cpp
*         Author : Shing Liu(eryar@163.com)
*           Date : 2015-12-05 21:00
*
*    Description : Learn Newton\'s Method for multiple variables
*                  function.
*/

#define WNT
#include <TColgp_Array1OfPnt.hxx>
#include <TColgp_Array2OfPnt.hxx>

#include <math_NewtonMinimum.hxx>

#include <GeomTools.hxx>
#include <BRepTools.hxx>

#include <GC_MakeSegment.hxx>

#include <Extrema_GlobOptFuncCS.hxx>

#include <GeomAPI_PointsToBSpline.hxx>
#include <GeomAPI_PointsToBSplineSurface.hxx>

#include <BRepBuilderAPI_MakeEdge.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, \"TKBRep.lib\")
#pragma comment(lib, \"TKGeomBase.lib\")
#pragma comment(lib, \"TKGeomAlgo.lib\")
#pragma comment(lib, \"TKTopAlgo.lib\")

void testNewtonMethod(void)
{
// approximate curve from the points
TColgp_Array1OfPnt aCurvePoints(1, 5);
aCurvePoints.SetValue(1, gp_Pnt(0.0, 0.0, -2.0));
aCurvePoints.SetValue(2, gp_Pnt(1.0, 2.0, 2.0));
aCurvePoints.SetValue(3, gp_Pnt(2.0, 3.0, 3.0));
aCurvePoints.SetValue(4, gp_Pnt(4.0, 3.0, 4.0));
aCurvePoints.SetValue(5, gp_Pnt(5.0, 5.0, 5.0));

GeomAPI_PointsToBSpline aCurveApprox(aCurvePoints);

// approximate surface from the points.
TColgp_Array2OfPnt aSurfacePoints(1, 5, 1, 5);
aSurfacePoints(1, 1) = gp_Pnt(-4,-4,5);
aSurfacePoints(1, 2) = gp_Pnt(-4,-2,5);
aSurfacePoints(1, 3) = gp_Pnt(-4,0,4);
aSurfacePoints(1, 4) = gp_Pnt(-4,2,5);
aSurfacePoints(1, 5) = gp_Pnt(-4,4,5);

aSurfacePoints(2, 1) = gp_Pnt(-2,-4,4);
aSurfacePoints(2, 2) = gp_Pnt(-2,-2,4);
aSurfacePoints(2, 3) = gp_Pnt(-2,0,4);
aSurfacePoints(2, 4) = gp_Pnt(-2,2,4);
aSurfacePoints(2, 5) = gp_Pnt(-2,5,4);

aSurfacePoints(3, 1) = gp_Pnt(0,-4,3.5);
aSurfacePoints(3, 2) = gp_Pnt(0,-2,3.5);
aSurfacePoints(3, 3) = gp_Pnt(0,0,3.5);
aSurfacePoints(3, 4) = gp_Pnt(0,2,3.5);
aSurfacePoints(3, 5) = gp_Pnt(0,5,3.5);

aSurfacePoints(4, 1) = gp_Pnt(2,-4,4);
aSurfacePoints(4, 2) = gp_Pnt(2,-2,4);
aSurfacePoints(4, 3) = gp_Pnt(2,0,3.5);
aSurfacePoints(4, 4) = gp_Pnt(2,2,5);
aSurfacePoints(4, 5) = gp_Pnt(2,5,4);

aSurfacePoints(5, 1) = gp_Pnt(4,-4,5);
aSurfacePoints(5, 2) = gp_Pnt(4,-2,5);
aSurfacePoints(5, 3) = gp_Pnt(4,0,5);
aSurfacePoints(5, 4) = gp_Pnt(4,2,6);
aSurfacePoints(5, 5) = gp_Pnt(4,5,5);

GeomAPI_PointsToBSplineSurface aSurfaceApprox(aSurfacePoints);

// construct the function.

math_Vector aStartPoint(1, 3, 0.2);
math_NewtonMinimum aSolver(aFunction, aStartPoint);
aSolver.Perform(aFunction, aStartPoint);

if (aSolver.IsDone())
{
aSolver.Dump(std::cout);
math_Vector aLocation = aSolver.Location();

GC_MakeSegment aSegmentMaker(aPoint1, aPoint2);

BRepBuilderAPI_MakeEdge anEdgeMaker(aSegmentMaker.Value());
BRepTools::Write(anEdgeMaker.Shape(), \"d:/tcl/min.brep\");
}

GeomTools::Dump(aCurveApprox.Curve(), std::cout);
GeomTools::Dump(aSurfaceApprox.Surface(), std::cout);

BRepBuilderAPI_MakeEdge anEdgeMaker(aCurveApprox.Curve());
BRepBuilderAPI_MakeFace aFaceMaker(aSurfaceApprox.Surface(), Precision::Approximation());

BRepTools::Write(anEdgeMaker.Shape(), \"d:/tcl/edge.brep\");
BRepTools::Write(aFaceMaker.Shape(), \"d:/tcl/face.brep\");

// need release memory for the adaptor surface pointer manually.
// whereas do not need release memory for the adaptor curve.
// because it is mamanged by handle.
}

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

return 0;
}```

Figure 3.1 The minimum between a curve and a surface

4.Conclusion

5. References

3. Shing Liu. OpenCASCADE Root-Finding Algorithm

5. 同济大学数学教研室. 高等数学. 高等教育出版社. 1996

6. 同济大学应用数学系. 线性代数. 高等教育出版社. 2003

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

8. 《运筹学》教材编写组. 运筹学. 清华大学出版社. 2012

9. 何坚勇. 最优化方法. 清华大学出版社. 2007

10. 杨庆之. 最优化方法. 科学出版社. 2015

11. 王宜举, 修乃华. 非线性最优化理论与方法. 科学出版社. 2012

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