Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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

gradH.m:

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.

share|improve this question

1 Answer 1

up vote 1 down vote accepted

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
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.