4by12

Fun with PDEs

May 1, 2007 – 4:58 am

I just finished working on a numerical simulation of a set of partial differential equations (PDE). I developed these equations for a physics research project I’m involved in. The equations did not seem to be solvable analytically, so I had to do it numerically. This was my first attempt at solving a PDE, and writing the simulation turned out to be much more involved than with ordinary differential equations. Here are a couple of interesting pitfalls I encountered.

PDE Primer

In case you’re not familiar with the terminology, I’ll first explain what a PDE is. A simple equation contains one or more unknowns which represent numbers. For example: . An ordinary differential equation (ODE) is similar, except the unknown is a function rather than a number. Such an equation involves derivatives of the function. Here is an example:

-
The solution of this particular equation is , where can be any number. Finally, a partial differential equation (PDE) involves a function that has two or more parameters, and includes partial derivatives of this function. For example, the following equation describes waves propagating through a medium:

-
PDEs are very important in physics. In fact, many of the basic laws of nature are described as PDEs. Examples include Maxwell’s equations, Shcrodinger’s equation, and Einstein’s field equations.

On to the simulation!

Pitfall 1: Exploding Waves

So, I built my model for the problem, derived the equations, and was ready to solve them. By ’solving’ I mean that I start out with the known function at time t=0, and I want to find out what that function is at a later time. My function initially looked like this:

Some thousands of time-steps later, it evolved into this:

So far so good, but then it completely exploded:

Going back a bit in time, I was able to trace the beginning of this explosion:

And zooming in on the ‘wavy’ part:

It looked as though waves were forming on my function, and then ‘exploding’.

Inherent Instabilities

I was certain I had a bug, but I couldn’t find it. While debugging, at one point I decreased the spatial resolution — using less points per unit of space to describe the function… and the problem was gone! So, decreasing the accuracy of my solution actually solved the instability… That was very weird.

Mentioning this to a Ph.D student at the lab, he said this problem sounded familiar to him. And as it turns out, this is a universal problem with PDEs: If the time step is too large compared with the spatial resolution, the amplitude of small waves with short wavelengths quickly increases with time until they dominate the solution. This is due to the way numerical derivatives are calculated. The difficulty here is that the time step needs to be incredibly small, making calculation unfeasible. For some equation, the situation is even worse, as they are unstable for any time step, no matter how small.

For simple PDEs, it is very easy to see this effect by taking the function f to be a wave, and watching what happens to the amplitude over time. You can see a derivation of this result here. For a more in-depth discussion, Numerical Recipes is your friend. This method of analyzing equations is called von Neumann stability analysis.

In Comes Lax

Okay, so I found out not alone, but what can be done to solve this problem? The first thing I tried was to calculate the derivative more accurately. There is a method called Savitzky-Golay, where you fit a polynomial to your function at each point, and calculate the polynomial’s derivative at that point. The brilliant thing is that this whole operation (fit + derive) can be done using a single convolution, which costs a meager O(n log n) of processing time.

So I implemented S-G, only to discover it doesn’t solve the problem. More on that in a future post.

As it turns out, there is an incredibly simple solution due to Lax, which says the following. When advancing the function value to the next time step, you do something like this for each position:

-
The Lax method says that the f[j] at the right-hand side should be replaced by an average of it’s neighboring cells:

-
And that’s it! This replacement causes a numerical diffusion that ’sedates’ the unruly waves, causing them to decay instead of explode. The time step used in the simulation still needs to be below some value, but now it decreases linearly with the spatial distance dx, which is much better than before. So Lax saved the day — and that was the end of my first pitfall. This is getting to be quite a long post, so I’ll describe the second problem in another post. Cheers!