4by12

Fun With PDEs – Part 2

Last time I described the first pitfall I encountered when solving a PDE — an inherent instability in the partial derivatives. This time I’ll talk about the second pitfall, which is simpler conceptually, but has wider implications for programming in general.

In my first implementation of the solution I used a simple method to calculate the derivative:

-
This is just the average of the right and left derivatives. I wanted to improve the accuracy of the derivative calculation because of various reasons, so I turned to Savitsky-Golay filters, which were also described in the previous post. Using this method you get a digital filter that you apply to your data. If you’re not familiar with digital filters, you can think of a filter as an array of numbers which are applied to your function f as follows:

-
(I’m probably mixing up some of the signs here, but the intention is clear I’m sure). To calculate a derivative using S-G, you first apply their filter to your data, and then divide by the spatial resolution dx — the distance between adjacent points on the grid. It is worth noting that there is a highly efficient way of doing this calculation using FFT. For more details, see Numerical Recipes.

I implemented S-G and tested it using a simple sin(x) function, and it worked flawlessly. But after I inserted it into my main program, the simulation started spewing out strange results. After some debugging I discovered something very strange: Using S-G, the derivative of a constant function wasn’t zero, and in one case reached . How odd! Debugging further, I found that non-zero derivatives only happened when the initial function f had very large (and constant) values, on the order of , and dx was very small, about .

To continue debugging, I dropped the fancy FFT convolution and switched to straightforward calculation — the one that’s shown in the last equation. This finally revealed the problem: When multiplying each value of the function f by the factor and summing, some of the least significant bits are lost, because the numbers don’t all have the same exponent. So even though the calculation is accurate, the limitation of the computer’s double precision causes a loss of data. When the convolution is done, you’re left with a small value that’s not zero. But then to get the derivative you divide by dx, a very small number, and this shoots that small error through the roof. The smaller dx, the larger the derivative of the constant function! So once again, decreasing dx actually caused a larger error.

What helped me solve this problem was noticing that the S-G coefficients are anti-symmetric, i.e. . Specifically, anti-symmetric coefficients have the same exponent. Therefore I changed the summing order to sum pairs of anti-symmetric factors:

-
This solved the problem: Constant functions now had a zero derivative, because each term in the sum was exactly zero, even in double precision. And the happy consequence was that this fix also solved the weird simulation results I started with.