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.
=
=
=
=
=
=
=
=
We can also easily solve this quadratic equation by factoring:
=
=
=
or
=
=
=
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
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
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:
=
=
=
=
=
=
=
=
=
=
=
=
=
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.