headmyshoulder / odeint-v2

odeint - solving ordinary differential equations in c++ v2
http://headmyshoulder.github.com/odeint-v2/
Other
337 stars 102 forks source link

add support for complex<long double> #243

Open potatomama opened 4 years ago

potatomama commented 4 years ago

Hi. I have a task where my state vector is

typedef boost::array< std::complex<long double> , 4 > state_type;

This task does not compile. It looks like odeint does not work with "complex long double" at all.

To illustrate the problem, I modified the "stuart_landau.cpp" example, where I replaced "double" -> "long double". It would not compile. But if you change line 21 (see attached file) from

#define my_float long double

to

 #define my_float double

, it would compile just fine.

Is it a bug or I am missing something? Thanks!

/*
 * stuart_landau.cpp
 *
 * This example demonstrates how one can use odeint can be used with state types consisting of complex variables.
 *
 * Copyright 2011-2012 Karsten Ahnert
 * Copyright 2011-2013 Mario Mulansky
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 */

#include <iostream>
#include <complex>
#include <boost/array.hpp>

#include <boost/numeric/odeint.hpp>

#define my_float long double

using namespace std;
using namespace boost::numeric::odeint;
//[ stuart_landau_system_function
typedef complex< my_float > state_type;

struct stuart_landau
{
    my_float m_eta;
    my_float m_alpha;

    stuart_landau( my_float eta = 1.0 , my_float alpha = 1.0 )
    : m_eta( eta ) , m_alpha( alpha ) { }

    void operator()( const state_type &x , state_type &dxdt , my_float t ) const
    {
        my_float one=1.0;
        const complex< my_float > I( 0.0 , 1.0 );
        dxdt = ( one + m_eta * I ) * x - ( one + m_alpha * I ) * norm( x ) * x;
    }
};
//]

/*
//[ stuart_landau_system_function_alternative
my_float eta = 1.0;
my_float alpha = 1.0;
void stuart_landau( const state_type &x , state_type &dxdt , my_float t )
{
    const complex< my_float > I( 0.0 , 1.0 );
    dxdt = ( 1.0 + m_eta * I ) * x - ( 1.0 + m_alpha * I ) * norm( x ) * x;
}
//]
*/

struct streaming_observer
{
    std::ostream& m_out;

    streaming_observer( std::ostream &out ) : m_out( out ) { }

    template< class State >
    void operator()( const State &x , my_float t ) const
    {
        m_out << t;
        m_out << "\t" << x.real() << "\t" << x.imag() ;
        m_out << "\n";
    }
};

int main( int argc , char **argv )
{
    //[ stuart_landau_integration
    state_type x = complex< my_float >( 1.0 , 0.0 );

    const my_float dt = 0.1;

    typedef runge_kutta4< state_type > stepper_type;

    my_float zero=0.0;
    my_float ten=10.0;
    integrate_const( stepper_type() , stuart_landau( 2.0 , 1.0 ) , x , zero , ten , dt , streaming_observer( cout ) );
    //]

    return 0;
}
ddemidov commented 4 years ago

Should probably be runge_kutta4< state_type, my_float >.