# Numerical Computing: Solving Quadratic Equation

This topic is 3403 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

First, THIS IS HOMEWORK! I AM NOT LOOKING FOR SOLUTIONS, JUST A PUSH IN THE RIGHT DIRECTION. Now, I am taking a class in numerical computing, and we are discussing floating point (IEEE 754), loss of significance, and reducing round-off errors. In a homework problem, we are asked to write a routine to compute the two roots of a quadratic a^2 + bx + c = 0 while reducing round-off error and writing 'efficient' code (relative to matlab's vectorization ability). Now, we all know that the two roots of a quadratic equation can be found using (-b +/- sqrt(b^2 - 4ac)) / 2a. What I am attempting to do is recognize possible areas where precision can be lost, and how to reduce them. I see two spots. First, if b^2-4ac are close, or are both very large, significant round-off errors may occur. Second, if b > 0 and we are looking at the first root (+), we might run into round-off errors. Finally, if b<0 and we are looking at the second root (-), we could run into the same issues. The last two seem fixable, but I can't seem to get rid of the first. I've been trying to find clever ways to multiply by 1 to translate the inner value into an addition ... but then I get hung-up in the cases if a or c are positive/negative. Ultimately, I think I might just be approaching this problem wrong. Does anyone have any experience with numerical computing that could help me out? I've been trying to look and see if using a Taylor expansion anywhere might help, but so far, no dice. Thanks Corey

##### Share on other sites
In another book, I found a specific example of this problem, but no general solution. The recommendations for solving the specific problem were,

1) Use an alternative formula. Note that the product of the two roots must equal c. So r2 = c / r1

2) When looking at -b +/- sqrt(b^2 - 4ac), note that we can re-write b as sqrt(b^2). sqrt(z+e)-sqrt(z) can be rewritten to e/(sqrt(z+e)+sqrt(z)) -- in this case, e should be 4ac.

3) Or, use a Taylor series. If f(x) = sqrt(x), then f(x+e)-f(x) = f'(x)e + 1/2 * f''(x)e^2 + ...

Now, my issue still remains with the (b^2 - 4ac). Are there no numerical issues with this part of the problem?

##### Share on other sites
If you really care that much about precision, using float's is wrong. Use fixed-point instead: http://en.wikipedia.org/wiki/Fixed-point_math

##### Share on other sites
You might have more luck in the Math and Physics forum.

One thing you didn't mention is whether you're working with complex numbers or not. If it's just real numbers, you have the case where there are no roots.

Off the top of my head, you're going to lose precision (that you might be able to avoid by rearrangement) when you divide (by a small/big number?) or when you subtract (and get a small number?). The +/- aside, that means you should be looking at what happens when you divide by 2a and what happens with the b^2-4ac. That's my intuition but nothing much else is occurring to me either. >_<

##### Share on other sites
Quote:
 Original post by BlobbersonIf you really care that much about precision, using float's is wrong. Use fixed-point instead: http://en.wikipedia.org/wiki/Fixed-point_math

You are missing the point of the exercise. The whole idea is that I am trying to use floats and avoid catastrophic collisions. Being homework, I can't just change the problem on a whim. It is a learning exercise.

Anyway, I found a solution online, looking through Numerical Recipes in C. It seems that you break it down by the sign of b. If sign(b) is positive, if you take the positive root. If b is negative, you take the negative root -- so you don't end up with collisions. This solution, however, still ignores the work done inside the square root ... so I don't know if it is a full answer. Anyway, I codified it as follows

b >= 0   t = -0.5 * (b + sqrt(b^2 - 4*a*c));   first_root = t/a;   second_root = c/t;else   t = -0.5 * (b - sqrt(b^2 - 4*a*c));   first_root = t/a;   second_root = c/t;end

##### Share on other sites
When finding roots you do usually have the option of using Newton Raphson to increase the precision of the result.

##### Share on other sites
The Numerical Recipes code mainly tries to avoid the case where 4ac is almost zero, so that you subtract b from almost b, ie b - sqrt(b^2-epsilon). When b^2 is approximately equal to 4ac, then how big of a problem is it compared to the case when 4ac is almost zero?

• 19
• 10
• 19
• 14
• 20