CoCalc Public FilesMIWsmi1Dv2.sagews
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)=[x_1(t)...x_{N*J*D}(t)]$ is determined in general in MIW by system of $N*J*D$ coupled non-linear second order partial differential equations, where $N$=number of worlds, $J$=number of particles, $D$=number of dimensions. $m_k*\ddot{x}_{kn}(t) = f_k(x_n(t)) + r_{kN}(x_n(t);Xp(t))$ where $k$ indexes the particle and $n$ the world.

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

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

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

where I have suppressed the non-linear functional dependencies of $F$ and $R$, and incorporated the mass scaling into the new function $Fc$ $Rc$ 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 ($x$ and $\dot{x}$) to be named $X$.

For our 1 D case, there will be $N$ equations. $X$ has $2N$ state variables (position, velocity) In a discrete state space representation, the system is described as $2N$ first order linear state transition equations, and is representable as a $2N$ 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)=A*X(n)+F(x_n)+R(x_n,X)$ with $A=expm(Ac*T)$; where $Ac$ is the continuous time state propagation matrix, $T$ 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 $T$ approaches 0, for any finite value of $T$, 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 ();

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]);

%sage salvus.file("MIW1D.mp4")