Tutorial 3: Solving Equations

Introduction To Calculus, Carleton College Math 111, Prof. Joshua R. Davis

revised 2008 February 14

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:

SetOptions[Plot, PlotStyle->Thickness[Large]];
Plot[x + E^x, {x, -1, 1}]

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

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

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 an unhelpful answer. 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 currentX, setting it to -1 since the function seems to have a zero near there (based on the graph above).

currentX = -1

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, following the Newton iteration formula at the bottom of page 335. We update currentX to this new value. Then we use the N command to show currentX to 8 significant figures.

currentY = y /. x -> currentX
currentDydx = dydx /. x -> currentX
currentX = currentX - currentY / currentDydx
N[currentX, 8]

You should get -1.11111111. Now that we've changed the value of currentX, we can run the same cell again, without alteration, and get yet another value for currentX. Do this; you should get about -1.0961358. Then run the same cell again. Keep repeating until the answer stabilizes -- meaning that it the answers agree, to however many decimal places you want.

In this example Mathematica is showing you what appear to be very complicated numbers, because it is keeping track of the exact value of currentX as a ratio of two very large integers. You don't really want to see that stuff, right? Then add some semicolons as follows. Commands that end in semicolons are performed like any other, but Mathematica does not show you the output.

currentY = y /. x -> currentX;
currentDydx = dydx /. x -> currentX;
currentX = currentX - currentY / currentDydx;
N[currentX, 8]

The answer seems to stabilize at -1.0958077.

Based on the graph, the function seems to have another zero around x = 0.4. Go back to where we set currentX equal to -1, change that to 0.4, and go through the iterations again until you stabilize at an answer. You should end up with about  0.3669856.

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

Mathematica's NSolve and FindRoot Functions

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. Calculus is used "behind the scenes" like this throughout modern computing.

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.

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

However, there's another numerical solving function called FindRoot that we can try. FindRoot requires a seed value to get it started. Let's try -0.1.

FindRoot[y, {x, -0.1}]

A numerical root-finding method that starts from a seed value? It sounds like Newton's method, right? The Mathematica implementation notes at

        http://reference.wolfram.com/mathematica/note/SomeNotesOnInternalImplementation.html#2118

indicate that FindRoot does indeed use a version of Newton's method, at least sometimes.