src/tools/runge_kutta.h

Utility to execute Runge-Kutta integration using Eigen. More…

Namespaces

Name
allpix
Helper class to hold support layers for a detector model.
allpix::tableau

Classes

Name
class allpix::RungeKutta
Class to perform arbitrary Runge-Kutta integration.
class allpix::RungeKutta::Step
Utility type to return both the value and the error at every step.

Functions

Name
const auto RK3((Eigen::Matrix< double, 5, 3 >()« 0, 0, 0, 1.0/2, 0, 0, -1, 2, 0, 1.0/6, 2.0/3, 1.0/6, 0, 0, 0).finished() )
Kutta’s third order method.
const auto RK4((Eigen::Matrix< double, 6, 4 >()« 0, 0, 0, 0, 1.0/2, 0, 0, 0, 0, 1.0/2, 0, 0, 0, 0, 1, 0, 1.0/6, 1.0/3, 1.0/3, 1.0/6, 0, 0, 0, 0).finished() )
Classic original Runge-Kutta method.
const auto RK5((Eigen::Matrix< double, 8, 6 >()« 0, 0, 0, 0, 0, 0, 1.0/4, 0, 0, 0, 0, 0, 3.0/32, 9.0/32, 0, 0, 0, 0, 1932.0/2197, -7200.0/2197, 7296.0/2197, 0, 0, 0, 439.0/216, -8, 3680.0/513, -845.0/4104, 0, 0, -8.0/27, 2, -3544.0/2565, 1859.0/4104, -11.0/40, 0, 16.0/135, 0, 6656.0/12825, 28561.0/56430, -9.0/50, 2.0/55, 25.0/216, 0, 1408.0/2565, 2197.0/4104, -1.0/5, 0).finished() )
Runge-Kutta-Fehlberg method Values from https://ntrs.nasa.gov/citations/19680027281, p.13, Table III.
const auto RKCK((Eigen::Matrix< double, 8, 6 >()« 0, 0, 0, 0, 0, 0, 1.0/5, 0, 0, 0, 0, 0, 3.0/40, 9.0/40, 0, 0, 0, 0, 3.0/10, -9.0/10, 6.0/5, 0, 0, 0, -11.0/54, 5.0/2, -70.0/27, 35.0/27, 0, 0, 1631.0/55296, 175.0/512, 575.0/13824, 44275.0/110592, 253.0/4096, 0, 37.0/378, 0, 250.0/621, 125.0/594, 0, 512.0/1771, 2825.0/27648, 0, 18575.0/48384, 13525.0/55296, 277.0/14336, 1.0/4).finished() )
Runge-Kutta-Cash-Karp method.
template <typename T ,int S,int D =3,class… Args>
RungeKutta< T, S, D >		 make_runge_kutta(const Eigen::Matrix< T, S+2, S > & tableau, Args &&… args)
Utility function to create RungeKutta class using template deduction.

Detailed Description

Utility to execute Runge-Kutta integration using Eigen.

Copyright: Copyright (c) 2017-2024 CERN and the Allpix Squared authors. This software is distributed under the terms of the MIT License, copied verbatim in the file “LICENSE.md”. In applying this license, CERN does not waive the privileges and immunities granted to it by virtue of its status as an Intergovernmental Organization or submit itself to any jurisdiction. SPDX-License-Identifier: MIT

Functions Documentation

function RK3

static const auto RK3(
    (Eigen::Matrix< double, 5, 3 >()<< 0, 0, 0, 1.0/2, 0, 0, -1, 2, 0, 1.0/6, 2.0/3, 1.0/6, 0, 0, 0).finished() 
)

Kutta’s third order method.

Warning: Without error function

function RK4

static const auto RK4(
    (Eigen::Matrix< double, 6, 4 >()<< 0, 0, 0, 0, 1.0/2, 0, 0, 0, 0, 1.0/2, 0, 0, 0, 0, 1, 0, 1.0/6, 1.0/3, 1.0/3, 1.0/6, 0, 0, 0, 0).finished() 
)

Classic original Runge-Kutta method.

Warning: Without error function

function RK5

static const auto RK5(
    (Eigen::Matrix< double, 8, 6 >()<< 0, 0, 0, 0, 0, 0, 1.0/4, 0, 0, 0, 0, 0, 3.0/32, 9.0/32, 0, 0, 0, 0, 1932.0/2197, -7200.0/2197, 7296.0/2197, 0, 0, 0, 439.0/216, -8, 3680.0/513, -845.0/4104, 0, 0, -8.0/27, 2, -3544.0/2565, 1859.0/4104, -11.0/40, 0, 16.0/135, 0, 6656.0/12825, 28561.0/56430, -9.0/50, 2.0/55, 25.0/216, 0, 1408.0/2565, 2197.0/4104, -1.0/5, 0).finished() 
)

Runge-Kutta-Fehlberg method Values from https://ntrs.nasa.gov/citations/19680027281, p.13, Table III.

function RKCK

static const auto RKCK(
    (Eigen::Matrix< double, 8, 6 >()<< 0, 0, 0, 0, 0, 0, 1.0/5, 0, 0, 0, 0, 0, 3.0/40, 9.0/40, 0, 0, 0, 0, 3.0/10, -9.0/10, 6.0/5, 0, 0, 0, -11.0/54, 5.0/2, -70.0/27, 35.0/27, 0, 0, 1631.0/55296, 175.0/512, 575.0/13824, 44275.0/110592, 253.0/4096, 0, 37.0/378, 0, 250.0/621, 125.0/594, 0, 512.0/1771, 2825.0/27648, 0, 18575.0/48384, 13525.0/55296, 277.0/14336, 1.0/4).finished() 
)

Runge-Kutta-Cash-Karp method.

function make_runge_kutta

template <typename T ,
int S,
int D =3,
class... Args>
RungeKutta< T, S, D > make_runge_kutta(
    const Eigen::Matrix< T, S+2, S > & tableau,
    Args &&... args
)

Utility function to create RungeKutta class using template deduction.

Parameters:

Return: Instantiation of RungeKutta class with the forwarded arguments

Source code


#ifndef ALLPIX_RUNGE_KUTTA_H
#define ALLPIX_RUNGE_KUTTA_H

#include <functional>

#include <Eigen/Core>
#include <Eigen/Geometry>

namespace allpix {

    template <typename T, int S, int D = 3> class RungeKutta {
    public:
        // FIXME Is this a appropriate return type or should pairs be preferred
        class Step {
        public:
            Eigen::Matrix<T, D, 1> value;
            Eigen::Matrix<T, D, 1> error;
        };

        using StepFunction = std::function<Eigen::Matrix<T, D, 1>(T, Eigen::Matrix<T, D, 1>)>;

        RungeKutta(Eigen::Matrix<T, S + 2, S> tableau,
                   StepFunction function,
                   T step_size,
                   Eigen::Matrix<T, D, 1> initial_y,
                   T initial_t = 0)
            : tableau_(std::move(tableau)), function_(std::move(function)), h_(std::move(step_size)),
              y_(std::move(initial_y)), t_(std::move(initial_t)) {
            error_.setZero();
        }

        void setTimeStep(T step_size) { h_ = std::move(step_size); }
        T getTimeStep() const { return h_; }

        void setValue(Eigen::Matrix<T, D, 1> y) { y_ = std::move(y); }

        Eigen::Matrix<T, D, 1> getValue() const { return y_; }
        Eigen::Matrix<T, D, 1> getError() const { return error_; }
        T getTime() const { return t_; }
        void advanceTime(double t) { t_ += t; }

        Step step() {
            // Initialize values
            Step step;
            Eigen::Matrix<T, D, 1> ys;
            Eigen::Matrix<T, D, 1> yse;
            ys.setZero();
            yse.setZero();

            // Compute step
            Eigen::Matrix<T, S, D> k;
            for(int i = 0; i < S; ++i) {
                Eigen::Matrix<T, D, 1> yt = y_;
                T tt = t_;
                for(int j = 0; j < i; ++j) {
                    yt += h_ * tableau_(i, j) * k.row(j);
                    tt += tableau_(i, j);
                }
                k.row(i) = function_(tt, yt);

                ys += h_ * tableau_(S, i) * k.row(i);
                yse += h_ * tableau_(S + 1, i) * k.row(i);
            }

            // Update values with new step
            y_ += ys;
            t_ += h_;
            error_ += ys - yse;

            // Return step information
            step.value = ys;
            step.error = ys - yse;
            return step;
        }

        Step step(int amount) {
            Step result;
            result.value.setZero();
            result.error.setZero();
            for(int i = 0; i < amount; ++i) {
                Step single = step();
                result.value += single.value;
                result.error += single.error;
            }
            return result;
        }

    private:
        const Eigen::Matrix<T, S + 2, S> tableau_;
        StepFunction function_;
        // Step size
        T h_;

        // Vector to integrate
        Eigen::Matrix<T, D, 1> y_;
        // Total error vector
        Eigen::Matrix<T, D, 1> error_;
        // Current time
        T t_;
    };

    // clang-format off
    namespace tableau {
        static const auto RK3((Eigen::Matrix<double, 5, 3>() <<
            0, 0, 0,
            1.0/2, 0, 0,
            -1, 2, 0,
            1.0/6, 2.0/3, 1.0/6,
            0, 0, 0).finished());
        static const auto RK4((Eigen::Matrix<double, 6, 4>() <<
            0, 0, 0, 0,
            1.0/2, 0, 0, 0,
            0, 1.0/2, 0, 0,
            0, 0, 1, 0,
            1.0/6, 1.0/3, 1.0/3, 1.0/6,
            0, 0, 0, 0).finished());
        static const auto RK5((Eigen::Matrix<double, 8, 6>() <<
            0, 0, 0, 0, 0, 0,
            1.0/4, 0, 0, 0, 0, 0,
            3.0/32, 9.0/32, 0, 0, 0, 0,
            1932.0/2197, -7200.0/2197, 7296.0/2197, 0, 0, 0,
            439.0/216, -8, 3680.0/513, -845.0/4104, 0, 0,
            -8.0/27, 2, -3544.0/2565, 1859.0/4104, -11.0/40, 0,
            16.0/135, 0, 6656.0/12825, 28561.0/56430, -9.0/50, 2.0/55,
            25.0/216, 0, 1408.0/2565, 2197.0/4104, -1.0/5, 0).finished());
        static const auto RKCK((Eigen::Matrix<double, 8, 6>() <<
            0, 0, 0, 0, 0, 0,
            1.0/5, 0, 0, 0, 0, 0,
            3.0/40, 9.0/40, 0, 0, 0, 0,
            3.0/10, -9.0/10, 6.0/5, 0, 0, 0,
            -11.0/54, 5.0/2, -70.0/27, 35.0/27, 0, 0,
            1631.0/55296, 175.0/512, 575.0/13824, 44275.0/110592, 253.0/4096, 0,
            37.0/378, 0, 250.0/621, 125.0/594, 0, 512.0/1771,
            2825.0/27648, 0, 18575.0/48384, 13525.0/55296, 277.0/14336, 1.0/4).finished());
    }
    // clang-format on

    template <typename T, int S, int D = 3, class... Args>
    RungeKutta<T, S, D> make_runge_kutta(const Eigen::Matrix<T, S + 2, S>& tableau, Args&&... args) {
        return RungeKutta<T, S, D>(tableau, std::forward<Args>(args)...);
    }
} // namespace allpix

#endif /* ALLPIX_RUNGE_KUTTA_H */

Updated on 2025-02-27 at 14:14:46 +0000