avf4_formula.h

Go to the documentation of this file.

/*
 * Copyright 2021 MusicScience37 (Kenta Kabashima)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#pragma once
 
#include "num_collect/base/exception.h"
#include "num_collect/base/index_type.h"
#include "num_collect/base/norm.h"
#include "num_collect/base/precondition.h"
#include "num_collect/constants/half.h"  // IWYU pragma: keep
#include "num_collect/constants/one.h"   // IWYU pragma: keep
#include "num_collect/constants/zero.h"  // IWYU pragma: keep
#include "num_collect/integration/gauss_legendre_integrator.h"
#include "num_collect/logging/log_tag_view.h"
#include "num_collect/logging/logging_macros.h"
#include "num_collect/ode/avf/impl/avf_integrand.h"
#include "num_collect/ode/concepts/differentiable_problem.h"
#include "num_collect/ode/evaluation_type.h"
#include "num_collect/ode/non_embedded_formula_wrapper.h"
#include "num_collect/ode/simple_solver.h"
 
namespace num_collect::ode::avf {

template <concepts::differentiable_problem Problem>
class avf4_formula {
public:
    using problem_type = Problem;

    using variable_type = typename problem_type::variable_type;

    using scalar_type = typename problem_type::scalar_type;

    using jacobian_type = typename problem_type::jacobian_type;

    static constexpr index_type stages = 1;

    static constexpr index_type order = 4;

    static constexpr auto log_tag =
        logging::log_tag_view("num_collect::ode::avf::avf4_formula");

    explicit avf4_formula(const problem_type& problem = problem_type())
        : integrand_(problem) {}

    void step(scalar_type time, scalar_type step_size,
        const variable_type& current, variable_type& estimate) {
        const auto dim = current.size();
 
        problem().evaluate_on(time, current,
            evaluation_type{.diff_coeff = true, .jacobian = true});
        static constexpr scalar_type coeff_jacobi =
            -static_cast<scalar_type>(1) / static_cast<scalar_type>(12);
        jacobian_type coeff = step_size *
            (jacobian_type::Identity(dim, dim) +
                coeff_jacobi * step_size * step_size * problem().jacobian() *
                    problem().jacobian());
 
        estimate = current + coeff * problem().diff_coeff();
 
        integrand_.time(time);
        integrand_.prev_var(current);
        variable_type prev_estimate;
        constexpr index_type max_loops = 10000;
        for (index_type i = 0; i < max_loops; ++i) {
            problem().evaluate_on(time,
                constants::half<scalar_type> * (current + estimate),
                evaluation_type{.diff_coeff = true, .jacobian = true});
            coeff = step_size *
                (jacobian_type::Identity(dim, dim) +
                    coeff_jacobi * step_size * step_size *
                        problem().jacobian() * problem().jacobian());
 
            integrand_.next_var(estimate);
            prev_estimate = estimate;
            estimate = current +
                coeff *
                    integrator_(integrand_, constants::zero<scalar_type>,
                        constants::one<scalar_type>);
            if (norm(estimate - prev_estimate) < tol_residual_norm_) {
                return;
            }
        }
    }

    [[nodiscard]] auto problem() -> problem_type& {
        return integrand_.problem();
    }

    [[nodiscard]] auto problem() const -> const problem_type& {
        return integrand_.problem();
    }

    void tol_residual_norm(scalar_type val) {
        NUM_COLLECT_PRECONDITION(val > static_cast<scalar_type>(0),
            "Tolerance of residual norm must be a positive value.");
        tol_residual_norm_ = val;
    }
 
private:
    impl::avf_integrand<problem_type> integrand_;

    static constexpr index_type integrator_degree = 5;

    integration::gauss_legendre_integrator<variable_type(scalar_type)>
        integrator_{integrator_degree};

    static constexpr auto default_tol_residual_norm =
        static_cast<scalar_type>(1e-8);

    scalar_type tol_residual_norm_{default_tol_residual_norm};
};

template <concepts::differentiable_problem Problem>
using avf4_solver = simple_solver<avf4_formula<Problem>>;

template <concepts::differentiable_problem Problem>
using avf4_auto_solver = non_embedded_auto_solver<avf4_formula<Problem>>;
 
}  // namespace num_collect::ode::avf