/* * two_dimensional_phase_lattice.cpp * * This example show how one can use matrices as state types in odeint. * * Copyright 2011-2012 Karsten Ahnert * Copyright 2011-2013 Mario Mulansky * 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 <map> #include <string> #include <fstream> #ifndef M_PI //not there on windows #define M_PI 3.1415927 //... #endif #include <boost/numeric/odeint.hpp> using namespace std; using namespace boost::numeric::odeint; //[ two_dimensional_phase_lattice_definition typedef boost::numeric::ublas::matrix< double > state_type; struct two_dimensional_phase_lattice { two_dimensional_phase_lattice( double gamma = 0.5 ) : m_gamma( gamma ) { } void operator()( const state_type &x , state_type &dxdt , double /* t */ ) const { size_t size1 = x.size1() , size2 = x.size2(); for( size_t i=1 ; i<size1-1 ; ++i ) { for( size_t j=1 ; j<size2-1 ; ++j ) { dxdt( i , j ) = coupling_func( x( i + 1 , j ) - x( i , j ) ) + coupling_func( x( i - 1 , j ) - x( i , j ) ) + coupling_func( x( i , j + 1 ) - x( i , j ) ) + coupling_func( x( i , j - 1 ) - x( i , j ) ); } } for( size_t i=0 ; i<x.size1() ; ++i ) dxdt( i , 0 ) = dxdt( i , x.size2() -1 ) = 0.0; for( size_t j=0 ; j<x.size2() ; ++j ) dxdt( 0 , j ) = dxdt( x.size1() -1 , j ) = 0.0; } double coupling_func( double x ) const { return sin( x ) - m_gamma * ( 1.0 - cos( x ) ); } double m_gamma; }; //] struct write_for_gnuplot { size_t m_every , m_count; write_for_gnuplot( size_t every = 10 ) : m_every( every ) , m_count( 0 ) { } void operator()( const state_type &x , double t ) { if( ( m_count % m_every ) == 0 ) { clog << t << endl; cout << "sp '-'" << endl; for( size_t i=0 ; i<x.size1() ; ++i ) { for( size_t j=0 ; j<x.size2() ; ++j ) { cout << i << "\t" << j << "\t" << sin( x( i , j ) ) << "\n"; } cout << "\n"; } cout << "e" << endl; } ++m_count; } }; class write_snapshots { public: typedef std::map< size_t , std::string > map_type; write_snapshots( void ) : m_count( 0 ) { } void operator()( const state_type &x , double t ) { map< size_t , string >::const_iterator it = m_snapshots.find( m_count ); if( it != m_snapshots.end() ) { ofstream fout( it->second.c_str() ); for( size_t i=0 ; i<x.size1() ; ++i ) { for( size_t j=0 ; j<x.size2() ; ++j ) { fout << i << "\t" << j << "\t" << x( i , j ) << "\t" << sin( x( i , j ) ) << "\n"; } fout << "\n"; } } ++m_count; } map_type& snapshots( void ) { return m_snapshots; } const map_type& snapshots( void ) const { return m_snapshots; } private: size_t m_count; map_type m_snapshots; }; int main( int argc , char **argv ) { size_t size1 = 128 , size2 = 128; state_type x( size1 , size2 , 0.0 ); for( size_t i=(size1/2-10) ; i<(size1/2+10) ; ++i ) for( size_t j=(size2/2-10) ; j<(size2/2+10) ; ++j ) x( i , j ) = static_cast<double>( rand() ) / RAND_MAX * 2.0 * M_PI; write_snapshots snapshots; snapshots.snapshots().insert( make_pair( size_t( 0 ) , string( "lat_0000.dat" ) ) ); snapshots.snapshots().insert( make_pair( size_t( 100 ) , string( "lat_0100.dat" ) ) ); snapshots.snapshots().insert( make_pair( size_t( 1000 ) , string( "lat_1000.dat" ) ) ); observer_collection< state_type , double > obs; obs.observers().push_back( write_for_gnuplot( 10 ) ); obs.observers().push_back( snapshots ); cout << "set term x11" << endl; cout << "set pm3d map" << endl; integrate_const( runge_kutta4<state_type>() , two_dimensional_phase_lattice( 1.2 ) , x , 0.0 , 1001.0 , 0.1 , boost::ref( obs ) ); // controlled steppers work only after ublas bugfix //integrate_const( make_dense_output< runge_kutta_dopri5< state_type > >( 1E-6 , 1E-6 ) , two_dimensional_phase_lattice( 1.2 ) , // x , 0.0 , 1001.0 , 0.1 , boost::ref( obs ) ); return 0; }