alesolano
commented
7 years ago Fixed. I don't know exactly why. I'm just learning how to use this tool, but I simplified the example code to only use the basic structure: 'rhs' as a void function instead of a class, and only integrate with the 'integrate observer'.
/*
Copyright 2010-2012 Karsten Ahnert
Copyright 2011-2013 Mario Mulansky
Copyright 2013 Pascal Germroth
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 <vector>
#include <boost/numeric/odeint.hpp>
//[ rhs_function
/* The type of container used to hold the state vector */
typedef std::vector< double > state_type;
const double gam = 0.15;
/* The rhs of x' = f(x) */
void harmonic_oscillator( const state_type &x , state_type &dxdt , const double /* t */ )
{
dxdt[0] = x[1];
dxdt[1] = -x[0] - gam*x[1];
}
//]
//[ integrate_observer
struct push_back_state_and_time
{
std::vector< state_type >& m_states;
std::vector< double >& m_times;
push_back_state_and_time( std::vector< state_type > &states , std::vector< double > × )
: m_states( states ) , m_times( times ) { }
void operator()( const state_type &x , double t )
{
m_states.push_back( x );
m_times.push_back( t );
}
};
//]
int main(int /* argc */ , char** /* argv */ )
{
using namespace std;
using namespace boost::numeric::odeint;
//[ state_initialization
state_type x(2);
x[0] = 1.0; // start at x=1.0, p=0.0
x[1] = 0.0;
//]
//[ integrate_observ
vector<state_type> x_vec;
vector<double> times;
size_t steps = integrate( harmonic_oscillator ,
x , 0.0 , 10.0 , 0.1 ,
push_back_state_and_time( x_vec , times ) );
/* output */
for( size_t i=0; i<=steps; i++ )
{
cout << times[i] << '\t' << x_vec[i][0] << '\t' << x_vec[i][1] << '\n';
}
//]
return 0;
}The new output is:
0 1 0
0.1 0.995029 -0.0990884
0.342911 0.942763 -0.32773
0.59639 0.832373 -0.537274
0.865439 0.662854 -0.714058
1.15445 0.436595 -0.839871
1.44346 0.184494 -0.892171
1.70128 -0.0446646 -0.875731
1.9591 -0.262234 -0.80313
2.21692 -0.454545 -0.681207
2.48815 -0.616993 -0.510476
2.77613 -0.733844 -0.296963
3.06411 -0.7866 -0.0685387
3.35209 -0.773719 0.155753
3.64008 -0.699015 0.358007
3.92806 -0.571123 0.522852
4.21604 -0.402601 0.638582
4.50402 -0.20875 0.697943
4.792 -0.0062679 0.698535
5.07998 0.188159 0.642795
5.36797 0.359199 0.537596
5.65595 0.494053 0.39351
5.94393 0.583381 0.223796
6.23191 0.621907 0.0432166
6.51989 0.608673 -0.133212
6.80787 0.546935 -0.291443
7.09586 0.44372 -0.419492
7.39491 0.303457 -0.510857
7.69397 0.143005 -0.553869
7.99302 -0.0228226 -0.546913
8.29208 -0.179394 -0.492786
8.59113 -0.313525 -0.398255
8.89019 -0.414555 -0.273299
9.18925 -0.475166 -0.130101
9.4883 -0.49187 0.0181101
9.78736 -0.465143 0.158247
10 -0.421907 0.246407
This is my output for examples/harmonic_oscillator.cpp using CodeBlocks IDE in Windows with MinGW compiler set to C++14. I didn't change anything of the code.
The initial condition is [1.0, 0.0], but my output shows that it starts at [0.117289, -0.198813]. So I did something wrong...
I have tested it in MATLAB and the output is:
(I'm in the middle of switching from MATLAB ode solvers to something open-source)
Thank you!