1
$\begingroup$

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 ?

$\endgroup$
1
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ When I substitue it back I find that one of the roots should be calculated by using $root_1 = (-b + sqrt(delta)) / 2 * a$ and the other one $root_2 = 2 * c / (-b - sqrt(delta))$. Only these two roots give me numbers that are cloe to 0. $\endgroup$ – Layla Jul 12 '19 at 14:02
  • $\begingroup$ Is this makes sense ?. Thats is also what it says in wiki I guess $\endgroup$ – Layla Jul 12 '19 at 14:02
  • $\begingroup$ @Reign Yes, only one of the roots is close to zero. $\endgroup$ – David Richerby Jul 12 '19 at 14:45
1
$\begingroup$

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

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Hmm yes you are right..I noticed that I missed a parentheses in my calculations which led this kind of a result. Thanks for noticing and sharing it to me $\endgroup$ – Layla Jul 14 '19 at 13:07
0
$\begingroup$
  1. In order to check it, compute roots yourself using pencil and paper
  2. Floating-point arithmetic: Accuracy problems
| cite | improve this answer | |
$\endgroup$
0
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ I tried your case and that is really really interesting. Theres a difference but I dont know can we consider its as huge. What numbers should I put for instance ? $\endgroup$ – Layla Jul 12 '19 at 18:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.