1 | /* Boost numeric test of the rosenbrock4 stepper test file |
2 | |
3 | Copyright 2012 Karsten Ahnert |
4 | Copyright 2012 Mario Mulansky |
5 | |
6 | Distributed under the Boost Software License, Version 1.0. |
7 | (See accompanying file LICENSE_1_0.txt or |
8 | copy at http://www.boost.org/LICENSE_1_0.txt) |
9 | */ |
10 | |
11 | // disable checked iterator warning for msvc |
12 | #include <boost/config.hpp> |
13 | #ifdef BOOST_MSVC |
14 | #pragma warning(disable:4996) |
15 | #endif |
16 | |
17 | #define BOOST_TEST_MODULE numeric_rosenbrock |
18 | |
19 | #include <iostream> |
20 | #include <cmath> |
21 | |
22 | #include <array> |
23 | |
24 | #include <boost/test/unit_test.hpp> |
25 | |
26 | #include <boost/numeric/odeint/stepper/rosenbrock4.hpp> |
27 | #include <boost/numeric/ublas/vector.hpp> |
28 | #include <boost/numeric/ublas/matrix.hpp> |
29 | |
30 | using namespace boost::unit_test; |
31 | using namespace boost::numeric::odeint; |
32 | |
33 | typedef double value_type; |
34 | typedef boost::numeric::ublas::vector< value_type > state_type; |
35 | typedef boost::numeric::ublas::matrix< value_type > matrix_type; |
36 | |
37 | // harmonic oscillator, analytic solution x[0] = sin( t ) |
38 | struct sys |
39 | { |
40 | void operator()( const state_type &x , state_type &dxdt , const value_type &t ) const |
41 | { |
42 | dxdt( 0 ) = x( 1 ); |
43 | dxdt( 1 ) = -x( 0 ); |
44 | } |
45 | }; |
46 | |
47 | struct jacobi |
48 | { |
49 | void operator()( const state_type &x , matrix_type &jacobi , const value_type &t , state_type &dfdt ) const |
50 | { |
51 | jacobi( 0 , 0 ) = 0; |
52 | jacobi( 0 , 1 ) = 1; |
53 | jacobi( 1 , 0 ) = -1; |
54 | jacobi( 1 , 1 ) = 0; |
55 | dfdt( 0 ) = 0.0; |
56 | dfdt( 1 ) = 0.0; |
57 | } |
58 | }; |
59 | |
60 | |
61 | BOOST_AUTO_TEST_SUITE( numeric_rosenbrock4 ) |
62 | |
63 | BOOST_AUTO_TEST_CASE( rosenbrock4_numeric_test ) |
64 | { |
65 | typedef rosenbrock4< value_type > stepper_type; |
66 | stepper_type stepper; |
67 | |
68 | const int o = stepper.order()+1; |
69 | |
70 | state_type x0( 2 ) , x1( 2 ); |
71 | x0(0) = 0.0; x0(1) = 1.0; |
72 | |
73 | double dt = 0.5; |
74 | |
75 | stepper.do_step( system: std::make_pair( x: sys() , y: jacobi() ) , x: x0 , t: 0.0 , xout&: x1 , dt ); |
76 | const double f = 2.0 * std::abs( x: std::sin(x: dt) - x1(0) ) / std::pow( x: dt , y: o ); |
77 | |
78 | std::cout << o << " , " << f << std::endl; |
79 | |
80 | while( f*std::pow( x: dt , y: o ) > 1E-16 ) |
81 | { |
82 | stepper.do_step( system: std::make_pair( x: sys() , y: jacobi() ) , x: x0 , t: 0.0 , xout&: x1 , dt ); |
83 | std::cout << "Testing dt=" << dt << std::endl; |
84 | BOOST_CHECK_SMALL( std::abs( std::sin(dt) - x1(0) ) , f*std::pow( dt , o ) ); |
85 | dt *= 0.5; |
86 | } |
87 | } |
88 | |
89 | BOOST_AUTO_TEST_SUITE_END() |
90 | |