CoCalc Public FilesMIWsmi1Dv2.sagewsOpen with one click!
Authors: David Cyganski, Bill Page
Views : 51
Description: Implementation of Wiseman many interacting worlds simulation for 1D, 1 Particle systems
Compute Environment: Ubuntu 18.04 (Deprecated)

###MIWsim1D###

Author: D Cyganski

Created: 2015-05-17

%auto %default_mode octave

####Implementation of Wiseman many interacting worlds simulation for 1D, 1 Particle systems####

  • v1 - demo of packet spreading and Fresnel diffraction of 1D packet with spreading
  • v2 - single particle with gaussian distribution and nonzero mean velocity, colliding with a potential barrier
  • v3 - single particle in two gaussian superposition with velocities, colliding with no interacting fields

The future track of the many worlds position state Xp(t)=[x1(t)...xNJD(t)]Xp(t)=[x_1(t)...x_{N*J*D}(t)] is determined in general in MIW by system of NJDN*J*D coupled non-linear second order partial differential equations, where NN=number of worlds, JJ=number of particles, DD=number of dimensions. mkx¨kn(t)=fk(xn(t))+rkN(xn(t);Xp(t)) m_k*\ddot{x}_{kn}(t) = f_k(x_n(t)) + r_{kN}(x_n(t);Xp(t)) where kk indexes the particle and nn the world.

Such a system of NJDN*J*D second order non-linear PDEs can always be written as a system of of 2NJD2*N*J*D first order linear PDEs:

which in turn can be put into matrix state-system form for state vector X=[x, y]'

X˙(t)=AcX(t)+Fc(t)+Rc(t) \dot{X}(t) = Ac*X(t)+Fc(t)+Rc(t)

where I have suppressed the non-linear functional dependencies of FF and RR, and incorporated the mass scaling into the new function FcFc RcRc and stacked all the particles and worlds into the single vector

We will now discretize this matrix differential equation and introduce the discrete time state vector of positions and velocities (xx and x˙\dot{x}) to be named XX.

For our 1 D case, there will be NN equations. XX has 2N2N state variables (position, velocity) In a discrete state space representation, the system is described as 2N2N first order linear state transition equations, and is representable as a 2N2N dimenstional matrix state space equation.

clear all; %Generate initial distribution of particle states, including position,x, and velocity, xdot % X=[x1,x2,....xN,x1dot,x2dot,...xNdot]' % X=InitialGaussian1D(N,0,1,1); %(worlds, center, width, velocity) %For a simple test, distribute 101 values equally over a window and set all velocities to zero %xvalues=[-1:.05:1]'; %Prob2 %xvalues=[-1:.025:1]'; %Prob3 xvalues=[-0.5:.025:0.5]'; %Prob4 %xuniform=[-0.95:0.1:0.95]'; %For the simple test
%octave %xuniform=[-0.95:0.025:0.95]'; %Use inverse error function to generate a Gaussian distribution of trajectories %xvalues=erfcinv(xuniform); %xvalues=erfinv(xuniform); N=size(xvalues,1); %Number of worlds X = [xvalues; zeros(N,1)]; %zero initial velocities %X = [xvalues; ones(N,1)]; %non-zero initial velocities
h = figure('visible', 'off'); %Plot the input trajectory density Pin=(N*(xvalues(2:N,1)-xvalues(1:(N-1),1))).^(-1); plot(xvalues(2:N,1),Pin); %xlim([-1,1]); ylim([0,1]); saveas(h,"figure1.png");
%sage salvus.file("figure1.png")

Now use the appropriate state propagation matrix to step this system forward in time. X(n+1)=AX(n)+F(xn)+R(xn,X) X(n+1)=A*X(n)+F(x_n)+R(x_n,X) with A=expm(AcT)A=expm(Ac*T); where AcAc is the continuous time state propagation matrix, TT is the time interval between time steps, and expm is octave's matrix exponentiation function.

Note: The HDW paper does not explicitly use the matrix exponentiation based conversion of the continuous time to discrete time conversion. While a naive conversion produces results that converge to the same solution as TT approaches 0, for any finite value of TT, the exponentiation approach will produce solutions closer to that of the un-discretized system of equations.

I discovered that when I applied the exponentiation to the Ac matrix, thanks to the fact that the differential equations involved did not involve any state variables in its linear part, the outcome of expm(Ac) yields the most simple (naive) summation equation. Since the matrix A that results is extremely sparse, it is computationally advantageous to simply use direct for loop implementation than matrix multiplication.

T=24000; %number of time steps to compute dT=0.08e-2; %time interval of a single step m =1; %mass of particle hbar=.1; %Planck constant %.02 works but nonguassian after 500 iterations %Reserve recording space Xrec=zeros(N,T); %Define a force field at from x=4.9 to 5.1 which repells the oncoming particles %F=inline("0"); %Simply allow the initial pulse to evolve %F=inline("ifelse((x>4.7)&&(x<5.3), -0.8, 0)"); %For the simple test F=inline("ifelse((x>4.7)&&(x<5.3), ifelse(x<5,-1.5,1.5), 0)"); %Triangular potential
t0 = clock (); %Perform time updates for t=1:T %Update the y (==xdot) values for n=(N+1):(2*N) R=MIWforce(N,n-N,X(1:N),m,hbar); %Equations 24,25 of HEW paper X(n)=X(n)+(F(X(n-N))+R)*dT; endfor; %Update the x values for n=1:N X(n)=X(n)+X(n+N)*dT; %Record the position state for plotting Xrec(:,t)=X(1:N,1); endfor; endfor; t elapsed_time = etime (clock (), t0)
t = 24000 elapsed_time = 694.00
trec=[1:T]'*dT; h = figure('visible', 'off'); plot(trec,Xrec); saveas(h,"figure2.png"); %ylim([-2,2]);
%sage salvus.file("figure2.png")
h = figure('visible', 'off'); %Plot the density at the end of the simulation Dist=((Xrec(2:N,T)-Xrec(1:(N-1),T))); Prob=(1/N)*Dist.^(-1); plot(Xrec(2:N,T),Prob); saveas(h,"figure3.png"); %xlim([-1.5,1.5]); ylim([0,1]);
%sage salvus.file("figure3.png")
h = figure('visible', 'off'); %Plot the density at 3/4 through the simulation Tsamp=3*T/4; Dist=((Xrec(2:N,Tsamp)-Xrec(1:(N-1),Tsamp))); Prob=(1/N)*Dist.^(-1);
Inline plot failed, consider trying another graphics toolkit
error: print: no axes object in figure to print
error: called from
    _make_figures>safe_print at line 122 column 7
    _make_figures at line 44 column 13

plot(Xrec(2:N,Tsamp),Prob); saveas(h,"figure4.png"); %xlim([-1.5,1.5]); ylim([0,1]);
%sage salvus.file("figure4.png")
h = figure('visible', 'off'); Dist=((Xrec(2:N,T)-Xrec(1:(N-1),T))); Prob=(1/N)*Dist.^(-1); plot(xvalues(2:N,1),Pin); hold on Dist=((Xrec(2:N,T)-Xrec(1:(N-1),T))); Prob=(1/N)*Dist.^(-1); plot(Xrec(2:N,T),Prob,"r"); saveas(h,"figure5.png"); hold off
%sage salvus.file("figure5.png")
figure("visible","off"); moviename="MIW1D.mp4"; movie("init",moviename) for Tsamp=1:10:T Dist=((Xrec(2:N,Tsamp)-Xrec(1:(N-1),Tsamp))); Prob=Dist.^(-1); plot(Xrec(2:N,Tsamp),Prob); xlim([-20,20]); ylim([0,20]); movie("add",moviename); endfor movie("close",moviename,40); %last arg is frame rate
GPL Ghostscript 9.18: **** Could not open the file 'MIW1D.d/000237.png'. GPL Ghostscript 9.18: Unrecoverable error, exit code 1 error: save: unable to open output file 'MIW1D.d/+frame-number+' error: called from movie at line 82 column 7 error: load: unable to find file MIW1D.d/+frame-number+ error: called from movie at line 102 column 9
%sage salvus.file("MIW1D.mp4")
Error in lines 1-1 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1009, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 838, in file file_uuid = self._conn.send_file(filename) File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 167, in send_file f = open(filename, 'rb') IOError: [Errno 2] No such file or directory: u'MIW1D.mp4'