2
$\begingroup$

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.

1 Answers 1

1

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