2

This question is motivated by the following implementation of the Runge Kutta formula for integration of first order ordinary differential equations (ODE).

template <typename Field, typename Vector>
auto
runge_kutta(const Vector& y,
           std::function<Vector(Field, const Vector&)> rhs,
           Field t, Field delta_t)
{
  Vector k1 = rhs(t, y);
  Vector k2 = rhs(t + delta_t/2, y + delta_t/2*k1);
  Vector k3 = rhs(t + delta_t/2, y + delta_t/2*k2);
  Vector k4 = rhs(t + delta_t,   y + delta_t * k3);
  return y + delta_t/6*(k1 + 2*k2 + 2*k3 + k4);
}

This may be used to integrate functions like the following

double f(double t, double y) { return y; }

But the Runge Kutta code above may also be used to solve higher order ODE or ODE with more than one dependent variable. For this the template parameter Vector has to be a container of numbers (or different kind of vectors) with suitably defined arithmetic operators.

Therefore I am looking for an implementation of a type for Vector with the following properties:

  1. The type is a container of numbers or different kind of vectors.
  2. Suitably overloaded operators operator+, operator-, operator* and operator/.
  3. No using declaration in the calling code necessary to use these operators.
  4. The standard library is reused.

The best, I came up with, is to create the operators for std::array in a new namespace

namespace StackOverflow {

  template<typename Field, typename Element, std::size_t Dim>
  auto
  operator*(Field lhs, const std::array<Element, Dim>& rhs);

  // etc., etc.
}

and have suitable using declaration in the calling code (violating point 2 above)

using StackOverflow::operator+; ...

Another possibility is to derive publicly from std::array. This works, but, as I have learned, this is explicitly forbidden by the C++ language standard.

Is it possible to reuse std::array - retaining the first four properties above - without rewriting std::array completely?

user0809348012
  • 21
  • 2
  • Maybe [std::valarray](https://en.cppreference.com/w/cpp/numeric/valarray) is what you're looking for? – Lukas-T Aug 02 '20 at 15:20
  • Well, I considered ```std::valarray```. – user0809348012 Aug 02 '20 at 16:02
  • Advantages of `std::valarray`: It's in the standard library and it (kind of) works. Disadvantages of `std::valarray`: length checks while adding vectors are not at compile time, it's not a container in the STL sense, it's hard to print and it's hard to test. The last two points are due to the fact, that for optimization purposes the return types of the `std::valarray` operators are rather *fluid*. So yes, `std::valarray` might be an answer, but I was hoping for something better. – user0809348012 Aug 02 '20 at 16:14
  • `boost::eigen` exists as a quasi-standard implementation of OO linear algebra, examples of ODE solvers using the `Eigen` library also exist. – Lutz Lehmann Aug 03 '20 at 10:33
  • Looks like there is no solution satisfying my constraints readily available. Usual practice seems to be to create a new type with the required functionality even if that means reimplementing part of standard library. – user0809348012 Aug 15 '20 at 15:16

0 Answers0