33 #include "autodiff/forward/dual.hpp"
34 #include "autodiff/forward/real.hpp"
36 #include <autodiff/forward/dual/dual.hpp>
43 template <
typename T >
55 double exp(
double x );
63 constexpr
double radToDeg(
const double alpha )
70 constexpr
double degToRad(
const double alpha )
79 return scalar >= 0 ? scalar : 0.0;
86 return scalar >= 0 ? 1 : 0;
91 return scalar > 0 ? 1 : 0;
96 template <
typename T >
97 constexpr
int sgn(
const T& val ) noexcept
99 return ( T( 0 ) < val ) - ( val < T( 0 ) );
102 double makeReal(
const double& value );
103 double makeReal(
const std::complex< double >& value );
104 double makeReal(
const autodiff::real& value );
106 template <
typename T,
typename G >
107 double makeReal(
const autodiff::detail::Dual< T, G >& number )
109 return double( number );
112 template <
typename T,
int... Rest >
113 Eigen::Matrix< double, Rest... >
makeReal(
const Eigen::Matrix< T, Rest... > mat )
115 Eigen::Matrix< double, Rest... > out;
116 const size_t m =
static_cast< size_t >( mat.rows() );
117 const size_t n =
static_cast< size_t >( mat.cols() );
119 for (
size_t i = 0; i < m; i++ ) {
120 for (
size_t j = 0; j < n; j++ ) {
121 out( i, j ) =
makeReal( mat( i, j ) );
126 template <
typename T >
127 Eigen::VectorXd
makeReal( Eigen::Vector< T, Eigen::Dynamic > in )
130 int inSize = in.size();
131 Eigen::VectorXd out( inSize );
132 for (
int i = 0; i < inSize; i++ ) {
133 out( i ) = double( in( i ) );
141 template <
int nRows,
int nCols >
142 Eigen::Matrix< double, nRows, nCols >
macaulyMatrix(
const Eigen::Matrix< double, nRows, nCols >& mat )
144 Eigen::Matrix< double, nRows, nCols > positivePart = Eigen::Matrix< double, nRows, nCols >::Zero();
145 for (
int i = 0; i < nRows; i++ ) {
146 for (
int j = 0; j < nCols; j++ ) {
147 positivePart( i, j ) =
macauly( mat( i, j ) );
155 template <
typename functionType,
typename yType,
typename... Args >
156 yType
explicitEuler( yType yN,
const double dt, functionType fRate, Args&&... fRateArgs )
158 return yN + fRate( yN, fRateArgs... ) * dt;
166 template <
int ySize,
typename functionType,
typename... Args >
170 Args&&... fRateArgs )
172 Eigen::Matrix< double, ySize, ySize > fS = Eigen::Matrix< double, ySize, ySize >::Zero();
173 Eigen::Matrix< double, ySize, ySize > Iy = Eigen::Matrix< double, ySize, ySize >::Identity();
175 Eigen::VectorXd leftX( ySize );
176 Eigen::VectorXd rightX( ySize );
178 for (
size_t i = 0; i < ySize; i++ ) {
184 fS.col( i ) = 1. / ( 2. * h ) * ( fRate( rightX, fRateArgs... ) - fRate( leftX, fRateArgs... ) );
187 return yN + ( Iy - dt * fS ).inverse() * dt * fRate( yN, fRateArgs... );
193 template <
int nRows,
int nCols,
typename functionType >
194 Eigen::Matrix< double, nRows, nCols >
centralDiff( functionType f,
const Eigen::Matrix< double, nCols, 1 >& X )
196 Eigen::Matrix< double, nRows, nCols > dXdY = Eigen::Matrix< double, nRows, nCols >::Zero();
198 Eigen::VectorXd leftX( nCols );
199 Eigen::VectorXd rightX( nCols );
201 for (
size_t i = 0; i < nCols; i++ ) {
207 dXdY.col( i ) = 1. / ( 2. * h ) * ( f( rightX ) - f( leftX ) );
217 template <
typename functionType,
typename yType,
typename... Args >
220 yType fN = fRate( yN, fRateArgs... );
221 yType u = yN + fN * dt;
222 yType v = yN + fN * dt / 2.;
223 yType w = v + fRate( v, fRateArgs... ) * dt / 2.;
232 template <
int ySize,
typename functionType,
typename... Args >
234 Eigen::Matrix< double, ySize, 1 > yN,
238 Args&&... fRateArgs )
241 typedef Eigen::Matrix< double, ySize, 1 > ySized;
242 ySized fN = fRate( yN, fRateArgs... );
243 ySized u = yN + fN * dt;
244 ySized v = yN + fN * dt / 2.;
245 ySized w = v + fRate( v, fRateArgs... ) * dt / 2.;
246 ySized yNew = 2. * w - u;
249 const double AERR = 1.0;
250 const double aI = AERR / TOL;
251 const double rI = 1.0;
253 ySized ESTVec = ySized::Zero();
255 for (
int i = 0; i < ySize; i++ ) {
256 scaling = aI + rI * std::max( abs( yNew( i ) ), abs( yN( i ) ) );
257 ESTVec( i ) = abs( w( i ) - u( i ) ) / abs( scaling );
260 const double EST = ESTVec.maxCoeff();
261 const double tauNew = dt * std::min( 2., std::max( 0.2, 0.9 * std::sqrt( TOL / EST ) ) );
263 return std::make_tuple( yNew, tauNew );