Ganesh Swami

One of the many things I learnt from my nonlinear physics project is the pseudospectral method. I’m going to try explaining this method the way I understand it.

Any function can be expressed as an infinite sum of the set of basis functions:

For the purpose of an example, if we take the polynomial basis of order 2, then the function can be approximated as

The coefficients can be solved for by building a system of equations. The three unknowns require three points on the interval where we know the exact value of the function. Let these three points be x0, x1 and x2 and the value of the function at these points be b0, b1 and b3.

The system of equations is now,

written compactly as

The system can be solved by inverting the matrix A, which is of complexity . The function p(x) is an interpolating polynomial in the interval of interest. As you can see, using a higher order polynomial quickly becomes prohibitively expensive.

Instead of picking the polynomial basis, you can pick any other basis. The Chebyshev basis is commonly used. A friend of mine picked the Fourier-Bessel basis. Depends on the application.

Fourier basis

The Fourier basis is fundamental to linear systems theory. Linear operators become diagonal in the Fourier basis. This fact comes handy in numerical schemes such as exponential integrators, where you’ll have to take the exponential of a matrix.

The second advantage falls out of the uncertainty principle. Increasing the number of modes in the Fourier domain is equivalent to shrinking the spacing between steps in the space domain. This leads to an exponential convergence in the number of points required as opposed to the usual polynomial convergence with finite differences.

The third key point in the use of the Fourier basis is that the FFT lets us do the equivalent of a matrix inversion in . This is very fast.

Summary

In summary, three random facts come together to form a beautiful theory:

• The FFT is much faster than a matrix inversion.
• Linear operators are diagonal in the Fourier basis.
• Exponential convergence of spectral methods.

Hope you enjoyed reading this post as much as I enjoyed writing it.

Dynamical systems following Hamilton mechanics can be formulated as

where M is the mass and q is the state vector. These equations arise in pretty much every single physics engine (astrophysics and molecular dynamics for example.) Traditionally, solvers such as the forward/backward Euler have been used, but these solvers do not respect the manifold of the configuration system. This led to the development of solvers such as the Verlet (and Velocity Verlet) which are symplectic in nature and follow the geometry of the problem. The Velocity Verlet specifically is the workhorse of molecular dynamics.

By geometry, I mean the invariants that arise due to symmetry in the system. For example, rotational and translational symmetry in the system give rise to conservation of angular and linear momentum. Above, the equations of motion for a simple pendulum were integrated with four different solvers. The symplectic nature of the problem does not fall out naturally from the solvers because these “local” methods still look at differential changes in momenta and position.

Variational methods, on the other hand, directly deal with equations arising out of Hamilton’s action principle. The Lagrangian is defined as

and the action functional is the integral of the Lagrangian along a path q(t).

First order variations to compute the stationary action leads to the Euler-Lagrange differential equation. A similar derivation can be done for discrete variables yielding the Discrete Euler Lagrange (DEL) equation.

This is very attractive because physical invariants that arise from the variational principle are guaranteed to be maintained in the discrete situation as well. This is also true for constraints on the system: constraint versions of Verlet and Velocity Verlet, SHAKE and RATTLE are natural in the DEL equation setting.

The only roadblock to the rapid adoption of these algorithms is the computational expense of solving the DEL. Each integration timestep requires the solution of a set of implicit nonlinear equations (usually by the use of Newton’s method.) This is a problem for molecular dynamics where you are trying to run the simulation for a few nanoseconds and your timestep is in femtoseconds.

I was almost going to choose this topic for my nonlinear physics course project. Though this is very interesting and exciting, I couldn’t find enough content to develop into a project.

Granular media, like sand, rice or cereal, also happen to form interesting patterns when vibrated on a membrane. This was first discovered and published in Nature by Umbanhowar, Melo & Swinney from the Center for Nonlinear Dynamics at the University of Texas at Austin. This was the same team that had created the very popular Faraday waves video posted on YouTube.

The physical phenomena is similar to Faraday waves seen in liquids almost a hundred and thirty years ago. The key difference is that while liquid faraday waves can be derived from a continnum model such as the Navier-Stokes equations, no such model exists for granular media as far as my limited knowledge goes. Therefore, the analytical derivations for the purposes of a project is limited.

The phase-plot is in terms of a dimensionless amplitude given by

where A is the driving amplitude, g is the acceleration due to gravity, and f is the driving frequency. Some of these have been experimentally verified by means of a molecular dynamics simulation. Once again, the computational part wasn’t sufficient to fulfill the requirements because the equations are “simple” non-linear coupled ordinary differential equations. Whereas systems we have studied in class were described by partial differential equations.

The term “oscillon” was coined by the paper mentioned above. In recent years, there has been a lot of research into these structures. The picture above is oscillons in a nonlinear field model and the author has some really cool videos on his site.

I spent most of last week looking for an idea for my nonlinear physics course project. We were given a few pointers to topics that could become potential projects: various kinds of pattern formation in Fluid Mechanics, Complex fluids, Biological systems, Chemical systems, Optics, Combustion and Hydrodynamics.

I came up with four topics (all four of them weren’t on the initial list of things suggested.) Somebody else in the class had exactly the same proposal as mine, including references, so that was out of the picture. One turned out to be computationally infeasible (couple of hours for a single step of a Finite Element Method), the other turned out to be extremely hard and outside the scope of the course.

I finally settled on “Global feedback methods for the subcritical Complex Ginzburg-Landau Equation.” If that didn’t make any sense, I’ve tried to explain what each word means in my proposal.

(picture courtesy of Shankar Venkataramani.)

Real life has been kicking my ass lately. My brain has been working over-time, so I’m taking some time off to relax. I’m fortunate to have a lot of friends with birthdays this week, so that helps. Anyways, here’s an update on what I’ve been upto…

Parabolic excitation

Working on a homework problem for my non-linear physics class has been a huge timesink. This week, we were asked to solve the general low viscosity Mathieu equation

for a parabolic excitation signal. My professor, Dr. Bechhoefer had published a paper to describe parametrically excited surface waves with delta and triangle excitation signals. We were asked to extend this. At the end, you get a complicated implicit equation for the threshold condition. I had to learn how to solve implicit equations numerically in Matlab (numerically because Maple couldn’t do it analytically.)

Geometry of Diffusion Tensors

As I’m working a lot with Diffusion Tensors, I wanted to get a better feel for the computations on them. I found this paper(pdf) by Dr. Joshi titled “Principal Geodesic Analysis on Symmetric Spaces: Statistics of Diffusion Tensors” really helpful. Ofcourse, I couldn’t understand some of the mathematical terms, so it’s more reading for me. I also realized that not all graduate students are on top of their game — they sometimes can mislead you. It’s better to mislead yourself than to be misled.

Geodesic Shooting

As part of my weekly reading, a friend and I are trying to understand this paper titled “Geodesic Shooting for Computational Anatomy.” Amidst all the complicated math, the basic idea they are trying to show is that flows on the deformation diffeomorphism conserve momentum. And because of this conservation law, you only need the velocity at t=0 to completely determine the flow. In the grand scheme of things, you do not and cannot average images just by averaging pixel/voxel intensities. This is because these images do not form a vector space and simple linear averaging does not respect the curved aspect of the space. The space of the velocity vectors form a vector space, and we can average those instead.

MRI scan

I also had my first MRI scan last week. I had to go all the way to the UBC hospital for that. A graduate friend of mine is doing some research on parameter optimization on the MR machine, but with about 400 controls you can tune, this is a very difficult problem. It was an interesting experience for me, though I slept through most of it. To take your mind off the heavy beating (and for claustrophobic people I guess,) they give you headphones. I wondered for a long time how electronic headphones worked in presence of a strong magnetic field (MR machines have a magnetic field strength of 1.5 to 3.0 tesla, in contrast to the earth’s field of 30 to 60 microteslas.) I’ll leave that question unanswered because it makes me feel too dumb.

That’s all for now. I’ll save the rest for later.