initial_step_size_calculator.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 <algorithm>
#include <cmath>  // IWYU pragma: keep
#include <tuple>
 
#include "num_collect/logging/log_tag_view.h"
#include "num_collect/logging/logging_macros.h"
#include "num_collect/logging/logging_mixin.h"
#include "num_collect/ode/concepts/formula.h"
#include "num_collect/ode/error_tolerances.h"
#include "num_collect/ode/evaluation_type.h"
#include "num_collect/ode/impl/get_least_known_order.h"
#include "num_collect/ode/step_size_limits.h"
 
namespace num_collect::ode {
 
constexpr auto initial_step_size_calculator_log_tag =
    logging::log_tag_view("num_collect::ode::initial_step_size_calculator");
 
template <concepts::formula Formula>
class initial_step_size_calculator : public logging::logging_mixin {
public:
    using formula_type = Formula;
 
    using problem_type = typename formula_type::problem_type;
 
    using variable_type = typename problem_type::variable_type;
 
    using scalar_type = typename problem_type::scalar_type;
 
    initial_step_size_calculator()
        : logging::logging_mixin(initial_step_size_calculator_log_tag) {}
 
    [[nodiscard]] auto calculate(problem_type& problem,
        const scalar_type& initial_time, const variable_type& initial_variable,
        const step_size_limits<scalar_type>& limits,
        const error_tolerances<variable_type>& tolerances) const
        -> scalar_type {
        problem.evaluate_on(initial_time, initial_variable,
            evaluation_type{.diff_coeff = true});
        const variable_type initial_diff = problem.diff_coeff();
 
        const auto [step_size_from_diff, initial_diff_norm] =
            calculate_step_size_from_diff(
                initial_variable, initial_diff, tolerances);
 
        const scalar_type step_size_from_second_diff =
            calculate_step_size_from_second_diff(problem, initial_time,
                initial_variable, initial_diff, step_size_from_diff,
                initial_diff_norm, tolerances);
 
        const scalar_type step_size_without_limit =
            std::min(1e+2 * step_size_from_diff, step_size_from_second_diff);
        NUM_COLLECT_LOG_TRACE(this->logger(),
            "Selection of step size without limits: {}",
            step_size_without_limit);
 
        const scalar_type final_step_size =
            limits.apply(step_size_without_limit);
        NUM_COLLECT_LOG_TRACE(this->logger(),
            "Final selection of step size: {}", step_size_without_limit);
 
        return final_step_size;
    }
 
private:
    [[nodiscard]] auto calculate_step_size_from_diff(
        const variable_type& initial_variable,
        const variable_type& initial_diff,
        const error_tolerances<variable_type>& tolerances) const
        -> std::tuple<scalar_type, scalar_type> {
        // d0 in Hairer1993
        const scalar_type initial_variable_norm =
            tolerances.calc_norm(initial_variable, initial_variable);
        NUM_COLLECT_LOG_TRACE(
            this->logger(), "Norm of variable: {}", initial_variable_norm);
 
        // d1 in Hairer1993
        const scalar_type initial_diff_norm =
            tolerances.calc_norm(initial_variable, initial_diff);
        NUM_COLLECT_LOG_TRACE(
            this->logger(), "Norm of first derivative: {}", initial_diff_norm);
 
        // h0 in Hairer1993
        const scalar_type step_size_from_diff =
            (initial_variable_norm >= static_cast<scalar_type>(1e-5) &&
                initial_diff_norm >= static_cast<scalar_type>(1e-5))
            ? static_cast<scalar_type>(1e-2) * initial_variable_norm /
                initial_diff_norm
            : static_cast<scalar_type>(1e-6);
        NUM_COLLECT_LOG_TRACE(this->logger(),
            "First estimate of step size using differential coefficient: {}",
            step_size_from_diff);
 
        return {step_size_from_diff, initial_diff_norm};
    }
 
    [[nodiscard]] auto calculate_step_size_from_second_diff(
        problem_type& problem, const scalar_type& initial_time,
        const variable_type& initial_variable,
        const variable_type& initial_diff,
        const scalar_type& step_size_from_diff,
        const scalar_type& initial_diff_norm,
        const error_tolerances<variable_type>& tolerances) const {
        // y1 in Hairer1993 (Explicit Euler method)
        const variable_type euler_updated_variable =
            initial_variable + step_size_from_diff * initial_diff;
 
        problem.evaluate_on(initial_time + step_size_from_diff,
            euler_updated_variable, evaluation_type{.diff_coeff = true});
        const variable_type& euler_updated_diff = problem.diff_coeff();
 
        // d2 in Hairer1993 (Evaluation of second derivative)
        const scalar_type second_diff_norm =
            tolerances.calc_norm(
                initial_variable, euler_updated_diff - initial_diff) /
            step_size_from_diff;
        NUM_COLLECT_LOG_TRACE(
            this->logger(), "Norm of second derivative: {}", second_diff_norm);
 
        const scalar_type larger_norm =
            std::max(initial_diff_norm, second_diff_norm);
        constexpr scalar_type exponent_of_order = static_cast<scalar_type>(1) /
            static_cast<scalar_type>(
                impl::get_least_known_order<Formula>() + 1);
        // h1 in Hairer1993
        const scalar_type step_size_from_second_diff =
            (larger_norm > static_cast<scalar_type>(1e-15))
            ? std::pow(static_cast<scalar_type>(1e-2) / larger_norm,
                  exponent_of_order)
            : std::max(static_cast<scalar_type>(1e-6),
                  static_cast<scalar_type>(1e-3) * step_size_from_diff);
        NUM_COLLECT_LOG_TRACE(this->logger(),
            "Second estimate of step size using second derivative: {}",
            step_size_from_second_diff);
 
        return step_size_from_second_diff;
    }
};
 
}  // namespace num_collect::ode