1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
|
/*
* resizing_lattice.cpp
*
* Demonstrates the usage of resizing of the state type during integration.
* Examplary system is a strongly nonlinear, disordered Hamiltonian lattice
* where the spreading of energy is investigated
*
* Copyright 2011-2012 Mario Mulansky
* Copyright 2012-2013 Karsten Ahnert
* 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 <utility>
#include <boost/numeric/odeint.hpp>
#include <boost/ref.hpp>
#include <boost/random.hpp>
using namespace std;
using namespace boost::numeric::odeint;
//[ resizing_lattice_system_class
typedef vector< double > coord_type;
typedef pair< coord_type , coord_type > state_type;
struct compacton_lattice
{
const int m_max_N;
const double m_beta;
int m_pot_start_index;
vector< double > m_pot;
compacton_lattice( int max_N , double beta , int pot_start_index )
: m_max_N( max_N ) , m_beta( beta ) , m_pot_start_index( pot_start_index ) , m_pot( max_N )
{
srand( time( NULL ) );
// fill random potential with iid values from [0,1]
boost::mt19937 rng;
boost::uniform_real<> unif( 0.0 , 1.0 );
boost::variate_generator< boost::mt19937&, boost::uniform_real<> > gen( rng , unif );
generate( m_pot.begin() , m_pot.end() , gen );
}
void operator()( const coord_type &q , coord_type &dpdt )
{
// calculate dpdt = -dH/dq of this hamiltonian system
// dp_i/dt = - V_i * q_i^3 - beta*(q_i - q_{i-1})^3 + beta*(q_{i+1} - q_i)^3
const int N = q.size();
double diff = q[0] - q[N-1];
for( int i=0 ; i<N ; ++i )
{
dpdt[i] = - m_pot[m_pot_start_index+i] * q[i]*q[i]*q[i] -
m_beta * diff*diff*diff;
diff = q[(i+1) % N] - q[i];
dpdt[i] += m_beta * diff*diff*diff;
}
}
void energy_distribution( const coord_type &q , const coord_type &p , coord_type &energies )
{
// computes the energy per lattice site normalized by total energy
const size_t N = q.size();
double en = 0.0;
for( size_t i=0 ; i<N ; i++ )
{
const double diff = q[(i+1) % N] - q[i];
energies[i] = p[i]*p[i]/2.0
+ m_pot[m_pot_start_index+i]*q[i]*q[i]*q[i]*q[i]/4.0
+ m_beta/4.0 * diff*diff*diff*diff;
en += energies[i];
}
en = 1.0/en;
for( size_t i=0 ; i<N ; i++ )
{
energies[i] *= en;
}
}
double energy( const coord_type &q , const coord_type &p )
{
// calculates the total energy of the excitation
const size_t N = q.size();
double en = 0.0;
for( size_t i=0 ; i<N ; i++ )
{
const double diff = q[(i+1) % N] - q[i];
en += p[i]*p[i]/2.0
+ m_pot[m_pot_start_index+i]*q[i]*q[i]*q[i]*q[i] / 4.0
+ m_beta/4.0 * diff*diff*diff*diff;
}
return en;
}
void change_pot_start( const int delta )
{
m_pot_start_index += delta;
}
};
//]
//[ resizing_lattice_resize_function
void do_resize( coord_type &q , coord_type &p , coord_type &distr , const int N )
{
q.resize( N );
p.resize( N );
distr.resize( N );
}
//]
const int max_N = 1024;
const double beta = 1.0;
int main()
{
//[ resizing_lattice_initialize
//start with 60 sites
const int N_start = 60;
coord_type q( N_start , 0.0 );
q.reserve( max_N );
coord_type p( N_start , 0.0 );
p.reserve( max_N );
// start with uniform momentum distribution over 20 sites
fill( p.begin()+20 , p.end()-20 , 1.0/sqrt(20.0) );
coord_type distr( N_start , 0.0 );
distr.reserve( max_N );
// create the system
compacton_lattice lattice( max_N , beta , (max_N-N_start)/2 );
//create the stepper, note that we use an always_resizer because state size might change during steps
typedef symplectic_rkn_sb3a_mclachlan< coord_type , coord_type , double , coord_type , coord_type , double ,
range_algebra , default_operations , always_resizer > hamiltonian_stepper;
hamiltonian_stepper stepper;
hamiltonian_stepper::state_type state = make_pair( q , p );
//]
//[ resizing_lattice_steps_loop
double t = 0.0;
const double dt = 0.1;
const int steps = 10000;
for( int step = 0 ; step < steps ; ++step )
{
stepper.do_step( boost::ref(lattice) , state , t , dt );
lattice.energy_distribution( state.first , state.second , distr );
if( distr[10] > 1E-150 )
{
do_resize( state.first , state.second , distr , state.first.size()+20 );
rotate( state.first.begin() , state.first.end()-20 , state.first.end() );
rotate( state.second.begin() , state.second.end()-20 , state.second.end() );
lattice.change_pot_start( -20 );
cout << t << ": resized left to " << distr.size() << ", energy = " << lattice.energy( state.first , state.second ) << endl;
}
if( distr[distr.size()-10] > 1E-150 )
{
do_resize( state.first , state.second , distr , state.first.size()+20 );
cout << t << ": resized right to " << distr.size() << ", energy = " << lattice.energy( state.first , state.second ) << endl;
}
t += dt;
}
//]
cout << "final lattice size: " << distr.size() << ", final energy: " << lattice.energy( state.first , state.second ) << endl;
}
|
