c++ - odeint：复杂且精度更高的类型

Question

任何人都可以为我提供一些帮助，帮助我创建一个使用复数向量的 c++ odeint 求解器的最小示例，如果可能的话，精度更高（boost.multiprecision 或 libquadmath __float128 ， __complex128 说）。

文档中有一个使用复杂标量的示例http://headmyshoulder.github.com/odeint-v2/doc/boost_numeric_odeint/tutorial/special_topics.html

它在那里提到：

The fact that we have to configure a different algebra is solely due to the fact 
that we use a non-vector state type and not to the usage of complex values.
So for, e.g. vector<  complex<double> >, this would not be required. 

我试图通过以下更改来修改它typedef vector<complex<double>> state_type：

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

#include <boost/numeric/odeint.hpp>

using namespace std;
using namespace boost::numeric::odeint;


typedef vector<complex< double >> state_type;

struct stuart_landau
{
double m_eta;
double m_alpha;

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

void operator()( const state_type &x , state_type &dxdt , double t ) const
{
    const complex<double> I(0.0,1.0);
    dxdt[0] = x[1];
    dxdt[1] = (1.0+m_eta*I)*x[0]-(1.0+m_alpha*I)*x[1];
}
};

 struct streaming_observer
{
std::ostream& m_out;

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

template< class State >
void operator()( const State &x , double t ) const
{
    m_out << t;
    m_out << "\t" << x[0].real() << "\t" << x[0].imag() ;
    m_out << "\n";*/
}
};

int main( int argc , char **argv )
{
//[ stuart_landau_integration
state_type x(2);
 x[0] = complex< double >( 1.0 , 0.0 );
 x[1] = complex< double >( 1.0 , 0.0 );

const double dt = 0.1;

typedef runge_kutta4< state_type , double , state_type , double ,
                      vector_space_algebra > stepper_type;

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

return 0;

}

但是，这会产生大量错误。

谁能给我一个复杂 ODE 的最小工作示例？甚至更好地实现高精度数据类型..看起来 odeint 可以支持相当任意的状态类型，但我在让它们工作时遇到了很多麻烦。

更新代码

#include <iostream>
#include <complex>
#include <boost/array.hpp>
#include <stdlib.h>
#include <math.h>

#include <boost/numeric/odeint.hpp>


extern "C" {
#include <quadmath.h>
}

using namespace std;
using namespace boost::numeric::odeint;

//[ stuart_landau_system_function
typedef std::vector<std::complex < __float128 > > state_type;

struct stuart_landau
{
__float128 m_eta;
__float128 m_alpha;

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

void operator()( const state_type &x , state_type &dxdt , double t ) const
{
    const complex< __float128 > I( 0.0 , 1.0 ); //define complex I     

    dxdt[0] = x[1];
    dxdt[1] = ( 1.0 + m_eta * I ) * x[0] - ( 1.0 + m_alpha * I )*x[1];
}
};
//]




struct streaming_observer
{
 std::ostream& m_out;

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

template < class State >
void operator()( const State &x , double t ) const
{
    m_out << t;
  /*  m_out << "\t" << x[0].real() << "\t" << x[0].imag() ;
    m_out << "\n";*/
 }
};




int main( int argc , char **argv )
{
//[ stuart_landau_integration
state_type x(2);
 x[0] = complex< __float128 >( 1.0L , 0.0L );
 x[1] = complex< __float128 >( 1.0L , 0.0L );

const double dt = 0.1;

typedef runge_kutta4< state_type, __float128 > stepper_type;

integrate_const( stepper_type() , stuart_landau( 2.0L , 1.0L ) , x , 0.0 , 10.0 , dt    
 /*, streaming_observer( cout ) */);
 //]

 return 0;
 }

错误登录编译：

vecPrec.cpp: In member function ‘void stuart_landau::operator()(const state_type&,    
state_type&, double) const’:
 vecPrec.cpp:33:35: error: no match for ‘operator+’ in ‘1.0e+0 + std::operator* [with 
 _Tp = __float128]((* &((const stuart_landau*)this)->stuart_landau::m_eta), (* & I))’
 vecPrec.cpp:33:35: note: candidates are:
 /usr/include/c++/4.6/bits/stl_iterator.h:327:5: note: template<class _Iterator>    std::reverse_iterator<_Iterator> std::operator+(typename 

 std::reverse_iterator<_Iterator>::difference_type, const   std::reverse_iterator<_Iterator>&)
 /usr/include/c++/4.6/bits/basic_string.h:2306:5: note: template<class _CharT, class  _Traits, class _Alloc> std::basic_string<_CharT, _Traits, _Alloc> std::operator+(const    std::basic_string<_CharT, _Traits, _Alloc>&, const std::basic_string<_CharT, _Traits,   _Alloc>&)
 /usr/include/c++/4.6/bits/basic_string.tcc:694:5: note: template<class _CharT, class  _Traits, class _Alloc> std::basic_string<_CharT, _Traits, _Alloc> std::operator+(const   _CharT*, const std::basic_string<_CharT, _Traits, _Alloc>&)
  /usr/include/c++/4.6/bits/basic_string.tcc:710:5: note: template<class _CharT, class _Traits, class _Alloc> std::basic_string<_CharT, _Traits, _Alloc> std::operator+(_CharT,     const std::basic_string<_CharT, _Traits, _Alloc>&)
  /usr/include/c++/4.6/bits/basic_string.h:2343:5: note: template<class _CharT, class  _Traits, class _Alloc> std::basic_string<_CharT, _Traits, _Alloc> std::operator+(const   std::basic_string<_CharT, _Traits, _Alloc>&, const _CharT*)
 /usr/include/c++/4.6/bits/basic_string.h:2359:5: note: template<class _CharT, class   _Traits, class _Alloc> std::basic_string<_CharT, _Traits, _Alloc> std::operator+(const   std::basic_string<_CharT, _Traits, _Alloc>&, _CharT)
  /usr/include/c++/4.6/complex:321:5: note: template<class _Tp> std::complex<_Tp>  std::operator+(const std::complex<_Tp>&, const std::complex<_Tp>&)
   /usr/include/c++/4.6/complex:330:5: note: template<class _Tp> std::complex<_Tp>  std::operator+(const std::complex<_Tp>&, const _Tp&)
    /usr/include/c++/4.6/complex:339:5: note: template<class _Tp> std::complex<_Tp>  std::operator+(const _Tp&, const std::complex<_Tp>&)
   /usr/include/c++/4.6/complex:440:5: note: template<class _Tp> std::complex<_Tp>  std::operator+(const std::complex<_Tp>&)
   /usr/include/c++/4.6/bits/stl_bvector.h:266:3: note: std::_Bit_iterator   std::operator+(std::ptrdiff_t, const std::_Bit_iterator&)
   /usr/include/c++/4.6/bits/stl_bvector.h:266:3: note:   no known conversion for argument 2 from ‘std::complex<__float128>’ to ‘const std::_Bit_iterator&’
    /usr/include/c++/4.6/bits/stl_bvector.h:352:3: note: std::_Bit_const_iterator std::operator+(std::ptrdiff_t, const std::_Bit_const_iterator&)
   /usr/include/c++/4.6/bits/stl_bvector.h:352:3: note:   no known conversion for   argument 2 from ‘std::complex<__float128>’ to ‘const std::_Bit_const_iterator&’
   vecPrec.cpp:33:66: error: no match for ‘operator+’ in ‘1.0e+0 + std::operator*  [with _Tp = __float128]((* &((const stuart_landau*)this)->stuart_landau::m_alpha), (* &   I))’
  vecPrec.cpp:33:66: note: candidates are:
   /usr/include/c++/4.6/bits/stl_iterator.h:327:5: note: template<class _Iterator>   std::reverse_iterator<_Iterator> std::operator+(typename    std::reverse_iterator<_Iterator>::difference_type, const std::reverse_iterator<_Iterator>&)
   /usr/include/c++/4.6/bits/basic_string.h:2306:5: note: template<class _CharT, class     _Traits, class _Alloc> std::basic_string<_CharT, _Traits, _Alloc> std::operator+   (const std::basic_string<_CharT, _Traits, _Alloc>&, const   std::basic_string<_CharT, _Traits, _Alloc>&)
  /usr/include/c++/4.6/bits/basic_string.tcc:694:5: note: template<class _CharT, class   _Traits, class _Alloc> std::basic_string<_CharT, _Traits, _Alloc> std::operator+  (const _CharT*, const std::basic_string<_CharT, _Traits, _Alloc>&)
   /usr/include/c++/4.6/bits/basic_string.tcc:710:5: note: template<class _CharT, class _Traits, class _Alloc> std::basic_string<_CharT, _Traits, _Alloc> std::operator+   (_CharT, const std::basic_string<_CharT, _Traits, _Alloc>&)
    /usr/include/c++/4.6/bits/basic_string.h:2343:5: note: template<class _CharT, class _Traits, class _Alloc> std::basic_string<_CharT, _Traits, _Alloc> std::operator+(const std::basic_string<_CharT, _Traits, _Alloc>&, const _CharT*)
   /usr/include/c++/4.6/bits/basic_string.h:2359:5: note: template<class _CharT, class _Traits, class _Alloc> std::basic_string<_CharT, _Traits, _Alloc> std::operator+(const std::basic_string<_CharT, _Traits, _Alloc>&, _CharT)
   /usr/include/c++/4.6/complex:321:5: note: template<class _Tp> std::complex<_Tp> std::operator+(const std::complex<_Tp>&, const std::complex<_Tp>&)
    /usr/include/c++/4.6/complex:330:5: note: template<class _Tp> std::complex<_Tp> std::operator+(const std::complex<_Tp>&, const _Tp&)
   /usr/include/c++/4.6/complex:339:5: note: template<class _Tp> std::complex<_Tp> std::operator+(const _Tp&, const std::complex<_Tp>&)
   /usr/include/c++/4.6/complex:440:5: note: template<class _Tp> std::complex<_Tp> std::operator+(const std::complex<_Tp>&)
   /usr/include/c++/4.6/bits/stl_bvector.h:266:3: note: std::_Bit_iterator std::operator+(std::ptrdiff_t, const std::_Bit_iterator&)
   /usr/include/c++/4.6/bits/stl_bvector.h:266:3: note:   no known conversion for argument 2 from ‘std::complex<__float128>’ to ‘const std::_Bit_iterator&’
   /usr/include/c++/4.6/bits/stl_bvector.h:352:3: note: std::_Bit_const_iterator std::operator+(std::ptrdiff_t, const std::_Bit_const_iterator&)
   /usr/include/c++/4.6/bits/stl_bvector.h:352:3: note:   no known conversion for argument 2 from ‘std::complex<__float128>’ to ‘const std::_Bit_const_iterator&’

Answer 1

为了使上面的示例正常工作，您需要使用range_algebra. 您正在使用vector_space_algebrawhich 期望或多或少地为state_type所有运算符 +*-/ 定义。情况并非如此std::vector<>。只需使用

typedef runge_kutta4< state_type > stepper_type;

此 typedef 静默使用range_algebra.

对于另一种情况，请尝试

typedef vector< complex< __float128 > > state_type;
typedef runge_kutta4< state_type , __float128 > stepper_type;

上述两行工作的前提是 __float128 的操作符-*/+ 已经定义。