//? eqm1.cpp
//? C+- by Ulrich Mutze. Status of work 2021-10-26.
//? Copyright (c) 2021 Ulrich Mutze
//? contact: see contact-info at www.ulrichmutze.de
//?
//? This program is free software: you can redistribute it and/or
//? modify it under the terms of the GNU General Public License as
//? published by the Free Software Foundation, either version 3 of
//? the License, or (at your option) any later version.
//?
//? This program is distributed in the hope that it will be useful,
//? but WITHOUT ANY WARRANTY; without even the implied warranty of
//? MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
//? GNU General Public License for
//? more details.
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
using namespace CpmRoot;
using namespace CpmRootX;
using namespace CpmFunctions;
using namespace CpmSystem;
using namespace CpmArrays;
using namespace CpmGraphics;
using namespace CpmImaging;
using namespace CpmApplication;
using namespace CpmAlgorithms;
using namespace CpmGeo;
using namespace CpmQuantumMechanics;
using namespace std;
/* Expansive Quantum Mechanics.
We consider a particularly simple and transparent quantum mechanical model.
It describes a non-relativistic particle which at the beginning is in free
motion. It then encounters a smooth downhill ramp and gets an acceleration.
Finally it finds itself on a horizontal path. The plateau gets extended whenever
the wave function of the particle reaches the rim of the plateau.
As ever, physical space is idealized as a finite chain of equi-distant
points. There is already a Mathematica implementation of the system
quantumToClassic5.nb (and related files)
My aim is to make a meaningful comparison of this Mathematica-based
implementation and a C++-based implementation.
For the parameters given here the computation time for test1() was 68.4971 s
whereas my Mathematica comparison implementation worked for 676.152 s. Therefore
C++ 10 times faster than Mathematica (for this special case which, however,
looks quite typical for the computations that I presently are interested in.
*/
namespace eqm1{ // expansive quantum mechanics in one dimension.
// I try to use standard C++ (C++20 actually) wherever this is convenient.
// In particular std::complex and std::valarray<>
using N = std::size_t;
using Z = long;
using R = double;
using C = std::complex;
using VR = std::valarray;
using VC = std::valarray;
// concatenation
VC concat(VC const& a, VC const& b){
N na=a.size();
N nb=b.size();
VC c(na+nb);
for (N i=0; i f_;
VR pot_; // gives the values of the potential at the allowed particle
// positions
public:
Bio(N n, R xL, R xU, R m, R v, function f):
n_{n}, xL_{xL}, xU_{xU}, d_{(xU_-xL_)/n_},m_(m),v_(v),
x_{VR(n_)},f_{f},pot_{VR(n_)}
{for(N i=0;i getf()const{ return f_;}
};
class State: public Bio{
VC dat_;
public:
State(Bio const& b):Bio(b),dat_{VC(n_)}{}
State(Bio const& b, VC const& vc):Bio(b),dat_{vc}{}
State make(C const& z)const{ return datChanged(VC(z,n_));}
State make( function g)const;
State random()const;
State datChanged(VC const& s)const;
VC& operator+=(VC const& s){ dat_+=s;return dat_;}
C operator[](N i)const{ return dat_[i];}
VC dat()const{ return dat_;}
State dot()const;
VC dotLean()const;
N dim()const{ return n_;}
R nor()const;
// actually norm, this name resulted in intransparent name conflicts
// so I changed the name without clarifying the conflict
R nor_();
// changes *this into a normalized vector and returns the
// norm of the original
R opNorm()const;
// here *this only brings in the biotope data, dat_ does not enter.
bool mark(Graph& gr,R fac)const;
bool markCf(Graph& gr,R fac)const;
};
R State::nor()const{
R res=0.;
for (N i=0;i g)const{
VC dat(n_);
for (N i=0;i(-1.,1.),default_random_engine());
VC dat(n_);
R sum=0.;
for (N i=0;i=b ? a : b;
}
R inf(R a, R b){
return a<=b ? a : b;
}
bool labile(std::vector const& rL, R res){
N rs=rL.size();
if (rs<=4) return true;
VR tail(3);
tail[0]=rL[rs-3];
tail[1]=rL[rs-2];
tail[2]=rL[rs-1]; // tail holds the 3 last elements of rL
R mean=(tail[0]+tail[1]+tail[2])/3.;
R min=inf(inf(tail[0],tail[1]),tail[2]);
R max=sup(sup(tail[0],tail[1]),tail[2]);
if (res<=max && res>=min && abs(res-mean)/(res+mean)<0.5) return false;
return true;
}
R State::opNorm()const{
State s=random();
R res=-1.;
std::vector rL;
while (labile(rL,res)==true){s=s.dot();res=s.nor_();rL.push_back(res);}
return res;
}
bool State::mark(Graph& gr, R fac)const{
if (fac>0.){
R el=0.5*m_*v_*v_;
VR xVal=x_;
VR potVal(n_);
VR reVal(n_);
VR imagVal(n_);
VR absVal(n_);
for (N i=0;i{BLUE,LIGHTBLUE,GREEN,RED,BLACK});
}
else{
VR xVal=x_;
VR reVal(n_);
VR imagVal(n_);
VR absVal(n_);
for (N i=0;i{GREEN,RED,BLACK});
}
return true;
}
bool State::markCf(Graph& gr, R fac)const{
Z mL=4;
Word loc("State::markCf(Graph& gr, R fac)constt");
CPM_MA
if (fac>0.){
R el=0.5*m_*v_*v_;
VR xVal=x_;
VR potVal(n_);
VR reVal(n_);
VR imagVal(n_);
VR absVal(n_);
for (N i=0;i{BLUE,LIGHTBLUE,GREEN,RED,BLACK});
}
else{
VR xVal=x_;
VR reVal(n_);
VR imagVal(n_);
VR absVal(n_);
for (N i=0;i{GREEN,RED,BLACK});
}
CPM_MZ
return true;
}
////////////////////////////////////////////////////////////////////////////////
class Sys{
R t_;
State psi_;
VC phi_;
R nr_;
R e_;
void ini_(){
nr_=psi_.nor();
e_=real(C(0.,1.)*sp(psi_.dat(),phi_));
}
public:
Sys(State const& st):t_{0.},psi_(st), phi_(psi_.dotLean()){ ini_();}
R e()const{ return e_;}
R nr()const{ return nr_;}
void step_(R dt, N dn, R eps){ // densified asynchronous leap frog (DALF)
R tau=dt*0.5;
R tauh=tau*0.5;
extend_(dn,eps);
t_+=tauh;
psi_+=phi_*tauh;
phi_+=(psi_.dotLean()-phi_)*2.;
psi_+=phi_*tau;
t_+=tau;
phi_+=(psi_.dotLean()-phi_)*2.;
psi_+=phi_*tauh;
t_+=tauh;
ini_();
}
State psi()const{ return psi_;}
VC phi()const{ return phi_;}
void extend_(N dn, R eps);
N dim()const{ return psi_.dim(); }
};
void Sys::extend_(N dn, R eps){
C z0(0.,0.);
Bio& b=psi_; // this is the orthodx method to access the base of State
// from the class Sys
N n=b.n();
R xL=b.xL();
R xU=b.xU();
R d=b.d();
R m=b.m();
R v=b.v();
VR x=b.x();
VR pot=b.pot();
function f=b.getf();
R dl=dn*d;
VC dpsi=VC(z0,dn);
VC psi=psi_.dat();
VC phi=phi_;
bool ready=true;
if (norm(psi_[1])>eps || norm(phi_[1])>eps ){ // squares of the components advisable since sum
// of all squares is always 1 for unit vectors
psi=concat(dpsi,psi);
phi=concat(dpsi,phi);
xL-=dl;
n+=dn;
ready=false;
}
N nCrit=psi_.n()-2;
if (norm(psi_[nCrit])>eps || norm(phi_[nCrit])>eps ){ // squares of the components advisable since sum
// of all squares is always 1 for unit vectors
psi=concat(psi,dpsi);
phi=concat(phi,dpsi);
xU+=dl;
n+=dn;
ready=false;
}
if (!ready){
Bio bMod(n,xL,xU,m,v,f);
psi_=State(bMod,psi);
phi_=phi;
}
// t_ was not modified, it is still upto date.
}
////////////////////////////////////////////////////////////////////////////////////////
R downHill(R x){
R p=0.810747;
if (x<=0.) return 1.;
if (x>=1.) return 0.;
R y= p*(2.*x - 1.)/(x*(x - 1.));
R z= exp(-abs(y));
return y<=0. ? z/(1.+z) : 1./(1.+z);
}
R smoothTransit(R x, R x1, R x2, R y1, R y2){
if (x<=x1) return y1;
if (x >= x2) return y2;
R xr = (x - x1)/(x2 - x1);
return y2 + downHill(xr)*(y1 - y2);
}
C fTest(R x){ return sin(x) + C(0.,1.)*cos(2.*x);}
void test1(){
cout<<"test1() started"<{"potential","energy level",
"wave function:","real part","imaginary part","absolute value"},gr,
V{BLUE,LIGHTBLUE,DARKGRAY,GREEN,RED,BLACK},0.55,-0.15);
}
else{
legend=FrmText(V{"wave function:","real part","imaginary part","absolute value"},gr,
V{DARKGRAY,GREEN,RED,BLACK},0.55,-0.15);
}
gr.setText(legend);
R fac0=1.1;
R fac1=0.6;
R_Vector steps(R2(fac0,fac1),still1);
for (N i=1;i<=still1;++i){
gr.setY(iv0.magDemag(steps[i]));
psi0.markCf(gr,fac);
gr.vis(movie);
cpmwait(0.1,2);
}
cpmwait(t0,2);
Sys sys(psi0);
R e0=sys.e();
cout<<"e0 = "< ev;
std::vector nrv;
ev.push_back(0.);
nrv.push_back(0.);
R opN=st.opNorm();
R dtCrit=2./opN;
R dt=dtCrit*dtRel;
cout<<"opN = "<{DARKGRAY},0.6,-0.15);
gr.setText(legend2);
gr.mark(rev);
for(N i=1;i<=still2;++i){gr.vis(movie);}
cpmwait(t2,2);
FrmText legend3(V{"relative norm error"},gr,V{DARKGRAY},0.6,-0.15);
gr.setText(legend3);
gr.mark(rnrv);
for(N i=1;i<=still2;++i){gr.vis(movie);cpmwait(0.2,2);}
cpmwait(t3,2);
cout<<"time for test2() was "<=2){
cout<<"args_[2] = "<=3){
cout<<"args_[3] = "<=4){
cout<<"args_[4] = "<