Output of Harmonic Oscillator example not reasonable
alesolano opened this issue · 1 comments
alesolano commented
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.
0 0.117289 -0.198813
0.1 0.0970056 -0.206492
0.427867 0.0270349 -0.216419
0.755734 -0.0423883 -0.203372
1.0836 -0.104073 -0.169808
1.41147 -0.151978 -0.120201
1.7666 -0.183363 -0.0552639
2.12174 -0.190834 0.0130209
2.47687 -0.174725 0.0761808
2.832 -0.138217 0.126806
3.18714 -0.0868043 0.159393
3.54227 -0.0275093 0.170918
3.8974 0.032053 0.161073
4.25254 0.0846297 0.132163
4.60767 0.124194 0.0886899
4.97682 0.147111 0.0345601
5.34597 0.149511 -0.0211835
5.71512 0.132187 -0.0711478
6.08427 0.0984694 -0.109095
6.45342 0.0536622 -0.1307
6.82257 0.00423529 -0.134009
7.19172 -0.0430842 -0.119576
7.56087 -0.0822076 -0.0902481
7.93002 -0.108449 -0.0506709
8.29917 -0.119062 -0.00658431
8.66832 -0.113498 0.0359904
9.03747 -0.0933576 0.0715544
9.40662 -0.062066 0.095824
9.80128 -0.0216132 0.106406
10 -0.000495785 0.105432
0 0.000354255 0.000398615
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:
0 1.0000 0
0.0001 1.0000 -0.0001
0.0005 1.0000 -0.0005
0.0025 1.0000 -0.0025
0.0125 0.9999 -0.0125
0.0625 0.9981 -0.0621
0.1793 0.9841 -0.1760
0.3471 0.9414 -0.3314
0.5609 0.8509 -0.5101
0.8230 0.6925 -0.6895
1.1192 0.4657 -0.8278
1.3839 0.2374 -0.8868
1.6144 0.0321 -0.8869
1.8012 -0.1309 -0.8531
1.9583 -0.2613 -0.8026
2.1625 -0.4161 -0.7098
2.4138 -0.5763 -0.5602
2.6989 -0.7074 -0.3560
2.9540 -0.7728 -0.1563
3.1238 -0.7879 -0.0212
3.2937 -0.7801 0.1112
3.4483 -0.7539 0.2260
3.6499 -0.6942 0.3635
3.8984 -0.5853 0.5068
4.2034 -0.4098 0.6330
4.4758 -0.2280 0.6928
4.7147 -0.0605 0.7021
4.8564 0.0382 0.6889
4.9981 0.1341 0.6622
5.1714 0.2448 0.6126
5.3907 0.3701 0.5260
5.6592 0.4936 0.3906
5.9345 0.5792 0.2290
6.1788 0.6165 0.0766
6.3250 0.6210 -0.0147
6.4712 0.6123 -0.1037
6.6356 0.5873 -0.1988
6.8463 0.5335 -0.3092
7.1049 0.4383 -0.4211
7.4033 0.2981 -0.5105
7.6700 0.1558 -0.5498
7.9027 0.0270 -0.5517
8.0927 -0.0762 -0.5314
8.2463 -0.1556 -0.5016
8.4471 -0.2510 -0.4463
8.6946 -0.3508 -0.3563
8.9819 -0.4353 -0.2297
9.2392 -0.4784 -0.1050
9.4133 -0.4892 -0.0189
9.5874 -0.4851 0.0655
9.7390 -0.4698 0.1357
9.9379 -0.4341 0.2205
10.0000 -0.4197 0.2448
(I'm in the middle of switching from MATLAB ode solvers to something open-source)
Thank you!
alesolano commented
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