8

I have a set of quadratic coefficients a, b, and c that need to be solved in order to determine the intersection between a parabola and a line on a plane. These coefficients have "arrived" from a different algorithm.

For context, the actual quadratic being solved is: $$a = 4.3162933554 \times 10^{-11}$$ $$b = 1.0361118521 \times 10^{6}$$ $$c = -4.5813932360 \times 10^{5}$$ where $f(x)=ax^2+bx+c$.

My problem, however, is that when my parabola is so "stretched out" that it starts behaving like a line around the area of interest, the "a" coefficient becomes very small, and actually disregarding it completely and instead solving linearly for b and c works perfectly for the required intersection.

In this case what is actually happening with the formula is that the $-b+\sqrt(b^2-4ac)$ term is being calculated as 0.0 by the computer! (I've confirmed that this isn't specific to my implementation as my calculator does the same thing! Also, all of this is involving double floating point precision)

My question is, how can I approach finding this root in an optimal way? Also, are there any other cases where the formula may fail?

I've thought about disregarding the $a$ term and simply solving the equation linearly, but then at what point do you decide to do that? By extension, when should you solve a cubic equation quadratically?

Gary Allen
  • 438

1 Answers1

13

You can switch over to the modified solution formula, according to binomial formulas, $$ x=\frac{-b+\sqrt{b^2-4ac}}{2a}=-\frac{2c}{b+\sqrt{b^2-4ac}}. $$ The cancellation case in the first formula will give a numerically stable computation in the second formula, giving a value close to $-\frac{c}{b}$, the solution of the linear approximation.

In general, most of the cases with numerical issues when computing the roots, in the case of real roots, are avoided by first determining the non-cancelling variant of the numerator $$ u=-b-{\rm sign}^+(b)\sqrt{b^2-4ac},~~~~{\rm sign}^+(0)=+1, $$ and then computing the two roots as $$ x_1=\frac{u}{2a},~~~ x_2=\frac{2c}{u}. $$ After that, the stabilization discussion turns to the square root, to be able to keep computing it in the case of floating-point over- or underflow.

Lutz Lehmann
  • 131,652
  • At what point do you switch over to this formula? Is there a value I can check for to avoid needless computations? – Gary Allen Jan 26 '21 at 07:37
  • 1
    The usual way is that that you take the square root with the same sign as $-b$ for the first root. That is, you compute $u=-b-{\rm sign}(b)\sqrt{b^2-4ac}$ to then get the roots as $\frac{u}{2a}$ and $\frac{2c}{u}$. – Lutz Lehmann Jan 26 '21 at 07:40
  • Thanks a lot! This post also helped clear it up for me: https://people.csail.mit.edu/bkph/articles/Quadratics.pdf – Gary Allen Jan 26 '21 at 07:42
  • 1
    This is great! One thing that just bit me was that you don't actually want the $sign$ function, as it gives $0$ for $b=0$, which is incorrect. Instead, you want $-1$ in the zero case. – Daskie Jan 23 '23 at 07:03