📄 newtonsolver.h
字号:
//////////////////////////////////////////////////////////////////////////// Copyright (c) 2000, Yusuke Miyao/// You may distribute under the terms of the Artistic License.////// <id>$Id: NewtonSolver.h,v 1.3 2003/04/24 08:29:36 kazama Exp $</id>/// <collection>Maximum Entropy Estimator</collection>/// <name>Solver.h</name>/// <overview>Equation solver</overview>/////////////////////////////////////////////////////////////////////////#ifndef Amis_NewtonSolver_h_#define Amis_NewtonSolver_h_#include <amis/configure.h>#include <amis/Real.h>AMIS_NAMESPACE_BEGIN///////////////////////////////////////////////////////////////////////// <classdef>/// <name>NewtonSolver</name>/// <overview>Equation solver</overview>/// <desc>/// Solver for polynomial equations with Newton's and bisection methods./// </desc>/// <body>template < class FeatureFreqArray >class NewtonSolver {public: typedef typename FeatureFreqArray::FeatureFreq FeatureFreq; static const int MAX_NEWTON_ITERATIONS = 200;protected: Real initialUpdate( const FeatureFreqArray& array, Real constant ) const { FeatureFreq order; Real freq; if ( array.highestOrder( order, freq ) ) { //cerr << "pow " << constant / freq << " " << 1.0 / order << std::endl; return pow( constant / freq, (Real)(1.0 / order) ); } else { AMIS_ABORT( "Cannot find the term of the highest order" ); } } /// Compute the initial value for newton's methodpublic: Real solveEquation( const FeatureFreqArray& array, Real constant, Real epsilon = 0.0, int max_iter = MAX_NEWTON_ITERATIONS ) const { Real lower_bound = 0.0; Real upper_bound = initialUpdate( array, constant ); Real current_x = upper_bound; for ( int i = 0; i < max_iter; ++i ) { Real polynomial = -constant; Real derivative = 0.0; array.computePolynomialDerivative( current_x, polynomial, derivative ); if ( fabs( polynomial ) <= epsilon ) break; current_x -= polynomial / derivative; // newton if ( current_x < lower_bound || upper_bound < current_x ) { //cerr << " out of range! " << current_x << std::endl; current_x = ( lower_bound + upper_bound ) * 0.5; // bisection } if ( ! finite( current_x ) ) break; if ( polynomial < 0.0 ) lower_bound = current_x; else upper_bound = current_x; if ( fabs( upper_bound - lower_bound ) <= epsilon ) break; } return current_x; } /// Solve the polynomial equation by newton's and bisection method};/// </body>/// </classdef>AMIS_NAMESPACE_END#endif // NewtonSolver_h_// end of NewtonSolver_h_⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -