Solving Quadratic Equations

A quadratic equation is an expression of the form $y=ax^2 + bx + c$. Solving one of these equations means finding the value of $x$ that returns $y$ equal to zero.

x = seq(-3, 3, len = 100)
a = 1
b = -2
c = -4
y = a * x^2 + b * x + c
y[1:10]
##  [1] 11.000 10.519 10.045  9.579  9.119  8.668  8.223  7.786  7.356  6.934

plot(x, y)
abline(h = 0)
plot of chunk unnamed-chunk-2

Let's zoom in a bit

x = seq(-2, -1, len = 100)
y = a * x^2 + b * x + c
plot(x, y, type = "l")
abline(h = 0, col = "blue", lwd = 3)
plot of chunk unnamed-chunk-3

Let's make a function to do some of this

quadratic = function(x, a, b, c) {
    return(a * x^2 + b * x + c)
}
quadratic(x, a, b, c)[1:10]
##  [1] 4.000 3.939 3.879 3.819 3.759 3.700 3.640 3.581 3.522 3.463

Now we can do...

x = seq(-10, 10, len = 100)
plot(x, quadratic(x, a, b, c), type = "l")
abline(h = 0, col = "blue", lwd = 2)
plot of chunk unnamed-chunk-5

And we can see this equation has two solutions where it crosses $y=0$. Let's compute them exactly.

The solutions of a quadratic equation is given by the well-known formula: $$ x = \frac{{- b \pm \sqrt {b^2 - 4ac}}}{{2a}} $$ so we will duplicate this in R.

d = sqrt(b^2 - 4 * a * c)
s1 = (-b + d)/(2 * a)
s2 = (-b - d)/(2 * a)
s1
## [1] 3.236
s2
## [1] -1.236
plot(x, quadratic(x, a, b, c), type = "l")
abline(h = 0, col = "blue", lwd = 2)
abline(v = c(s1, s2), col = "green", lwd = 2)
plot of chunk unnamed-chunk-6

Again, we should put this into a function.

quadsolve = function(a, b, c) {
    d = sqrt(b^2 - 4 * a * c)
    s1 = (-b + d)/(2 * a)
    s2 = (-b - d)/(2 * a)
    return(c(s1, s2))
}
quadsolve(a, b, c)
## [1]  3.236 -1.236
quadsolve(-3, 1, 5)
## [1] -1.135  1.468

We can test our function by feeding the answer from quadsolve into quadratic and seeing if we get zero.

quadratic(quadsolve(a, b, c), a, b, c)
## [1] 0 0
quadratic(quadsolve(-3, 1, 5), -3, 1, 5)
## [1] 1.776e-15 8.882e-16

It's not exactly zero. But close enough.

So now we can solve any quadratic equation. Or can we?

quadsolve(1, 2, 3)
## Warning: NaNs produced
## [1] NaN NaN

What has happened here? Let's plot it, with some different limits on the y axis:

x = seq(-5, 5, len = 100)
plot(x, quadratic(x, 1, 2, 3), ylim = c(0, 10))
abline(h = 0, col = "blue", lwd = 2)
plot of chunk unnamed-chunk-10

Now we can see that the curve never crosses the $x$ axis, and so there are no real solutions to the equation. If d is negative in the expression, R returns NaN for "not a number". We can test for this and create an error. We'll also make another improvement by vectorising the plus-or-minus calculations.

quadsolve = function(a, b, c) {
    d2 = b^2 - 4 * a * c
    if (d2 < 0) {
        stop("No real solutions")
    }
    d = sqrt(d2)
    s = (-b + c(-1, 1) * d)/(2 * a)
    return(s)
}
quadsolve(1, 2, 3)  # expect error
## Error: No real solutions

Although there are no real solutions in this case, we can find solutions with complex numbers. In this case we should probably test and return a warning. We can persuade sqrt to not return NaN by converting its argument to a complex number.

quadsolve = function(a, b, c) {
    d2 = b^2 - 4 * a * c
    if (d2 < 0) {
        warning("No real solutions")
        d2 = as.complex(d2)
    }
    d = sqrt(d2)
    s = (-b + c(-1, 1) * d)/(2 * a)
    return(s)
}

This should now return complex numbers for quadratics that don't cross the $x$ axis.

quadsolve(1, 2, 3)
## Warning: No real solutions
## [1] -1-1.414i -1+1.414i

and we can check the solutions produce zero values of the function. In this case we get near-zero complex values.

quadratic(quadsolve(1, 2, 3), 1, 2, 3)
## Warning: No real solutions
## [1] -4.441e-16+0i -4.441e-16+0i

Further Thoughts

  • Now write a function to evaluate the cubic function, $y=ax^3+bx^2+cx+d$ and plot a few curves.
  • How would you write a function to compute a polynomial of any order?
  • Read the help for the curve function and plot some of the graphs in the previous examples