Numerically Stable Solution For The Quadratic Equation

AllPages RecentChanges Links to this page Edit this page Search Entry portal Advice For New Users

http://ow.ly/i/24XQl

I said

• Of course, that formula can't be used
  in industrial strength programming
  because of instabilities.

• Yeah...I'm going to have to take your word
  on that one. I'm looking as confused as my
  students right now ...

The classic formula is this:

This is to find the value(s) of x (if any) such that $f(x)=ax^2+bx+c$ where $a{\ne}0.$

So far so good.

But when we're using this in the real world, we usually don't know the absolute exact values for a, b and c. Usually these are the results of measurements, or possibly even the results of a model we have, which itself might not be exact. So we really ought to allow for the values of a, b and c not to be exact.

Why is that a problem?

Well, usually it isn't, but we can explore the idea that it might be. We can do that by jiggling the values of a, b and c by just a little, and seeing how much the answers change. That makes for a fun exercise. If you're allowed to change a, b and c by up to 0.001, how much can you make the values of x change?

Here's how we do that.

Suppose firstly that a is really, really small. Then when we divide by 2a we're dividing by a small number, which is the same as multiplying by a big number. That will magnify any instabilities in the top line, so if we're after instabilities, that's a good thing.

At this point you might think, "What? Why are instabilities a good thing? Well, instabilities are definitely a bad thing, but what we're trying to do here is show that the usual formula for the quadratic equation is not suitable as it stands for industrial strength applications. The reason it's not suitable is because it is prone to numerical instabilities. To demonstrate this point we need to show how and why it's prone to numerical instabilities, and so we want to find numerical instabilities. So finding them is, for now, good. Clear? Phew ...

The next question is how to introduce instabilities to the top line.

The discriminant is $b^2-4ac$ and we're going to take the square root of that. If a is small, and c is small, then $4ac$ is certainly small. If b is reasonably large, then we're subtracting a very small number from a very large number, and that has problems in computer arithmetic. Random real numbers are represented only approximately in a computer anyway, so subtracting a very small number from a very big number might not ever get noticed.

So the discriminant is very nearly $b^2,$ and when we subtract a very small number, it's still very nearly $b^2.$ Then we take the square root, and get something that's very nearly equal to b.

So on the top line we have the sum of two things that are very nearly equal, and that's OK. But we also have the difference of two things that are very nearly equal, and that's not OK. Remember, computer arithmetic is not exact, and besides, you may have got the wrong value in the 12th decimal place. So when you take the difference of two very nearly equal quantities you get small rubbish.

Now divide by a very small number, and you get big rubbish.

So if a is very small, and 4ac is even smaller, one of your solutions is fine, but the other depends critically on the exact values of a, b, and c, and may turn out to be complete rubbish.

Now awaiting feedback ...

Shortly I'll add to this page how we get around this problem, and some physical descriptions of what's going on. In short, we have two nearly parallel lines and we're looking for an intersection. Move one of the lines by a tiny amount and the intersection can move a very long way.

For reference, here is how to find the stable solutions:

Also, these from Michael Borcherds:

• Would be nice with an example eg 94906265.625x^2 − 189812534x + 94906268.375
• http://www.cs.berkeley.edu/~wkahan/Qdrtcs.pdf
• Also better to find q=-1/2[b+sgn(b)sqrt(b^2-4ac)] and then find x1=q/a and x2=c/q I think
• http://www.nr.com/