The algorithm itself is nothing special: Newton's method. Implementing it in a general and yet efficient manner in C++, however, presents a design choice.
I was tempted to accept a template argument and call that class's operator() member function.
template < typename Real, typename Functor >
int fsolve(matrix &x, const Functor &F, Real precision, int maxIterations) {
..
while () {
matrix fx = F(x); // operator()
}
..
}
This, however, would require a separate class for each constraint function. In my case of coupled systems, I need to determine roots for two distinct constraint functions that shared both constant state and analytically computed state. Additionally, the zeros of the first function would be needed as input to the second function.
Instead, I opted for the delegate design pattern. The functor includes both the object as well as references to member functions.
template <typename Real, typename Functor>
int fsolve(matrix<Real> &x, const Functor &F, matrix<Real> (Functor::*ptr)(matrix<Real>), Real precision, int maxIterations) {
..
while () {
matrix<Real> fx = (F.*ptr)(x); // delegate
}
..
}
..
BakerValve instance;
fsolve<double, BakerValve>(x1, instance, &BakerValve::constraintF1);
..
fsolve<double, BakerValve>(x2, instance, &BakerValve::constraintF2);
Some additional features of my solver are the Jacobian can either be estimated by a default solver or computed analytically by a user-supplied function.
The implementation seems to work well and always converge. I haven't implemented an aesthetically pleasing visualization of its output, as this was mostly scratching an itch for implementing the constraint solver.
No comments:
Post a Comment