Fun With PDEs – Part 2
September 17, 2007 – 6:42 amLast 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:
![\frac{\partial f}{\partial x}[j] = \frac{f[j+1]-f[j-1]}{2 \, dx} \frac{\partial f}{\partial x}[j] = \frac{f[j+1]-f[j-1]}{2 \, dx}](/latexrender/pictures/5a99aae890fe9ef5482881ed643463ad.png)
-
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:
![F[j] = \sum_{n=-l}^{l} f[j-n] * c_n F[j] = \sum_{n=-l}^{l} f[j-n] * c_n](/latexrender/pictures/9407c374a365816609eb54c168b72267.png)
-
(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:
![F[j] = \sum_{n=1}^{l} (f[j-n] \cdot c_n + f[j+n] \cdot c_{-n}) + f[j] \cdot c_0 F[j] = \sum_{n=1}^{l} (f[j-n] \cdot c_n + f[j+n] \cdot c_{-n}) + f[j] \cdot c_0](/latexrender/pictures/030368e3f1eca0972f29ac06368059f6.png)
-
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.
prisoners are planning an escape. They bribed a guard to let them all meet once. Each prisoner was given a task to complete, such as digging or stealing some key. Once all the tasks are complete, the prisoners will notify the corrupt guard who will turn a blind eye while they escape. The problem is, the prisoners can’t talk to each other because they are all in solitary confinement.
. 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:
, 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:





![f[j] = f[j] + \frac { \partial f } { \partial t } [j] * dt f[j] = f[j] + \frac { \partial f } { \partial t } [j] * dt](/latexrender/pictures/3de923578d6c440441bbc686d59bc5bf.png)
![f[j] = \frac{f[j-1] + f[j+1]}{2} + \frac { \partial f } { \partial t } [j] \; dt f[j] = \frac{f[j-1] + f[j+1]}{2} + \frac { \partial f } { \partial t } [j] \; dt](/latexrender/pictures/08f7b12d31b3ef545da95e66ae44b89a.png)
be a real function such that
has a limit at each point. Does
be a function such that

is continuous.
. Intuitively,
is
,
will be close to
will also be close to
. Then there exists
such that


such that
. Let’s count the points of
and define the set:
is finite.
is compact we can find a series
such that
and
for some
. But for large enough 
must be finite.
. Then from the lemma we have that A is an enumerable set. A also contains all the points for whom
.



is enumerable, so
, so there exists a point
.
Subscribe by email