使用 boost::numeric::odeint 对非线性函数求积分 f'(x, y, z) = a + b*I

Question

我想集成一个函数，该函数使用自适应步骤方案将 3D 点（由 t 参数化）映射到 2D 点（复平面）。我的函数的导数没有封闭形式，它是非线性的。

我已经尝试了以下方法以查看代码是否有效。它编译但结果是错误的。被积分的测试函数（t从0到1；i是复数）是

Exp[ -Norm[ {1.1, 2.4, 3.6}*t ] * i ]

预期结果是

-0.217141 - 0.279002 i

#include <iostream>
#include <complex>
#include <boost/numeric/ublas/vector.hpp>
namespace ublas = boost::numeric::ublas;
#include <boost/numeric/odeint.hpp>
namespace odeint = boost::numeric::odeint;

typedef std::complex<double> state_type;

class Integrand {
    ublas::vector<double> point_;
public:
    Integrand(ublas::vector<double> point){
        point_ = point;
    }
    void operator () (const state_type &x, state_type &dxdt, const double t){
        point_ *= t;
        const std::complex<double> I(0.0, 1.0);
        dxdt = std::exp( -norm_2(point_)*I );
    }
};

std::complex<double> integral(ublas::vector<double> pt) {
    typedef odeint::runge_kutta_cash_karp54< state_type > error_stepper_type;
    double err_abs = 1.0e-10;
    double err_rel = 1.0e-6;
    state_type x = std::complex<double>(1.0, 0.0);
    return odeint::integrate_adaptive(
               odeint::make_controlled<error_stepper_type>(err_abs, err_rel),
               Integrand(pt), x, 0.0, 1.0, 0.001);
}

int main() {
    ublas::vector<double> pt(3);
    pt(0) = 1.1;
    pt(1) = 2.4;
    pt(2) = 3.6;
    std::cout << integral(pt) << std::endl;
    return 0;
}

代码输出

5051 + 0 i

我怀疑问题出在我对状态向量 x 的定义中。我不知道它应该是什么。

Answer 1

我怀疑你的问题是因为你每次调用 Integrand::operator().

时都在修改 point_

而不是：

point_ *= t;
dxdt = exp(-norm_2(point_)*I);

你的意思可能是：

dxdt = exp(-norm_2(point_ * t) * I);

当您不更改要更改的成员变量时，您的 Integrand::operator() 应标记为 const 函数，这将有助于在将来捕获这些错误。

查看 odeint 的文档后，integrate_adaptive returns 执行的步骤数。输入参数 x 实际上保存了最终结果所以你想做：

odeint::integrate_adaptive(
           odeint::make_controlled<error_stepper_type>(err_abs, err_rel),
           Integrand(pt), x, 0.0, 1.0, 0.001);
return x;

运行这会打印出 (0.782859,-0.279002)，这仍然不是您要查找的答案。您正在寻找的答案是从 x 开始于 0 而不是 1.

state_type x = std::complex<double>(0.0, 0.0);
odeint::integrate_adaptive(
           odeint::make_controlled<error_stepper_type>(err_abs, err_rel),
           Integrand(pt), x, 0.0, 1.0, 0.001);
return x;

使用 boost::numeric::odeint 对非线性函数求积分 f'(x, y, z) = a + b*I

Using boost::numeric::odeint to integrate a non-linear function f'(x, y, z) = a + b*I

c++

boost

odeint