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; }