diff options
author | Thomas Capricelli <orzel@freehackers.org> | 2009-09-28 02:55:30 +0200 |
---|---|---|
committer | Thomas Capricelli <orzel@freehackers.org> | 2009-09-28 02:55:30 +0200 |
commit | 87be19de4aec6f0f04ff3a5f0304ca999cecbd13 (patch) | |
tree | 10da9b37e1d17f19b3225a1bb5ceab952e372fb8 /unsupported/Eigen/src/NumericalDiff | |
parent | 206b5e39723c51eaa9550db04d459ebd07604415 (diff) |
central sheme for numerical diff
Diffstat (limited to 'unsupported/Eigen/src/NumericalDiff')
-rw-r--r-- | unsupported/Eigen/src/NumericalDiff/NumericalDiff.h | 51 |
1 files changed, 40 insertions, 11 deletions
diff --git a/unsupported/Eigen/src/NumericalDiff/NumericalDiff.h b/unsupported/Eigen/src/NumericalDiff/NumericalDiff.h index 223cf9e0f..276b315f8 100644 --- a/unsupported/Eigen/src/NumericalDiff/NumericalDiff.h +++ b/unsupported/Eigen/src/NumericalDiff/NumericalDiff.h @@ -28,7 +28,13 @@ namespace Eigen { -template<typename Functor> class NumericalDiff : public Functor +enum NumericalDiffMode { + Forward, + Central +}; + + +template<typename Functor, NumericalDiffMode mode=Forward> class NumericalDiff : public Functor { public: typedef typename Functor::Scalar Scalar; @@ -62,14 +68,23 @@ public: int nfev=0; const int n = _x.size(); const Scalar eps = ei_sqrt((std::max(epsfcn,epsilon<Scalar>() ))); - ValueType val, fx; + ValueType val1, val2; InputType x = _x; // TODO : we should do this only if the size is not already known - val.resize(Functor::values()); - fx.resize(Functor::values()); + val1.resize(Functor::values()); + val2.resize(Functor::values()); - // compute f(x) - Functor::operator()(x, fx); + switch(mode) { + case Forward: + // compute f(x) + Functor::operator()(x, val1); nfev++; + break; + case Central: + // do nothing + break; + default: + assert(false); + }; /* Function Body */ @@ -78,11 +93,25 @@ public: if (h == 0.) { h = eps; } - x[j] += h; - Functor::operator()(x, val); - nfev++; - x[j] = _x[j]; - jac.col(j) = (val-fx)/h; + switch(mode) { + case Forward: + x[j] += h; + Functor::operator()(x, val2); + nfev++; + x[j] = _x[j]; + jac.col(j) = (val2-val1)/h; + break; + case Central: + x[j] += h; + Functor::operator()(x, val2); nfev++; + x[j] -= 2*h; + Functor::operator()(x, val1); nfev++; + x[j] = _x[j]; + jac.col(j) = (val2-val1)/(2*h); + break; + default: + assert(false); + }; } return nfev; } |