4by12

Fun With PDEs – Part 2

September 17, 2007 – 6:42 am

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:

$\frac{\partial f}{\partial x}[j] = \frac{f[j+1]-f[j-1]}{2 \, dx}$

-
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 $c_n, \, n=-l,\cdots,l$ which are applied to your function f as follows:

$F[j] = \sum_{n=-l}^{l} f[j-n] * c_n$

-
(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 $10^{44}$ . 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 $10^{22}$ , and dx was very small, about $10^{-10}$ .

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 c_n 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. $c_{-n} = -c_n$ . 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$

-
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.

Posted in Mathematics, Physics | No Comments »

No More Time Waste

August 30, 2007 – 3:57 am

Lately I’ve been spending way too much time on tech news sites like reddit. Whenever I had some free time and thought about doing something productive like reading a book, I would end up gravitating toward one of these time suckers instead. Hours would pass. From time to time I’d learn something new, but most of it was spent on mind-numbing entertainment like lolcat. (my personal favorite btw).
Read the rest of this entry »

Posted in Personal | 5 Comments »

August 19, 2007 – 8:12 pm

It’s not that I’m so smart, it’s just that I stay with problems longer.

– Albert Einstein

Posted in Quotes | No Comments »

Prisoners’ Escape – The Solution

August 13, 2007 – 12:30 pm

(you can find the original riddle here)

Let’s begin with a draft solution, one that doesn’t work but will start us in the right direction. And to keep things simple, let’s say all the tasks are done the minute the prisoners leave the meeting, so all they have left is to notify someone that they’re done.

Prisoner 1 is responsible for telling the guard once all the tasks are complete. We’ll call him the Leader. He keeps a counter called ‘num_done‘ which starts at 1 (denoting that prisoner 1 himself is done). When the leader enters the room, this is what he does:

if light is on:
    num_done = num_done + 1
    turn light off
else:
    do nothing
-

When any other prisoner enters the room, this is what he does:

if light is off and not turned_light_on:
    turn light on
    turned_light_on = true
else:
    do nothing
-

And each such prisoner maintains a flag called turned_light_on, initialized to ‘false’.

At the end of this, num_done should show the number of prisoners whose task is complete. And when num_done = n, prisoner 1 will notify the guard.

The problem is that this algorithm has an off-by-one error. If the light was originally on before any prisoner entered the room, then when the leader enters the room he will count that on-light as one prisoner. Therefore he’ll reach a count of n while only (n-1) prisoners reported in.

To solve this problem we need to modify the algorithm a bit, to handle the possible initial states of the light switch. We want to give the leader a chance to see the initial state before anyone else touches it. So the other prisoners start by doing nothing, until they get some sort of signal from the leader. This signal will be ‘flicking the switch’: They first want to see the light turned on, and then, when the time is right, they will turn it on themselves.

So now, each prisoner other than the leader maintains another flag called saw_light_on, which starts as ‘false’. The algorithm is:

if light is on:
    saw_light_on = true
else if saw_light_on and not turned_light_on:
    // light is OFF
    turn light on
    turned_light_on = true
else:
    do nothing
-

Next, the leader maintains an extra flag: saw_initial_state initialized to ‘false’. Her algorithm is:

if saw_initial_state:
	if light is on:
	    turn light off
	    num_done = num_done + 1
	else
	    // flick the switch to signal the others
	    turn light on
	    num_done = num_done - 1
else:
	// This is the first time the leader enters the room
	turn light off (if it's on)
	saw_initial_state = true
-

Note that the leader needs to continue flicking the switch forever. Otherwise, latecomers will never turn the light on. Also note that the leader must discount for his own flicking by decrementing num_done. Finally, as before, when num_done = n, the leader notifies the guard and the prisoners escape. There you have it!

Posted in Riddles | No Comments »

Prisoners’ Escape

August 1, 2007 – 9:53 pm

In an effort to make this blog more interactive, I’m switching the riddles to a new format. I’ll first post only the riddle, without the solution, so you guys can have a crack it. After a week or so, I’ll post the solution. Here goes…

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.

Once a day, a prisoner chosen at random enters a storage room that has a single light bulb in it. The prisoner can either flick the switch (turning the light bulb on or off), or leave it as it is. At the time of the initial meeting, as they plan their strategy, the state of the switch is unknown.

Using only switch flicking, how can the prisoners coordinate their escape? Specifically, one of the prisoners must know at some point that all the tasks are done, so he can notify the guard and lead them all to freedom. If the guard is notified before all tasks are complete, the plan will fail and the prisoners will be caught.

(ready for the solution? here it is)

Posted in Riddles | No Comments »

aMule Hacking

May 31, 2007 – 3:19 am

aMule is an open source eMule-like client. It enforces a strict upload/download ratio at the client side: You set the maximum upload rate, and your down rate is limited accordingly. For example, if you allow aMule to use 5 kB/s for upload, your maximum download rate is 20 kB/s. 9 kB/s upload gives you 36 kB/s. 10 kB/s gives you unlimited download bandwidth.

My problem is (or rather, was ) this: 36 kB/s for download is very little compared with my ADSL line’s capacity of 180 kB/s. But to get the full speed I need to allow aMule to use 10 kB/s out of about 12 kB/s of upload capacity, which makes web surfing extremely painful (ping google.com gives > 2 secs). I don’t mind sharing my bandwidth, but this is too much.

I tried being a good citizen by doing some traffic shaping with wondershaper at my router, giving aMule traffic low priority compared with web surfing. This isn’t trivial, because aMule uses varying ports at both ends of a connection. I used netstat --program to find the peers and then lowered their priorities by host. This helped matters a bit (reducing latency to about half a second), but surfing still felt sluggish.

So I pulled out the big(ger) guns. I downloaded the aMule source code and found this nice little function in Preferences.cpp:

// Here we slightly limit the users' ability to be a bad citizen: for very low upload rates
// we force a low download rate, so as to discourage this type of leeching.
// We're Open Source, and whoever wants it can do his own mod to get around this, but the
// packaged product will try to enforce good behavior.
//
// Kry note: of course, any leecher mod will be banned asap.
void CPreferences::CheckUlDlRatio()
{
	// Backwards compatibility
	if ( s_maxupload == 0xFFFF )
		s_maxupload = UNLIMITED;

	// Backwards compatibility
	if ( s_maxdownload == 0xFFFF )
		s_maxdownload = UNLIMITED;

	if ( s_maxupload == UNLIMITED )
		return;

	// Enforce the limits
	if ( s_maxupload < 4  ) {
		if ( ( s_maxupload * 3 < s_maxdownload ) || ( s_maxdownload == 0 ) )
			s_maxdownload = s_maxupload * 3 ;
	} else if ( s_maxupload < 10  ) {
		if ( ( s_maxupload * 4 < s_maxdownload ) || ( s_maxdownload == 0 ) )
			s_maxdownload = s_maxupload * 4;
	}
}

So, commented out the 'Enforce the limits' section, compiled, and viola -- your basic leech mod. Now I was free to set my upload rate to 1 kB/s without affecting my download rate, but I didn't want to do that because:

I didn't want to be a leech. As a citizen, I try to follow a set of rules that I believe is scalable. What I mean by 'scalable' is this: If everyone followed my set of rules, the community would prosper. My rules aren't necessarily the official ones. For example, I always vote at elections because I believe that if everyone voted it would be good for society. There's no rule that says you have to vote, but it's scalable so I do it. Anyway, being a leech on a P2P network is certainly not scalable.
aMule has a 'credit' system: When you are on someone's queue waiting to download a chunk, uploading him some data gives you credit. You use this credit to jump ahead in his queue. The queue wait is typically very long, so jumping ahead is significant. Hence your effective upload rate still affects your download rate, even with the mod.
Leeching clients are banned when they are discovered. As far as I can tell, this is not done automatically, but rather manually by people who notice skewed upload/download ratios.

Fortunately, there is good middle ground: Setting my upload rate to 7 kB/s. At this rate web surfing is smooth, I still share most of my bandwidth, and chances of getting detected are slim to none. Best of all, this is scalable behavior: If everyone set their upload rate to 3 kB/s under their max rate, and got max down rate in return, things would be terrific.

Posted in Software | 2 Comments »

Fwd: Spread this number

May 1, 2007 – 5:02 pm

09 F9 11 02 9D 74 E3 5B D8 41 56 C5 63 56 88 C0

This is the Processing Key that will allow viewing HD-DVD movies on Linux. What’s happening now with this key is similar to what happened a couple of years ago with the DVD key and DeCSS.

For more, see here and here. The second link is an account of how the key was recovered by the guy who did it. Spread the word!

Posted in Software | No Comments »

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: x^2-x-1 = 0 . 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:

$f\prime(x) = f(x)$

-
The solution of this particular equation is f(x) = c e^x , 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:

$\frac { \partial^2 f } { \partial t^2 } = v^2 \, \frac { \partial^2 f } { \partial x^2 }$

-
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:

$f[j] = f[j] + \frac { \partial f } { \partial t } [j] * dt$

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

$f[j] = \frac{f[j-1] + f[j+1]}{2} + \frac { \partial f } { \partial t } [j] \; dt$

-
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!

Posted in Mathematics, Physics, Uncategorized | 1 Comment »

Test your Math skills

March 17, 2007 – 7:04 pm

Here’s a question that was asked in a recent oral exam for a Master’s degree in Mathematics.

Let $f : [0,1] \rightarrow \mathbb R$ be a real function such that has a limit at each point. Does have at least one continuity point?

Everything you need to solve this question is covered in the first year or so of undergraduate math. Don’t continue reading if you want to try solving it yourself…

Solution. We will show that has a continuity point. has a limit at each point, so let $F : [0,1] \rightarrow \mathbb R$ be a function such that

$\forall x \in [0,1]$ $F(x) = lim_{y \rightarrow x} f(y)$

Lemma.

is continuous.

Proof. Choose some $x_0 \in [0,1]$ . Intuitively, F(x_0)

’s limit, so in a small enough surrounding of x_0

will be close to F(x_0)

. Hence, the limits F(x)

will also be close to F(x_0)

, and thus

is continuous at x_0

.

Formally, let $\epsilon > 0$ . Then there exists $\delta > 0$ such that

$\forall x \in (x_0-\delta,x_0+\delta) \setminus \{x_0\} \; |f(x)-F(x_0)|<\epsilon$

Which means that for all such x

$|lim_{y \rightarrow x} f(y)-F(x_0)|<\epsilon$

$|F(x)-F(x_0)|<\epsilon$
And the Lemma is proved.

Our purpose now is to show that there is a point

such that

. Let’s count the points of

that are ‘far away’ from the limit

. Choose some $\epsilon>0$ and define the set:

$A_\epsilon = \{ \, x \, | \, f(x) \, > \, F(x) + \epsilon \, \}$

Lemma. $A_{\epsilon}$ is finite.

Proof. Suppose $A_{\epsilon}$ is infinite, then because [0,1]

is compact we can find a series $\{x_n\} \subset A_\epsilon$ such that $x_n \longrightarrow x_0$ and $x_n \neq x_0$ for some $x_0 \in [0,1]$ .

Therefore, $f(x_n) \longrightarrow F(x_0)$ . But for large enough , we must have

$f(x_n) \, > \, F(x_n) + \epsilon \, > \, F(x_0) + {\epsilon \over 2}$
where we have used

’s continuity at x_0

. Thus we reach a contradiction, and $A_\epsilon$ must be finite.

Now let’s take $A = \bigcup_{n=1}^{\infty} A_{1/n}$ . Then from the lemma we have that A is an enumerable set. A also contains all the points for whom $f(x) \, > \, F(x)$ .

Likewise we can define

$B_\epsilon = \{ \, x \, | \, f(x) \, < \, F(x)-\epsilon \, \}$
$B = \bigcup_{n=1}^{\infty} B_{1/n}$

And together with A we find that

$A \cup B = \{\, x \, | \, f(x) \, \neq \, F(x) \, \}$
But $A \cup B$ is enumerable, so $[0,1] \setminus (A \cup B) \neq \emptyset$ , so there exists a point $x_0 \in [0,1]$ such that f(x_0) = F(x_0)

is a continuity point for

.

QED

Posted in Mathematics | No Comments »

Recommended Todo web app: Vitalist

February 21, 2007 – 5:44 am

I’ve been looking for a good application to manage todo lists — or Getting Things Done (GTD) as it’s called these days. I have been using Remember The Milk but found it annoying. Yesterday I discovered Vitalist and, well, I was blown away. It’s as though they had found a way to read my mind, searched through all the junk, collected any thought I ever had about The Perfect Todo Application, and then went ahead and materialized the whole thing in the form of a sleek web app, adding some brilliant ideas in the process. Brilliant.

So what’s the big deal about a todo app? For starters, It has a simple and effective workflow concept behind it: New todos usually go into the Inbox. Adding a todo is quick and painless. It can be done using the app, or you can send in the item by email. So for example if you’re somewhere listening to the radio, and you hear a song you’d like to download later, you can use your phone to send in an item that says ‘download that song’ — and it’ll show up in you Inbox. Neat, huh?

The next step is to catalog the items in the Inbox. Here you have several options:

Contexts: Catalog items by the context they should be performed in. For example, create a @Phone context that includes all items that require making a phone call. Then, when you’re at a phone, you can quickly look these up.
Projects: Just a simple grouping of items into projects. Contexts and Projects are orthogonal: Any item can belong to both a project and a context.
Waiting: Items that you delegate to someone else go into the ‘Waiting’ state until that person gets back to you.
Someday: Stuff you want to do someday but not right now. Vitalist will remind you every once in a while about these.
Reference: Little bits of information you want to store somewhere. For example, when I order heating gas for my apartment they give me a confirmation number for the order. I never know where to put it. Now I can just tap the number into my phone and email it to my Inbox. The next time I’m at a computer, I’ll move the item into the Reference list for keeps.

After GMail, this is definitely my ‘best web app around’.

Posted in Personal, Software | 1 Comment »