# MATLAB Fermi-Pasta-Ulam lattice model can't force periodic boundary condition

I am building MATLAB code to implement a symplectic integrator using a second order split-step method for the Fermi-Pasta-Ulam problem. The Hamiltonian is in the form:

H = sum from i=1 to N {0.5*p_i^2 + 0.5*(q_[i+1] - q_i)^2+(a/3)*(q_[i+1] - q_i)^3 + (b/4)*(q_[i+1] - q_i)^4}

I assume a boundary condition: q_[N+1]=q_1 and p_[N+1]=p_1. So I have my oscillators on a string in a circle. The trouble I have is I can't seem to force the periodic boundary condition! Here is my main program and the grad of Hamiltonian in a separate file, I think I have x0 wrong but not sure how to improve on it! The sine wave gets all distorted on the ends. Hamiltonian.m:

clear

numStep = 300;  % number of steps
dt = 0.01;  % time step
t = (1:numStep)*dt; % integration time
a = 0; % set the high order terms to zero for now
b = 0; % set the high order terms to zero for now
N = 32; % number of oscillators on the string
delta = 0.2 % displacement of oscillator from equilibrium position
length = 1; % distance between adjacent oscillators

% This is wrong I think!
x0 = [sin(linspace(0,2*pi,N))*delta+(1:N)*1,zeros(1,N)];

xSymplectic = zeros(2*N,numStep);
xSymplectic(:,1)=x0';

J = [zeros(N),eye(N);-eye(N),zeros(N)]; %symplectic matrix
A = [zeros(N),eye(N);zeros(N),zeros(N)]; % top part of J
B = [zeros(N),zeros(N);-eye(N),zeros(N)]; % bottom part of J

x = x0';
for j=1:numStep
% Symplectic integrator using the split-step method
x = x + A*gradH(x,a,b)*dt/2; % only use kinetic part
x = x + B*gradH(x,a,b)*dt; % only use potential part
x = x + A*gradH(x,a,b)*dt/2; % only use kinetic part
xSymplectic(:,j+1) = x;
InitCond = xSymplectic(1:N,j+1)-(1:N)'*1;
plot(InitCond)
pause(0.1)
end


function y = gradH(x,a,b)
N = length(x)/2;
y(2*N,1) = 0;
% The first half is for the q's (position)
y(1) = -x(2)+2*x(1)-x(N)+N-a*(x(2)-x(1)).^2+a*(x(1)-x(N)+N).^2-b*(x(2)-x(1)).^3+b*(x(1)-x(N)+N).^3;
y(2:N-1) = -x(3:N)+2*x(2:N-1)-x(1:N-2)-a*(x(3:N)-x(2:N-1)).^2+a*(x(2:N-1)-x(1:N-2)).^2-b*(x(3:N)-x(2:N-1)).^3+b*(x(2:N-1)-x(1:N-2)).^3;
y(N) = -x(1)+2*x(N)-N-x(N-1)-a*(x(1)-x(N)+N).^2+a*(x(N)-x(N-1)).^2-b*(x(1)-x(N)+N).^3+b*(x(N)-x(N-1)).^3;
% The second half is for p's (momentum)
y(N+1:2*N) = x(N+1:2*N);
end


Any advice is appreciated, thanks! Aina.

-

I have discovered the error after a whole day of bashing my head: in my sine wave I skip one step at each cycle. This is the code that works: Hamiltonian.m:

clear

numStep = 300;  % number of steps
dt = 0.01;  % time step
t = (1:numStep)*dt; % integration time
a = 0; % set the high order terms to zero for now
b = 0; % set the high order terms to zero for now
N = 32; % number of oscillators on the string
delta = 0.2 % displacement of oscillator from equilibrium position
length = 1; % distance between adjacent oscillators

%%%%%%% This is now correct! %%%%%%%
mesh = linspace(0,2*pi,N+1);
x0 = [sin(mesh(2:end))*delta+(1:N)*1,zeros(1,N)];

xSymplectic = zeros(2*N,numStep);
xSymplectic(:,1)=x0';

J = [zeros(N),eye(N);-eye(N),zeros(N)]; %symplectic matrix
A = [zeros(N),eye(N);zeros(N),zeros(N)]; % top part of J
B = [zeros(N),zeros(N);-eye(N),zeros(N)]; % bottom part of J

x = x0';
for j=1:numStep
% Symplectic integrator using the split-step method
x = x + A*gradH(x,a,b)*dt/2; % only use kinetic part
x = x + B*gradH(x,a,b)*dt; % only use potential part
x = x + A*gradH(x,a,b)*dt/2; % only use kinetic part
xSymplectic(:,j+1) = x;
InitCond = xSymplectic(1:N,j+1)-(1:N)'*1;
plot(InitCond)
pause(0.1)
end

-