BVP.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 * Copyright (c) 2025, Davide Stocco and Enrico Bertolazzi.                                      *
 *                                                                                               *
 * The Sandals project is distributed under the BSD 2-Clause License.                            *
 *                                                                                               *
 * Davide Stocco                                                               Enrico Bertolazzi *
 * University of Trento                                                     University of Trento *
 * e-mail: davide.stocco@unitn.it                             e-mail: enrico.bertolazzi@unitn.it *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef SANDALS_BVP_HH
#define SANDALS_BVP_HH
 
#include <Sandals/RungeKutta.hh>
 
namespace Sandals
{
 
  /*\
   |   ______     ______
   |  | __ ) \   / /  _ \
   |  |  _ \\ \ / /| |_) |
   |  | |_) |\ V / |  __/
   |  |____/  \_/  |_|
   |
  \*/

  template <typename Real, Integer N, Integer M = 0>
  class BVP
  {
  public:
    using VectorF = typename Implicit<Real, N, M>::VectorF; 
    using Pointer = std::shared_ptr<Implicit<Real, N, M>>; 
    using NewtonB = Optimist::RootFinder::Newton<Real, N>; 
 
  private:
    mutable NewtonB m_newton; 
 
    RungeKutta<Real, N, M> m_rk; 
    std::string m_name; 
 
  public:
    BVP(std::string t_name) : m_name(t_name) {}

    RungeKutta<Real, N, M> & rk() {return this->m_rk;}

    RungeKutta<Real, N, M> const & rk() const {return this->m_rk;}

    virtual VectorF b(VectorF const & x_ini, VectorF const & x_end) const = 0;

    virtual MatrixJF Jb_x_ini(VectorF const & x_ini, VectorF const & x_end) const = 0;

    virtual MatrixJF Jb_x_end(VectorF const & x_ini, VectorF const & x_end) const = 0;

    bool solve_single_shooting(VectorX const & t_mesh, VectorN const & ics, Solution<Real, N, M> & sol) const
    {
      using Eigen::last;
 
      #define CMD "Sandals::BVP::solve_single_shooting(...): "
 
      // Create the lambda function that evaluates the boundary conditions after integration
      auto b = [this](VectorF const & x_ini, VectorF & b_fun) -> VectorF {
        // Integrate the system using the Runge-Kutta method
        if (!this->m_rk.solve(t_mesh, x_ini, sol)) {
          SANDALS_ERROR(CMD "failed to integrate the system " << this->m_name << ".");
          return VectorF::Zero(N);
        }
        return this->b(x_ini, sol.x.col(last));
      };
 
      // Istantiate a Newton solver for the boundary conditions
      return this->m_newton.solve(b, Jb, x_ini, x_sol);
 
      #undef CMD
    }
 
  }; // class BVP
 
} // namespace Sandals
 
#endif // SANDALS_BVP_HH