Tutorial 2: Solving Equations

Introduction To Calculus, Carleton College Math 111 01-02, Fall 2007, Prof. Joshua R. Davis

revised Friday 2007 October 12, 11:45 AM

The purpose of this tutorial is to explore how Mathematica can be used to solve equations, both exactly and numerically. Its numerical solution methods are based on (modern versions of) Newton's method for finding the zeros of a function, which we study in Section 4.8.

Exact, Algebraic Solutions

Try the following cells to see how low-degree polynomials are solved. The first two examples should be familiar to you; they have been known since antiquity. The third and fourth are more complicated; they were only discovered in 1545 by an Italian mathematician named Cardano.

Clear[x, a, b]
Solve[a x + b == 0, x]

Clear[x, a, b, c]
Solve[a x^2 + b x + c == 0, x]

Clear[x, a, b, c, d]
Solve[a x^3 + b x^2 + c x + d == 0, x]

Clear[x, a, b, c, d, h]
Solve[a x^4 + b x^3 + c x^2 + d x + h== 0, x]

The bad news is that there is no algebraic formula for solving general degree-5 polynomials. The impossibility of finding an algebraic formula for degrees greater than 4 was proved by Nils Abel and Evariste Galois in the early 1800s. This surprising theorem from higher mathematics is mentioned on page 210 of your textbook.

Approximate, Numerical Solutions

Non-algebraic equations, involving functions like sin x and e^x, are almost always impossible to solve exactly using algebra. Even when they can be solved exactly, the answer often involves special mathematical constants that are difficult to get ahold of anyway. Consider this:

Clear[x]
Solve[x + E^x == 0, x]

The answer is a weird constant that Mathematica calls -ProductLog[1]. To find an approximate numerical value for it, use the N function:

N[Solve[x + E^x == 0, x]]

An alternative strategy is to graph the function and zoom in on its zero:

Plot[x + E^x, {x, -1, 1}, PlotStyle -> Thick]

Plot[x + E^x, {x, -0.6, -0.5}, PlotStyle -> Thick]

Plot[x + E^x, {x, -0.57, -0.565}, PlotStyle -> Thick]

Newton came up with another method for finding zeros of a function, called Newton's method. The idea is to pick a starting x-value, find the tangent line to the curve at that x, and find the x-intercept of that tangent line. Call this x-intercept the "new x". Then do the same thing with the new x: Find the tangent line there, intersect the line with the x-axis, and call the answer the "new x". If you keep doing this, then eventually the x-values may home in on a zero of the function. (It doesn't work for all functions, or for all starting x-values.)

Newton's Method Example

For an example, let's work with this function:

y = -x^5 + 4x^2 + 4x - 2
Plot[y, {x, -2, 2}]

Since it's a degree-5 polynomial, we don't have a formula to solve it exactly. Let's try Mathematica anyway:

Solve[y == 0, x]

Mathematica returns a weird answer that is, in fact, vacuous. It's not helpful. So let's do Newton's method. First we take the derivative. Did you know that Mathematica's D function computes derivatives for you?

dydx = D[y, x]

Now let's set up a starting x-value called oldX, setting it to -1 since the function seems to have a zero near there (based on the graph above). Then we plug this value into y using the weird /. syntax that we saw in the first tutorial. We also plug it into dy/dx. Then we compute a new x-value called newX, following the Newton iteration formula at the bottom of page 335. We use the N function here to get a decimal number for the new x-value.

oldX = -1
yOfOldX = y /. x -> oldX
dydxOfOldX = dydx /. x -> oldX
newx = N[oldX - yOfOldX / dydxOfOldX]

Now that we have a new x-value, we can copy and paste it into the same cell above for oldX, and get a new new x-value. Do this; you should get about -1.09614. Then copy and paste again. Keep repeating until the answer stabilizes -- meaning that you keep getting the same answer, as far as you can tell.

Mathematica defaults to reporting numerical answers to five decimal digits. If you need more digits shown, then you can add a second argument to the N function telling it how many digits to show:

oldX =-1
yOfOldX = y /. x -> oldX
dydxOfOldX = dydx /. x -> oldX
newx = N[oldX - yOfOldX / dydxOfOldX, 8]

The answer seems to stabilize at -1.0958077.

The function seems to have another zero around 0.4, based on the graph. Let's repeat the process starting there, in an attempt to find this second zero:

oldX =0.36698565332479416`
yOfOldX = y /. x -> oldX
dydxOfOldX = dydx /. x -> oldX
newx = N[oldX - yOfOldX / dydxOfOldX, 8]

The answer seems to stabilize around 0.3669856.

There seems to be another zero near x = 1.8. Find that one too.

Mathematica's NSolve Function

The task of finding numerical solutions to equations is so common that Mathematica automates the process above for you, in its NSolve function.

y = -x^5 + 4x^2 + 4x - 2
NSolve[y == 0, x]

This immediately reports the three real answers we found above, along with two non-real complex answers. The NSolve function is so easy that it appears magical. The important thing to remember is that it is really just doing a lot of numerical approximations, using something very similar to Newton's method.

Here is the other example, from earlier:

y = x + E^x
NSolve[y == 0, x]

Here is another example. Do you agree with the solutions it found? Can you find more solutions yourself? (Yes, you should be able to.)

y = x * Cos[x]
NSolve[y == 0, x]

Here is a more complicated example; Mathematica just seems to give up, which I don't understand.

y = x * Tan[x] + E^x
NSolve[y == 0, x]