//? 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"<{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 = "<{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 "<