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

Question

I would like to integrate a function that maps a 3D point (parametrized by t) to a 2D point (complex plane) using an adaptative step scheme. There is no closed-form for the derivative of my function and it is non-linear.

I've tried the following to see if the code would work. It compiles but the result is wrong. The test function being integrated (t from 0 to 1; i is the complex number) is

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

The expected result is

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

The code outputs

5051 + 0 i

I suspect the problem is in my definition of x, the state vector. I don't know what it should be.

Answer 1

I suspect your problem is because you are modifying point_ everytime you call Integrand::operator().

Instead of:

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

You probably meant:

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

Your Integrand::operator() should be marked as a const function when you don't member variables to change, that would help catch these errors in the future.

After looking at the docs for odeint, integrate_adaptive returns the number of steps performed. The input parameter x actually holds the final result so you want to do:

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

Running this prints (0.782859,-0.279002), which is still not the answer you're looking for. The answer you're looking for comes as a result of starting x at 0 instead of 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;