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)


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)


Let's make a function to do some of this

quadratic = function(x, a, b, c) {
return(a * x^2 + b * x + c)
}

##  [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)


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)


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))
}

## [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)


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