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.