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
|
/*
* gauss_packet.cpp
*
* Schroedinger equation with potential barrier and periodic boundary conditions
* Initial Gauss packet moving to the right
*
* pipe output into gnuplot to see animation
*
* Implementation of Hamilton operator via MTL library
*
* Copyright 2011-2013 Mario Mulansky
* Copyright 2011-2012 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 <complex>
#include <boost/numeric/odeint.hpp>
#include <boost/numeric/odeint/external/mtl4/mtl4.hpp>
#include <boost/numeric/mtl/mtl.hpp>
using namespace std;
using namespace boost::numeric::odeint;
typedef mtl::dense_vector< complex< double > > state_type;
struct hamiltonian {
typedef mtl::compressed2D< complex< double > > matrix_type;
matrix_type m_H;
hamiltonian( const int N ) : m_H( N , N )
{
// constructor with zero potential
m_H = 0.0;
initialize_kinetic_term();
}
//template< mtl::compressed2D< double > >
hamiltonian( mtl::compressed2D< double > &V ) : m_H( num_rows( V ) , num_cols( V ) )
{
// use potential V in hamiltonian
m_H = complex<double>( 0.0 , -1.0 ) * V;
initialize_kinetic_term();
}
void initialize_kinetic_term( )
{
const int N = num_rows( m_H );
mtl::matrix::inserter< matrix_type , mtl::update_plus< complex<double> > > ins( m_H );
const double z = 1.0;
// fill diagonal and upper and lower diagonal
for( int i = 0 ; i<N ; ++i )
{
ins[ i ][ (i+1) % N ] << complex< double >( 0.0 , -z );
ins[ i ][ i ] << complex< double >( 0.0 , z );
ins[ (i+1) % N ][ i ] << complex< double >( 0.0 , -z );
}
}
void operator()( const state_type &psi , state_type &dpsidt , const double t )
{
dpsidt = m_H * psi;
}
};
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 << "p [0:" << mtl::size(x) << "][0:0.02] '-'" << endl;
for( size_t i=0 ; i<mtl::size(x) ; ++i )
{
cout << i << "\t" << norm(x[i]) << "\n";
}
cout << "e" << endl;
}
++m_count;
}
};
static const int N = 1024;
static const int N0 = 256;
static const double sigma0 = 20;
static const double k0 = -1.0;
int main( int argc , char** argv )
{
state_type x( N , 0.0 );
// initialize gauss packet with nonzero velocity
for( int i=0 ; i<N ; ++i )
{
x[i] = exp( -(i-N0)*(i-N0) / ( 4.0*sigma0*sigma0 ) ) * exp( complex< double >( 0.0 , k0*i ) );
//x[i] += 2.0*exp( -(i+N0-N)*(i+N0-N) / ( 4.0*sigma0*sigma0 ) ) * exp( complex< double >( 0.0 , -k0*i ) );
}
x /= mtl::two_norm( x );
typedef runge_kutta4< state_type > stepper;
// create potential barrier
mtl::compressed2D< double > V( N , N );
V = 0.0;
{
mtl::matrix::inserter< mtl::compressed2D< double > > ins( V );
for( int i=0 ; i<N ; ++i )
{
//ins[i][i] << 1E-4*(i-N/2)*(i-N/2);
if( i < N/2 )
ins[ i ][ i ] << 0.0 ;
else
ins[ i ][ i ] << 1.0 ;
}
}
// perform integration, output can be piped to gnuplot
integrate_const( stepper() , hamiltonian( V ) , x , 0.0 , 1000.0 , 0.1 , write_for_gnuplot( 10 ) );
clog << "Norm: " << mtl::two_norm( x ) << endl;
return 0;
}
|
