积分器

space-ast2026/1/21大约 4 分钟

积分器

概述

ODE（常微分方程）积分器模块提供了多种数值积分方法，用于求解常微分方程初值问题。该模块参考了 hipparchus-ode 项目的类结构，实现了多种经典的 Runge-Kutta 方法。

核心类结构

ODEIntegrator 基类

ODE 积分器的抽象基类，定义了积分器的通用接口。

class ODEIntegrator {
public:
    virtual ~ODEIntegrator() {};
    
    // 初始化积分器
    virtual errc_t initialize(ODE& ode) = 0;
    
    // 积分ODE
    virtual errc_t integrate(ODE& ode, double t0, double tf, const double* y0, double* yf) = 0;
    
    // 积分ODE一步
    virtual errc_t integrateStep(ODE& ode, double& t, double tf, const double* y0, double* y) = 0;
    
    // 执行一步积分
    virtual errc_t singleStep(ODE& ode, double t0, double step, const double* y0, double* yf) = 0;
};

ODEFixedStepIntegrator

定步长积分器基类，适用于步长固定的积分方法。

ODEVarStepIntegrator

变步长积分器基类，支持自适应步长控制，根据误差估计自动调整步长。

积分器实现

定步长积分器

RK4 - 4阶龙格库塔方法

阶数: 4阶
类型: 定步长
精度: 中等精度，适用于一般精度要求的计算
特点: 经典的4阶Runge-Kutta方法，计算稳定，实现简单

RK4 integrator;
errc_t result = integrator.initialize(ode);
result = integrator.integrate(ode, t0, tf, y0, yf);

RK8 - 8阶龙格库塔方法

阶数: 8阶
类型: 定步长
精度: 高精度，适用于高精度要求的计算
特点: 高阶方法，精度高但计算量较大

RKV8 - 8阶Verner方法

阶数: 8阶
类型: 定步长
精度: 高精度
特点: Verner提出的8阶方法，具有较好的数值稳定性

变步长积分器

RKF45 - Runge-Kutta-Fehlberg 4(5)方法

阶数: 4(5)阶（主阶4阶，误差估计阶5阶）
类型: 变步长
精度: 中等精度，带误差控制
特点: 经典的Fehlberg方法，通过比较4阶和5阶结果来估计误差

RKF45 integrator;
integrator.setTolerance(1e-8);  // 设置容差
integrator.setMinStep(1e-6);    // 设置最小步长
integrator.setMaxStep(1.0);     // 设置最大步长

RKF56 - Runge-Kutta-Fehlberg 5(6)方法

阶数: 5(6)阶
类型: 变步长
精度: 较高精度，带误差控制
特点: 5阶主方法，6阶误差估计

RKF78 - Runge-Kutta-Fehlberg 7(8)方法

阶数: 7(8)阶
类型: 变步长
精度: 高精度，带误差控制
特点: 高阶方法，适用于高精度要求的计算

RKCK - Cash-Karp方法

阶数: 4(5)阶
类型: 变步长
精度: 中等精度，带误差控制
特点: Cash和Karp提出的方法，具有较好的数值特性

使用方法

基本使用流程

定义ODE方程: 继承 ODE 类并实现 evaluate 方法
选择积分器: 根据精度和性能需求选择合适的积分器
配置参数: 设置积分器参数（容差、步长等）
执行积分: 调用积分方法进行求解

使用示例

示例1：基本积分操作

演示如何使用不同的ODE积分器求解简单的常微分方程，包括指数衰减方程和简谐振动方程。

#include "AstMath/ODE.hpp"
#include <iostream>
#include <cmath>
#include <clocale>
#include <string>

AST_USING_NAMESPACE


/// @brief 简单的指数衰减ODE：dy/dt = -y
class ExponentialDecayODE : public ODE {
public:
    errc_t evaluate(const double* y, double* ydot, double t) override {
        ydot[0] = -y[0];  // dy/dt = -y
        return eNoError;
    }
    
    int getDimension() const override {
        return 1;
    }
};

/// @brief 简谐振动ODE：d²x/dt² = -ω²x
class HarmonicOscillatorODE : public ODE {
private:
    double omega_;  // 角频率
    
public:
    HarmonicOscillatorODE(double omega = 1.0) : omega_(omega) {}
    
    errc_t evaluate(const double* y, double* ydot, double t) override {
        // y[0] = x, y[1] = dx/dt
        // dx/dt = v
        // dv/dt = -ω²x
        ydot[0] = y[1];           // dx/dt = v
        ydot[1] = -omega_ * omega_ * y[0];  // dv/dt = -ω²x
        return eNoError;
    }
    
    int getDimension() const override {
        return 2;
    }
};

int main() {
    std::setlocale(LC_ALL, ".UTF-8");
    std::cout << "=== ODE积分器基本使用示例 ===" << std::endl;
    
    // 示例1：指数衰减方程
    std::cout << "\n1. 指数衰减方程 (dy/dt = -y)" << std::endl;
    std::cout << "   解析解: y(t) = y0 * exp(-t)" << std::endl;
    
    ExponentialDecayODE expODE;
    
    // 使用RK4积分器
    RK4 rk4;
    rk4.initialize(expODE);
    
    double t0 = 0.0;
    double tf = 2.0;
    double y0[1] = {1.0};  // y(0) = 1
    double yf_rk4[1]{y0[0]};
    double t = t0;
    errc_t result = rk4.integrate(expODE, yf_rk4, t, tf);
    if (result == eNoError) {
        double exact = y0[0] * std::exp(-tf);
        double error = std::abs(yf_rk4[0] - exact);
        std::cout << "   RK4结果: y(" << tf << ") = " << yf_rk4[0] << std::endl;
        std::cout << "   解析解: y(" << tf << ") = " << exact << std::endl;
        std::cout << "   误差: " << error << std::endl;
    }
    
    // 使用RKF45积分器
    RKF45 rkf45;
    rkf45.setMaxAbsErr(1e-8);
    rkf45.setMaxRelErr(1e-8);
    rkf45.initialize(expODE);
    
    double yf_rkf45[1]{y0[0]};
    t = t0;
    result = rkf45.integrate(expODE, yf_rkf45, t, tf);
    if (result == eNoError) {
        double exact = y0[0] * std::exp(-tf);
        double error = std::abs(yf_rkf45[0] - exact);
        std::cout << "   RKF45结果: y(" << tf << ") = " << yf_rkf45[0] << std::endl;
        std::cout << "   解析解: y(" << tf << ") = " << exact << std::endl;
        std::cout << "   误差: " << error << std::endl;
    }
    
    // 示例2：简谐振动方程
    std::cout << "\n2. 简谐振动方程 (d²x/dt² = -ω²x)" << std::endl;
    std::cout << "   解析解: x(t) = A*cos(ωt + φ)" << std::endl;
    
    double omega = 2.0;  // 角频率
    HarmonicOscillatorODE harmonicODE(omega);
    
    // 初始条件：x(0)=1, v(0)=0
    double y0_harmonic[2] = {1.0, 0.0};
    double yf_harmonic[2]{y0_harmonic[0], y0_harmonic[1]};
    
    // 使用RKF45积分
    rkf45.initialize(harmonicODE);
    t = t0;
    result = rkf45.integrate(harmonicODE, yf_harmonic, t, tf);
    
    if (result == eNoError) {
        double exact_x = std::cos(omega * tf);  // x(t) = cos(ωt)
        double exact_v = -omega * std::sin(omega * tf);  // v(t) = -ωsin(ωt)
        
        std::cout << "   RKF45结果:" << std::endl;
        std::cout << "   位置 x(" << tf << ") = " << yf_harmonic[0] << std::endl;
        std::cout << "   速度 v(" << tf << ") = " << yf_harmonic[1] << std::endl;
        std::cout << "   解析解:" << std::endl;
        std::cout << "   位置 x(" << tf << ") = " << exact_x << std::endl;
        std::cout << "   速度 v(" << tf << ") = " << exact_v << std::endl;
        std::cout << "   位置误差: " << std::abs(yf_harmonic[0] - exact_x) << std::endl;
        std::cout << "   速度误差: " << std::abs(yf_harmonic[1] - exact_v) << std::endl;
    }
    
    std::cout << "\n=== 示例完成 ===" << std::endl;
    
    return 0;
}

示例2：变步长积分器自适应控制

展示变步长积分器的自适应步长控制和误差估计功能，包括不同容差设置的影响。

#include "AstMath/ODE.hpp"
#include "AstGlobal.h"
#include <iostream>
#include <cmath>
#include <vector>
#include <memory>
#include <clocale>
#include <string>

AST_USING_NAMESPACE

/// @brief 刚性方程示例：dy/dt = -1000(y - sin(t)) + cos(t)
class StiffODE : public ODE {
public:
    errc_t evaluate(const double* y, double* ydot, double t) override {
        ydot[0] = -1000.0 * (y[0] - std::sin(t)) + std::cos(t);
        return eNoError;
    }
    
    int getDimension() const override {
        return 1;
    }
};



int main() {
    std::setlocale(LC_ALL, ".UTF-8");
    std::cout << "=== 变步长积分器自适应步长控制示例 ===" << std::endl;
    
    StiffODE stiffODE;
    
    // 测试不同的变步长积分器
    std::vector<std::pair<std::string, ODEVarStepIntegrator*>> integrators = {
        {"RKF45", new RKF45()},
        {"RKF78", new RKF78()},
        {"RKCK", new RKCK()}
    };
    
    for (auto& integrator_pair : integrators) {
        std::string name = integrator_pair.first;
        ODEVarStepIntegrator* integrator = integrator_pair.second;
        
        std::cout << "\n=== " << name << " 积分器 ===" << std::endl;
        
        // 配置积分器参数
        integrator->setMaxAbsErr(1e-14);
        integrator->setMaxRelErr(1e-12);
        integrator->initialize(stiffODE);
        
        // 初始条件
        double t0 = 0.0;
        double tf = 0.5;
        double y0[1] = {2.0};  // y(0) = 2
        
        // 使用积分器进行积分
        double t = t0;
        double yf[1]{y0[0]};
        errc_t result = integrator->integrate(stiffODE, yf, t, tf);
        
        if (result == eNoError) {
            // 获取积分统计信息
            double largest_step = integrator->getLargestStepSize();
            double smallest_step = integrator->getSmallestStepSize();
            int num_steps = integrator->getNumSteps();
            
            // 输出积分结果统计
            std::cout << "积分统计:" << std::endl;
            std::cout << "  总步数: " << num_steps << std::endl;
            std::cout << "  最小步长: " << smallest_step << std::endl;
            std::cout << "  最大步长: " << largest_step << std::endl;
            
            // 输出最终结果
            double exact = std::sin(tf) + (y0[0] - std::sin(t0)) * std::exp(-1000.0 * (tf - t0));
            double error = std::abs(yf[0] - exact);
            
            std::cout << "  数值解: y(" << tf << ") = " << yf[0] << std::endl;
            std::cout << "  解析解: y(" << tf << ") = " << exact << std::endl;
            std::cout << "  误差: " << error << std::endl;
        } else {
            std::cout << "积分失败" << std::endl;
        }
        
        delete integrator;
    }
    
    // 演示不同容差设置的影响
    std::cout << "\n=== 不同容差设置对比 ===" << std::endl;
    
    std::vector<double> tolerances = {1e-4, 1e-6, 1e-8};
    
    for (double tol : tolerances) {
        RKF45 integrator;
        integrator.setMaxAbsErr(tol);
        integrator.setMaxRelErr(tol);
        integrator.initialize(stiffODE);
        
        double t0 = 0.0;
        double tf = 0.5;
        double y0[1] = {2.0};
        double yf[1]{y0[0]};
        double t = t0;
        
        errc_t result = integrator.integrate(stiffODE, yf, t, tf);
        
        if (result == eNoError) {
            double exact = std::sin(tf) + (y0[0] - std::sin(t0)) * std::exp(-1000.0 * (tf - t0));
            double error = std::abs(yf[0] - exact);
            
            std::cout << "容差=" << tol << ": 误差=" << error 
                      << ", 步数统计=" << integrator.getNumSteps() << std::endl;
        }
    }
    
    std::cout << "\n=== 示例完成 ===" << std::endl;
    
    return 0;
}

示例3：性能比较分析

比较不同积分器在精度和计算效率方面的表现，包括定步长和变步长积分器的对比。

#include "AstMath/ODE.hpp"
#include "AstGlobal.h"
#include <iostream>
#include <cmath>
#include <chrono>
#include <vector>
#include <iomanip>
#include <clocale>
#include <string>

AST_USING_NAMESPACE


/// @brief 范德波尔振荡器：d²x/dt² - μ(1-x²)dx/dt + x = 0
class VanDerPolOscillator : public ODE {
private:
    double mu_;  // 非线性参数
    
public:
    VanDerPolOscillator(double mu = 1.0) : mu_(mu) {}
    
    errc_t evaluate(const double* y, double* ydot, double t) override {
        // y[0] = x, y[1] = dx/dt
        ydot[0] = y[1];  // dx/dt = v
        ydot[1] = mu_ * (1.0 - y[0] * y[0]) * y[1] - y[0];  // dv/dt = μ(1-x²)v - x
        return eNoError;
    }
    
    int getDimension() const override {
        return 2;
    }
};

/// @brief 性能测试结果结构
struct PerformanceResult {
    std::string integrator_name;
    double final_value;
    double error;
    long long time_ns;
    int steps;
    
    PerformanceResult(const std::string& name, double value, double err, long long time, int step_count)
        : integrator_name(name), final_value(value), error(err), time_ns(time), steps(step_count) {}
};

int main() {
    std::setlocale(LC_ALL, ".UTF-8");
    std::cout << "=== ODE积分器性能比较示例 ===" << std::endl;
    
    // 测试方程：范德波尔振荡器
    double mu = 5.0;  
    VanDerPolOscillator vdpODE(mu);
    
    // 积分参数
    double t0 = 0.0;
    double tf = 10.0;
    double y0[2] = {2.0, 0.0};  // x(0)=2, v(0)=0
    
    // 参考解（使用高精度积分器计算）
    RKF78 reference_integrator;
    reference_integrator.setMaxAbsErr(1e-12);
    reference_integrator.setMaxRelErr(1e-12);
    reference_integrator.setInitialStepSize(0.5);
    reference_integrator.initialize(vdpODE);
    
    double y_ref[2]{y0[0], y0[1]};
    double t = t0;
    errc_t err = reference_integrator.integrate(vdpODE, y_ref, t, tf);
    
    printf("err = %d\n", err);
    std::cout << "参考解: x(" << tf << ") = " << y_ref[0] << ", v(" << tf << ") = " << y_ref[1] << std::endl;
    std::cout << std::endl;
    
    // 测试的积分器列表
    std::vector<std::pair<std::string, ODEIntegrator*>> integrators;
    
    // 定步长积分器
    RK4* rk4 = new RK4();
    rk4->setStepSize(1e-2);
    integrators.emplace_back("RK4", rk4);
    RK8* rk8 = new RK8();
    rk8->setStepSize(1e-2);
    integrators.emplace_back("RK8", rk8);
    RKV8* rkv8 = new RKV8();
    rkv8->setStepSize(1e-2);
    integrators.emplace_back("RKV8", rkv8);
    
    // 变步长积分器
    RKF45* rkf45 = new RKF45();
    rkf45->setMaxAbsErr(1e-6);
    rkf45->setMaxRelErr(1e-6);
    integrators.emplace_back("RKF45", rkf45);
    
    RKF56* rkf56 = new RKF56();
    rkf56->setMaxAbsErr(1e-6);
    rkf56->setMaxRelErr(1e-6);
    integrators.emplace_back("RKF56", rkf56);
    
    RKF78* rkf78 = new RKF78();
    rkf78->setMaxAbsErr(1e-6);
    rkf78->setMaxRelErr(1e-6);
    rkf78->setInitialStepSize(0.5);
    integrators.emplace_back("RKF78", rkf78);
    
    RKCK* rkck = new RKCK();
    rkck->setMaxAbsErr(1e-6);
    rkck->setMaxRelErr(1e-6);
    integrators.emplace_back("RKCK", rkck);
    
    // 性能测试结果
    std::vector<PerformanceResult> results;
    
    for (auto& integrator_pair : integrators) {
        std::string name = integrator_pair.first;
        ODEIntegrator* integrator = integrator_pair.second;
        
        std::cout << "测试积分器: " << name << std::endl;
        
        // 初始化积分器
        integrator->initialize(vdpODE);
        
        // 测量执行时间
        auto start_time = std::chrono::high_resolution_clock::now();
        
        double yf[2]{y0[0], y0[1]};
        t = t0;
        errc_t result = integrator->integrate(vdpODE, yf, t, tf);
        
        auto end_time = std::chrono::high_resolution_clock::now();
        auto duration = std::chrono::duration_cast<std::chrono::nanoseconds>(end_time - start_time);
        
        if (result == eNoError) {
            // 计算误差
            double error_x = std::abs(yf[0] - y_ref[0]);
            double error_v = std::abs(yf[1] - y_ref[1]);
            double max_error = std::max(error_x, error_v);
            
            // 获取步数信息（对于变步长积分器）
            int step_count = 0;
            ODEVarStepIntegrator* var_integrator = dynamic_cast<ODEVarStepIntegrator*>(integrator);
            if (var_integrator != nullptr) {
                step_count = var_integrator->getNumSteps();
            } else {
                // 对于定步长积分器，估算步数
                step_count = static_cast<int>((tf - t0) / 0.01);  // 假设步长为0.01
            }
            
            results.emplace_back(name, yf[0], max_error, duration.count(), step_count);
            
            std::cout << "  结果: x=" << yf[0] << ", v=" << yf[1] << std::endl;
            std::cout << "  误差: " << max_error << std::endl;
            std::cout << "  时间: " << duration.count() / 1e6 << " ms" << std::endl;
            std::cout << "  步数: " << step_count << std::endl;
        } else {
            std::cout << "  积分失败" << std::endl;
        }
        
        std::cout << std::endl;
    }
    
    // 输出性能比较表格
    std::cout << "=== 性能比较结果 ===" << std::endl;
    std::cout << std::setw(10) << "积分器" 
              << std::setw(15) << "位置误差" 
              << std::setw(12) << "时间(ms)" 
              << std::setw(10) << "步数" 
              << std::setw(15) << "效率指标" << std::endl;
    std::cout << std::string(70, '-') << std::endl;
    
    for (const auto& result : results) {
        // 效率指标：误差越小、时间越短越好
        double efficiency = result.error * result.time_ns / 1e6;  // 误差×时间
        
        std::cout << std::setw(10) << result.integrator_name
                  << std::setw(15) << std::scientific << std::setprecision(2) << result.error
                  << std::setw(12) << std::fixed << std::setprecision(3) << result.time_ns / 1e6
                  << std::setw(10) << result.steps
                  << std::setw(15) << std::scientific << std::setprecision(2) << efficiency
                  << std::endl;
    }
    
    // 清理内存
    for (auto& integrator_pair : integrators) {
        delete integrator_pair.second;
    }
    
    // 不同步长对定步长积分器的影响
    std::cout << "\n=== 定步长积分器步长影响 ===" << std::endl;
    
    std::vector<double> step_sizes = {0.1, 0.05, 0.01, 0.005, 0.001};
    
    for (double step : step_sizes) {
        RK4 rk4;
        rk4.initialize(vdpODE);
        
        auto start_time = std::chrono::high_resolution_clock::now();
        
        // 手动执行定步长积分
        double t = t0;
        double y_current[2] = {y0[0], y0[1]};
        int step_count = 0;
        
        while (t < tf) {
            double actual_step = step;
            if (t + step > tf) {
                actual_step = tf - t;
            }
            
            rk4.singleStep(vdpODE, y_current, t, actual_step);
            
            t += actual_step;
            step_count++;
        }
        
        auto end_time = std::chrono::high_resolution_clock::now();
        auto duration = std::chrono::duration_cast<std::chrono::nanoseconds>(end_time - start_time);
        
        double error_x = std::abs(y_current[0] - y_ref[0]);
        double error_v = std::abs(y_current[1] - y_ref[1]);
        double max_error = std::max(error_x, error_v);
        
        std::cout << "步长=" << step << ": 误差=" << max_error 
                  << ", 时间=" << duration.count() / 1e6 << " ms"
                  << ", 步数=" << step_count << std::endl;
    }
    
    std::cout << "\n=== 示例完成 ===" << std::endl;
    
    return 0;
}

性能比较

积分器	阶数	类型	精度	计算量	适用场景
RK4	4	定步长	中等	中等	一般精度要求，计算量适中
RK8	8	定步长	高	大	高精度要求，可接受较大计算量
RKV8	8	定步长	高	大	高精度要求，数值稳定性好
RKF45	4(5)	变步长	中等	中等	一般精度，带误差控制
RKF56	5(6)	变步长	较高	较大	较高精度，带误差控制
RKF78	7(8)	变步长	高	大	高精度要求，带误差控制
RKCK	4(5)	变步长	中等	中等	一般精度，数值特性好

注意事项

步长选择: 定步长积分器需要合理选择步长，步长过大会导致精度不足，步长过小会增加计算量
容差设置: 变步长积分器的容差设置会影响计算精度和效率
内存使用: 高阶方法需要更多的中间变量存储，内存使用量较大

参考文献

Fehlberg, E. (1969). Low-order classical Runge-Kutta formulas with stepsize control and their application to some heat transfer problems.
Classical Fifth-, Sixth-, Seventh-, and Eighth-Order Runge-Kutta Formulas with Stepsize Control
Cash, J. R., & Karp, A. H. (1990). A variable order Runge-Kutta method for initial value problems with rapidly varying right-hand sides.
Verner, J. H. (1978). Explicit Runge-Kutta methods with estimates of the local truncation error.
Some Explicit Runge-Kutta Methods of High Order
Shanks, E. Baylis "Solutions of Differential Equations by Evaluations of Functions"

积分器

目录