Brian Silverman | 7c33ab2 | 2018-08-04 17:14:51 -0700 | [diff] [blame^] | 1 | /* Boost check_gmp.cpp test file |
| 2 | |
| 3 | Copyright 2010-2011 Karsten Ahnert |
| 4 | Copyright 2011 Mario Mulansky |
| 5 | |
| 6 | This file tests the odeint library with the gmp arbitrary precision types |
| 7 | |
| 8 | Distributed under the Boost Software License, Version 1.0. |
| 9 | (See accompanying file LICENSE_1_0.txt or |
| 10 | copy at http://www.boost.org/LICENSE_1_0.txt) |
| 11 | */ |
| 12 | |
| 13 | #define BOOST_TEST_MODULE odeint_gsl |
| 14 | |
| 15 | #include <gsl/gsl_vector.h> |
| 16 | #include <boost/numeric/odeint/stepper/euler.hpp> |
| 17 | #include <boost/numeric/odeint/external/gsl/gsl_wrapper.hpp> |
| 18 | |
| 19 | #include <boost/test/unit_test.hpp> |
| 20 | |
| 21 | using namespace boost::unit_test; |
| 22 | using namespace boost::numeric::odeint; |
| 23 | |
| 24 | typedef gsl_vector *state_type; |
| 25 | |
| 26 | const double sigma = 10.0; |
| 27 | const double R = 28.0; |
| 28 | const double b = 8.0 / 3.0; |
| 29 | |
| 30 | void lorenz( const state_type x , state_type dxdt , double t ) |
| 31 | { |
| 32 | gsl_vector_set( dxdt , 0 , sigma * ( gsl_vector_get(x , 1 ) - gsl_vector_get( x , 0 ) ) ); |
| 33 | gsl_vector_set( dxdt , 1 , R * gsl_vector_get( x , 0 ) - gsl_vector_get( x , 1 ) - gsl_vector_get( x , 0 ) * gsl_vector_get( x , 2) ); |
| 34 | gsl_vector_set( dxdt , 2 , gsl_vector_get( x , 0 ) * gsl_vector_get( x , 1 ) - b * gsl_vector_get( x , 2) ); |
| 35 | } |
| 36 | |
| 37 | BOOST_AUTO_TEST_CASE( gsl ) |
| 38 | { |
| 39 | euler< state_type > euler; |
| 40 | |
| 41 | state_type x = gsl_vector_alloc( 3 ); |
| 42 | |
| 43 | // check resizing |
| 44 | state_type y = 0; |
| 45 | boost::numeric::odeint::resize( y , x ); |
| 46 | BOOST_CHECK( 0 != y ); |
| 47 | |
| 48 | gsl_vector_set( x , 0 , 1.0); |
| 49 | gsl_vector_set( x , 1 , 1.0); |
| 50 | gsl_vector_set( x , 2 , 2.0); |
| 51 | |
| 52 | euler.do_step( lorenz , x , 0.0 , 0.1 ); |
| 53 | |
| 54 | //cout << gsl_vector_get( x , 0 ) << " " << gsl_vector_get( x , 1 ) << " " << gsl_vector_get( x , 2 ) << endl; |
| 55 | |
| 56 | gsl_vector_free( x ); |
| 57 | |
| 58 | } |