# Simulation, Dynamics, and the FPU Problem

<– Back to main page…

I recently wandered by Powell’s Technical Books (one of the many perks of working in Portland is having this nearby) and left with the book “The Genesis of Simulation in Dynamics: Pursuing the Fermi-Pasta-Ulam Problem”. I think it’s a pretty interesting book relative to others that I’ve picked up on dynamics and chaotic systems. The book revolves around one problem and follows its development since it was first posed in the 1950s as the field of dynamics unfolded around it. I figured I’d give a quick description of the problem that the book discusses here and possibly interest others to grab the book.

The model that the book discusses is known as the Fermi-Pasta-Ulam (or FPU) problem. It is quite simple. Consider a 1D string that is represented as a sequence of line segments connected at n vertices. The connections between vertices are basically springs that obey Hooke’s law. In the model, the focus is on the displacement of these vertices over time. For the case where all vertices are equally spaced, the string doesn’t move at all since it is in an equilibrium state. When they are not uniformly spaced, interesting things start happening.

This video shows the case where the vertices are displaced based on a sine wave.

(Note: the slowing and speeding up of the video is artificial and not part of the model, but is due to variations in compute time per frame on my computer since I recorded this as a screen capture).

The model is defined by tracking the second derivative of the position of the vertices. The unperturbed linear model is described as:

$\ddot{q}_j = (q_{j+1}-q_j) + (q_{j-1}-q_j)$

Adding in a perturbation that is nonlinear is what introduces the interesting behavior:

$\ddot{q}_j = (q_{j+1}-q_j) + (q_{j-1}-q_j) + \alpha \left[(q_{j+1}-q_j)^2 - (q_j - q_{j-1})^2\right]$

The parameter alpha allows us to control the amount that this nonlinear term contributes to the evolution of the system. In my example movies, I let alpha be somewhere around 0.5. In Matlab, this is easily written as:

function fpu(n,alpha,s,dt,iters)
% initial positions for n points on the line
% and two additional points for the fixed ends
q = (0:n+1)./(n+1);
qrest = (0:n+1)./(n+1);

% velocities start as zero
qvel = zeros(1,n+2);

% perturbation added with a sinusoid
q = q + 0.1.*sin(s.*q .* pi);

% make sure ends are fixed at 0 and 1
q(n+2) = 1; q(1) = 0;

for i=1:iters
% first term of q"
qterm1 = q(3:end)-q(2:end-1);
% second term of q"
qterm2 = q(1:end-2)-q(2:end-1);

% qacc is q".  note that (a-b)^2 = (-(b-a))^2
qacc = qterm1 + qterm2 + alpha .* (qterm1.^2 - qterm2.^2);

% velocity is velocity + acc
qvel(2:end-1) = qvel(2:end-1) + qacc;

% position is updated by velocity * time step
q(2:end-1) = q(2:end-1) + qvel(2:end-1).*dt;
end
end

Adding a few strategic plot commands at the end of that loop lets you visualize it and see the plot I showed above.

The goal of the model originally was to look at how the energy in the system moved between the different vibrational modes of the system. As can be seen pretty clearly in the video above, over time, the single low frequency sinusoid that the system starts with evolves into a more complex motion where higher frequency modes gain energy and contribute more to the state of the system.

We can look at the energy in each mode by writing the Hamiltonian of the system as:

$H_0(a,\dot{a}) = \frac{1}{2}\sum_{k=1}^n(\dot{a}_k^2 + \omega_k^2a_k^2)$

where:

$\omega_k = 2 sin \left(\frac{k \pi}{2(n+1)}\right)$

So, if we compute each component of the sum and look at them individually (instead of looking at the overall energy of the system represented by $$H_0$$), we can see how the different modes contribute to the overall energy of the system. The code above can be easily modified to add this computation after the update to the q vector:

    for k=1:n
a(i,k) = (sqrt(2/(n+1))) * sum(q(2:end-1).*sin((1:n)*k*pi/(n+1)));
ap(i,k) = (sqrt(2/(n+1))) * ...
sum((j/n+1).*pi.*q(2:end-1).*cos((1:n)*k*pi/(n+1)));
end
h(i,:) = 0.5 .* (ap(i,:).^2 + omega.^2 .* a(i,:).^2);

Here is a look at the energy of modes 2-8 for the first 30000 iterations of the Matlab code: