GaussLegendre6.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 *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#pragma once
 
#ifndef SANDALS_GAUSSLEGENDRE6_HH
#define SANDALS_GAUSSLEGENDRE6_HH
 
#include <Sandals.hh>
#include <Sandals/RungeKutta.hh>
 
namespace Sandals
{
  /*
  Symplectic Runge-Kutta Methods of High Order Based on W-Transformation, Kaifeng Xia, Yuhao Cong
  and Geng Sun. Journal of Applied Analysis ans Computation, Volume 7, Number 3, 2017(8), 1185-1199
  http://www.jaac-online.com/data/article/jaac/preview/pdf/20170325.pdf
  */

  template <typename Real>
  class GaussLegendre6Tableau : public Tableau<Real, 3>
  {
  public:
    using typename Tableau<Real, 3>::Type;
    using typename Tableau<Real, 3>::Vector;
    using typename Tableau<Real, 3>::Matrix;

    GaussLegendre6Tableau<Real>() {
      this->name  = "GaussLegendre6";
      this->type  = Type::IRK;
      this->order = 6;
      Real t_1{std::sqrt(15.0)/10.0};
      Real t_2{std::sqrt(15.0)/15.0};
      Real t_3{std::sqrt(15.0)/24.0};
      Real t_4{std::sqrt(15.0)/30.0};
      Real w{5.0/36.0};
      Real z{2.0/9.0};
      this->A << w,           z-t_2,   w-t_4,
                 w+t_3,       z,       w-t_3,
                 w+t_4,       z+t_2,   w;
      this->b << 5.0/18.0,    4.0/9.0, 5.0/18.0;
      this->c << 1.0/2.0-t_1, 1.0/2.0, 1.0/2.0+t_1;
    }
  }; // class GaussLegendre6Tableau

  template <typename Real, Integer N, Integer M = 0>
  class GaussLegendre6 : public RungeKutta<Real, 3, N, M>
  {
  public:
    using System = typename Implicit<Real, N, M>::Pointer; 

    GaussLegendre6() : RungeKutta<Real, 3, N, M>(GaussLegendre6Tableau<Real>()) {}

    GaussLegendre6(System t_system) : RungeKutta<Real, 3, N, M>(GaussLegendre6Tableau<Real>(), t_system) {}
 
  }; // class GaussLegendre6
 
} // namespace Sandals
 
#endif // SANDALS_GAUSSLEGENDRE6_HH