algorithm quadratic equation MATLAB

Question

What condition should I put in code of matlab so that get the exactly solutions of a quadratic with these formulas:

x1=(-2*c)/(b+sqrt(b^2-4*a*c))
x2=(-2*c)/(b-sqrt(b^2-4*a*c))

Directly implementing these formulas I don't get the correct solution in certain cases such x^2-1000001x+1

Thank you very much for your help

Answer 1

The correct set of formulas is

w = b+sign(b)*sqrt(b^2-4*a*c)

x1 = -w/(2*a)

x2 = -(2*c)/w

where sign(b)=1 if b>=0 and sign(b)=-1 if b<0.

Your formulas as well as the standard formulas lead to catastrophic cancellation in one root of b is large wrt. a and c.

---

If you want to go to the extremes, you can also guard against over- and underflow in the computation of the term under the square root.

Let m denote the maximum size of |a|, |b| and |c|, for instance the maximum of the exponent in their floating point representation, or of their absolute value... Then

w = b+sign(b)*m*sqrt( (b/m)*(b/m)-4*(a/m)*(c/m) )

has a term between -10 and 10 below the root. And if this term is zero, then that is not caused by underflow.

Answer 2

You're dealing with floating point arithmetic in matlab, so exact solutions are not guaranteed. (I.e. it's possible that every single floating point value will cause round-off error that gives non-zero answer when you plug into your original quadratic equation). A better way to check whether you have found a solution to an equation in floating point is to use a tolerance, and check whether the absolute value of your answer is less than the tolerance.