**Solving Quadratic Equations**

These are equations of the form

=

where a is not 0.

A typical solution method is to use the quadratic formula, which gives the following solutions

(roots of the equation):

=

=

If the discriminant

is zero, then there is really only one root (of multiplicity 2).

The roots can also be non-real complex numbers, but in this document we only investigate

the real solutions.

Although the quadratic formula is a standard method to use by hand, it can sometimes get us

into trouble if we try to use it on the computer. Consider the following examples:

Example 1:

=

Thus a = 3, b = 5, and c = -2.

The quadratic formula works fine:

=

=

=

=

=

=

=

=

We can also easily solve this quadratic equation by factoring:

=

=

=

or

=

=

=

The Fortran program quadeq1.f90, which uses the standard quadratic formula, produces

the same answers, though of course it has to approximate the 1/3. Remember that

single precision reals only give 7 decimal digits of precision.

Example 2:

=

Thus a = 1, b = -10^6, and c = 1. Note that we have used ^ for exponentiation here,

though Fortran would use ** instead.

The quadratic formula, if used in a Fortran program with single precision reals, gives:

=

=

=

=

=

=

Note that the discriminant is

=

This is because we only have 7 significant digits, resulting in the 4 not affecting the result.

Note that the quadeq1.f90 program, which directly uses the quadratic formula as we just

used it by hand, gives the same answers:

* Program to find real solutions of a quadratic equation*

* Enter (nonzero) coefficient a: *

*1.0*

* Enter coefficient b: *

*-1.0E6*

* Enter coefficient c: *

*1.0*

* Two real solutions: 1000000. and 0.*

How good are our 2 solutions? One way to tell is to check the answers in the original equation.

Let f(x) =

We then calculate

=

=

=

and

=

=

In both cases we get 1, not the desired 0, so our solutions are only approximate.

Question: Can we get better solutions?

We can probably do better in Fortran by using double precision reals, but we could then

find a different example for which even double precision gives answers that are further off

than we would like. We could also use a hand calculator or a software package such as

Mathematica. If we want to stick to single precision reals in a Fortran program, we need a

smarter algorithm!

So. how do we work smarter?

The problem with the quadratic formula approach was found in our second root:

=

=

Since the square root of the discriminant was very close to 10^6, closer than could be captured

with only 7 digits of precision, the difference in the numerator was reported as 0. This problem

did not occur in the calculation for the first root as the 2 terms in the numerator had the same sign.

Therefore, we can use the first root with some confidence that it is reasonably close to the exact

root, but the second root really is not the reported result of 0.

What can we do about that second root?

It turns out that there is a helpful equation for the product of the 2 roots. We can easily derive

this from the formulas for the roots. If we multiply the 2 roots we get:

=

=

=

=

=

Thus

=

Since we already have a pretty good (though not perfect) answer for the first root, we find the

second as

=

In general, find one root using the quadratic formula, being careful to calculate the one with

two terms on the top with the same sign. Then the second root is calculated from the first

using the above simple formula.

For our example problem this gives:

=

=

=

Is this a better answer than 0? Remember that f(0) = 1, not 0 as desired.

Now we try

=

=

=

This gives a result very close to the desired value of zero! Thus we do indeed have a

better answer.

Note that this better answer can also be found by approximately factoring the quadratic:

=

This matches very closely the desired

Only the middle coefficient is different and only by a very tiny percentage.

Example 3:

=

This is really a linear equation, not a quadratic equation, since a = 0. However, we could

try it in the quadratic formula to see what happens. Of course, the result is that we try to

divide by 0, so we cannot get an answer this way. If you use the quadeq1.f90 program, it

produces the following:

*Program to find real solutions of a quadratic equation*

* Enter (nonzero) coefficient a: *

*0.0*

* Enter coefficient b: *

*1.0*

* Enter coefficient c: *

*2.0*

* Two real solutions: NAN and -INF*

=

=

=

=

NAN (not a number)

=

=

=

=

-INF (negative infinity)

Example 4:

=

If we try this monster directly with the quadratic formula, as in quadeq1.f90, we get:

=

=

The problem here is that numbers around 10^60 cannot be in single-precision reals.

On our system these reals can only hold numbers up to about 10^38. Thus, the calculation

under the root sign cannot be carried out in our Fortran program, though it does report

the following as its answers:

*Program to find real solutions of a quadratic equation*

* Enter (nonzero) coefficient a: *

*1.0E30*

* Enter coefficient b: *

*3.0E30*

* Enter coefficient c: *

*2.0E30*

* Two real solutions: INF and -INF*

Of course, we know how to solve the above quadratic equation. We would just divide both

sides by 10^30 and get:

=

=

=

or

=

In a Fortran program, we would typically scale (or normalize) the coefficients a, b, and c by

replacing each by its old value divided by the maximum of the absolute values of a, b, and c.

Let

=

We then use

in place of the original a, b, and c, respectively.

In our example, this would give:

=

The new values of a, b, and c are then:

=

=

=

The quadratic formula then gives:

=

=

=

=

=

=

=

=

=

=

This gives the correct answers.

Conclusions:

In using the quadratic formula to solve quadratic equations on the computer we should use

both the normalization procedure shown above and the method where only the better root

is found using the quadratic formula, with the second one found from the first using:

=

Do the normalization first.