stepsize.hpp

Go to the documentation of this file.

#ifndef INCLUDE_RICCATI_STEPSIZE_HPP
#define INCLUDE_RICCATI_STEPSIZE_HPP
 
#include <riccati/utils.hpp>
 
namespace riccati {
template <typename SolverInfo, typename FloatingPoint>
inline FloatingPoint choose_nonosc_stepsize(SolverInfo&& info, FloatingPoint x0,
                                            FloatingPoint h,
                                            FloatingPoint epsilon_h) {
  auto ws = omega(info, riccati::scale(info.xp(), x0, h));
  if (ws.maxCoeff() > (1.0 + epsilon_h) / std::abs(h)) {
    return choose_nonosc_stepsize(info, x0, static_cast<FloatingPoint>(h / 2.0),
                                  epsilon_h);
  } else {
    return h;
  }
}
 
template <typename SolverInfo, typename FloatingPoint>
inline auto choose_osc_stepsize(SolverInfo&& info, FloatingPoint x0,
                                FloatingPoint h, FloatingPoint epsilon_h) {
  auto t = eval(info.alloc_, riccati::scale(info.xp_interp(), x0, h));
  auto s = eval(info.alloc_, riccati::scale(info.xp(), x0, h));
  // TODO: Use a memory arena for these
  auto ws = omega(info, s).eval();
  auto gs = gamma(info, s).eval();
  auto omega_analytic = eval(info.alloc_, omega(info, t));
  auto omega_estimate = info.L() * ws;
  auto gamma_analytic = eval(info.alloc_, gamma(info, t));
  auto gamma_estimate = info.L() * gs;
  FloatingPoint max_omega_err
      = (((omega_estimate - omega_analytic).array() / omega_analytic.array())
             .abs())
            .maxCoeff();
  FloatingPoint max_gamma_err
      = (((gamma_estimate - gamma_analytic).array() / gamma_analytic.array())
             .abs())
            .maxCoeff();
  FloatingPoint max_err = std::max(max_omega_err, max_gamma_err);
  if (max_err <= epsilon_h) {
    if (info.p_ != info.n_) {
      auto xn_scaled = eval(info.alloc_, riccati::scale(info.xn(), x0, h));
      ws = omega(info, xn_scaled);
      gs = gamma(info, xn_scaled);
    }
    return std::make_tuple(h, ws, gs);
  } else {
    auto h_scaling = static_cast<FloatingPoint>(std::min(
        0.7, 0.9 * std::pow(epsilon_h / max_err, (1.0 / (info.p_ - 1.0)))));
    return choose_osc_stepsize(info, x0, h * h_scaling, epsilon_h);
  }
}
 
}  // namespace riccati
#endif