IVPSolverLib – Библиотека решения задачи Коши 2-го порядка

IVPSolverLib – это библиотека для численного решения задачи Коши (initial value problem) для обыкновенных дифференциальных уравнений второго порядка.

Она позволяет решать уравнения вида $y'' = f (x, y, y')$ с начальными условиями $y(x_0) = q_0$, $y'(x_0) = q_1$ на отрезке $[a, b]$.

В библиотеке реализовано несколько классических как явных, так и неявных методов интегрирования, имеется возможность анализировать погрешность решения при различном шаге сетки и визуализировать полученные численные решения (как функцию, так и ее первую производную) для сравнения методов.

Преимуществом библиотеки является высокая производительность. Так, все методы выполнены с использованием контейнеров с выделением исключительно выровненной (по 64 бита) памяти aligned_alloc() и с доступом к ней assume_aligned , для наибольшей производительности. Так же везде, где в этом есть смысл, реализован параллелизм на уровне данных SIMD-векторизация.

Поддерживаемые численные методы

В текущей версии IVPSolverLib реализованы следующие численные методы решения ОДУ:

• Метод Эйлера – простой одношаговый явный метод первого порядка точности.
• Метод Эйлера с пересчётом – метод предиктор-корректор (модифицированный Эйлер) второго порядка точности.
• Метод Рунге–Кутты 2-го порядка – явная двухслойная схема второго порядка.
• Метод Рунге–Кутты 4-го порядка – классическая четырехстадийная схема Рунге–Кутты.
• Явный метод Адамса 3-го порядка – трехшаговая схема Адамса–Башфорта третьего порядка точности.

Каждый из этих методов может быть выбран для решения задачи через единый интерфейс библиотеки. Это позволяет сравнивать их между собой по точности и производительности на различных задачах.

Анализ погрешности численного решения

Одной из возможностей IVPSolverLib является оценка точности методов путем сравнения решений на разных шагах сетки. Пользователь может вычислить решение при шаге h, а затем повторить расчет при шаге 2h (более грубая сетка) и h/2 (более плотная сетка). Библиотека позволит автоматически сравнить полученные результаты и оценить абсолютные и относительные отклонения решения на шаге h относительно решений на шагах 2h и h/2.

Так же есть метод интегрирования с адаптивным подбором шага по заданной точности integrate_adaptive, функция сама подберет необходимый шаг по требуемой точности.

Общая сигнатура методов

Методы можно вызывать через общий хаб, как для адаптивного шага integrate_adaptive, так и для фиксированного integrate_freeze

Основные методы интегрирования

template<std::size_t ExpDim, class Span, class Rhs, class Kernel>
    ReturnBuffer<Span> integrate_adaptive(
                                    const Rhs&  F,
                                    typename Span::Scalar a,
                                    typename Span::Scalar b,
                                    typename Span::Scalar q0,
                                    typename Span::Scalar q1,
                                    typename Span::Scalar eps_tol,
                                    Kernel&& kernel
                                )

Здесь:

• ExpDim - предполагаемое количество узлов для достижение требуемой точности (чем ближе к оптимальному числу, тем быстрее алгоритм)
• Span - предпочитаемый view-контейнер
• Rhs - формат правой части уравнения (например лямбда функция)
• Kernel - ядро интегрирования (EulerStepper, EulerRecountStepper, RK2Stepper, RK4Stepper, AdamsStepper)

template<class Span, class Rhs, class Kernel>
    ReturnBuffer<Span> integrate_freeze(
                                    const Rhs&  F,
                                    typename Span::Scalar a,
                                    typename Span::Scalar b,
                                    typename Span::Scalar q0,
                                    typename Span::Scalar q1,
                                    Kernel&& kernel,
                                    std::size_t grid_size
                                )

Здесь:

• Span - предпочитаемый view-контейнер
• Rhs - формат правой части уравнения (например лямбда функция)
• Kernel - ядро интегрирования (EulerStepper, EulerRecountStepper, RK2Stepper, RK4Stepper, AdamsStepper)

Помимо методов вычислений, библиотека предлагает методы анализа полученных решений

• write_norms_table() Вывод таблицы норм разности решений
• write_data_to_json_file() Запись решения в json файл, для дальнейшей визуализации при помощи matplotlib,seaborn, etc

Пример использования

#include <cmath>
#include "odesolver.hpp"
#include "utils.hpp"
using namespace std;

int main(){
    //ДУ
    auto rhs = [](double x, double y1, double y2){return y2 * cos(x) + y1 * sqrt(x + 1) + pow(x,2) - 1;};
    double q0 = -4.7, q1 = 4.2;
    double a = 0.0, b = 1.0;

    //Требуемая точность
    double tol = 1e-5;

    using KernelTRK2 = ode::RK2Stepper<ali_span<double>, decltype(rhs)>;

    //Создание ядра интегрирования
    KernelTRK2 rk2;
    
    //Результат формата (x, y, y', success)
    auto Result_rk2 = ode::integrate_adaptive<101, ali_span<double>>(rhs, a, b, q0, q1, tol, rk2);
    
    ode::compare_results<ali_span<double>>(rhs, a, b, q0, q1, rk2, opt_N_rk2, "rk2", "norms_table.csv", false);
    
    write_data_to_json_file(
        "data.json", 
        make_pair(Result_rk2.x(), "rk2_x"),
        make_pair(Result_rk2.y1(), "rk2_y1"), 
        make_pair(Result_rk2.y2(), "rk2_y2"),
    );
    
    return 0;
}