integrator.hpp

Go to the documentation of this file.

// -------------------------------------------------------//
//
// SHAMROCK code for hydrodynamics
// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
//
// -------------------------------------------------------//
 
#pragma once

 
#include <functional>
 
namespace shammath {
 
    template<class T, class Lambda>
    inline constexpr T integ_riemann_sum(T start, T end, T step, Lambda &&fct) {
        T acc = {};
 
        for (T x = start; x < end; x += step) {
            acc += fct(x) * step;
        }
        return acc;
    }

    template<class T, class Lambda>
    inline constexpr std::pair<std::vector<T>, std::vector<T>> euler_ode(
        T start, T end, T step, Lambda &&ode, T x0, T u0) {
        std::vector<T> U = {u0};
        std::vector<T> X = {x0};
 
        T u_prev = u0;
        T u      = u0;
        for (T x = x0 + step; x < end; x += step) {
            u = u_prev + ode(u_prev, x) * step;
            X.push_back(x);
            U.push_back(u);
            u_prev = u;
        };
        u_prev = u0;
        std::vector<T> X_backward, U_backward;
        for (T x = x0 - step; x > start; x -= step) {
            u = u_prev - ode(u_prev, x) * step;
            X_backward.push_back(x);
            U_backward.push_back(u);
            u_prev = u;
        }
        X.insert(X.begin(), X_backward.rbegin(), X_backward.rend());
        U.insert(U.begin(), U_backward.rbegin(), U_backward.rend());
        return {X, U};
    }
 
} // namespace shammath