123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141 |
- /*
- * 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;
- }
|