Quadratic equations using 2 different approaches

Question

I am reading Mark Newman's Computational Physics and at chapter 4 page 133 in Exercise 4.2 he asks

a) Write a program that takes as input three numbers, a, b, and c, and prints out the two solutions to the quadratic equation $ax^2 + bx + c = 0$ using the standard formula $x = −b± (b^2 − 4ac)^{1/2}/2a$ . Use your program to compute the solutions of $0.001x^2 + 1000x + 0.001 = 0$.

b) There is another way to write the solutions to a quadratic equation. Multiplying top and bottom of the solution above by $-b∓ (b^2 − 4ac)^{1/2} $, show that the solutions can also be written as $x = 2c/−b∓(b^2 − 4ac)^{1/2}$. Add further lines to your program to print these values in addition to the earlier ones and again use the program to

I tried both ways and a) gives me

[-9.99989425e-13 -1.00000000e+00] and

b) [-1.00000000e-06 -1.00001058e+06]

how can I understand which one is correct ? Or why is this happening ?

Answer 1

You can understand which, if either, solution is correct by substituting it back into the quadratic and seeing which one works.

As to why you get different answers, it's to do with the accuracy of floating point representations. Roughly speaking, numbers are represented to some number of significant figures and this can cause problems when you add numbers with different orders of magnitude or subtract numbers that are approximately equal. For example, evaluated left-to-right to one significant figure, $0.3+0.3+1 = 0.6+1 = 2$, whereas $1+0.3+0.3=1+0.3=1$.

This is a particular issue in the calculation of $b^2-4ac$ in your example, where $b=1000$ and $a=c=0.001$: the correct answer is $999\,999.999\,996$, which is very susceptible to even slight rounding errors. You then have to calculate the square-root of that, which gives a number very close to $1\,000$, and compute $1\,000$ minus that number. All of these steps cause problems.

The concepts you need to read up on are loss of significance and numerical stability.

Answer 2

There seems to be something wrong in the first result you obtained - your case (a) - irrespective of questions to do with floating-point precision and cancellations.

$$x=\frac{-b-\sqrt{b^2-4ac}}{2a}$$ should not cause any problems at all, since there is no cancellation. This value is very nearly $\frac{-2b}{2a}$. The fact that you did not get 1E+06 or something very near to it means that your programming must have been wrong in this case.

$$x=\frac{-b+\sqrt{b^2-4ac}}{2a}$$ is indeed expected to cause problems because you are calculating $-b+{\rm{almost\ }} b$, which largely cancels out. So you would expect a peculiar answer in this case. I wouldn't be surprised, though, to find that once you corrected whatever the problem was with your first calculation, the second one got a little bit less strange as a consequence.

Your second result does not have this programming problem. You have got it right.

---

As for your question [in the comments] about what to do in general to get accurate roots, I would suggest taking the non-cancelling variant - thus, $-b-\sqrt{b^2-4ac}$ - to get one root. You can then use the fact that the product of the roots equals $-\frac{c}{a}$ to get the other root.

You would, however, run into different difficulties if $b^2$ were close enough to $4ac$ to give a lot of cancellation under the square root sign. In that case your best bet would be to substitute $y=x+\frac{b}{2a}$, and solve for $y$, since the roots of the quadratic would be very close to $y=0$.

Answer 3

1. In order to check it, compute roots yourself using pencil and paper
2. Floating-point arithmetic: Accuracy problems

Answer 4

This exercise shows you an extremely important concept: Floating point arithmetic always has rounding errors, but the magnitude of those errors can hugely depend on the exact way you arrange your calculation.

First check that both formulas are mathematically the same (common trick: (a-b) = (a^2-b^2) / (a+b)).

Then note that if you calculate both a and b with a reasonably small error, and they are close together, then a-b has a HUGE relative error. Which can totally destroy your calculations.

And think about the shape of the curve 0.001x^2 + 1000x + 0.001. Its close to 0 when x = 0. For positive x it goes to infinity. But for negative x, it goes negative very very quickly until the term 0.001x^2 becomes as large as 1000x with opposite sign, around x =-1000000.