UC SCI 111, 18 sep 2002
Notebook 3: Solving Equations

Needs["Graphics`ImplicitPlot`"]

Mathematica Prerequisites.

Solving equations with Solve, NSolve, and FindRoot.

Equations can be easy to solve, as in (x^2) = 9, or they can be quite difficult, as in (x^4) - 3  (x^2) + 6  x = 17. It all depends on the equation in question. Mathematica can solve many (but certainly not all) equations with its Solve and NSolve commands. Note the double equal sign in the syntax below. It is Mathematica's syntax for writing a mathematical equality (we cannot use = since that is already used for assigning names to quantities).

Solve[(x^2) == 9] // TableForm

Solve[(x^4) - 3  (x^2) + 6  x == 17] // TableForm

The last example shows why one might seek approximate numerical answers. The NSolve command does just that.

NSolve[(x^4) - 3  (x^2) + 6  x == 17] // TableForm

Note that Mathematica will report both real and complex solutions. Note also that NSolve will sometimes return a solution for an equation which Solve fails. If a particularly nasty equation foils both Solve and NSolve, one can seek numerical solutions using the Newton-Raphson method discussed in this Notebook. You will be pleased to know that Mathematica has a built-in command for implementing an algorithm based in part on the Newton-Raphson method, which is handy in such situations. It is your tool of last resort in solving an equation, and it is called FindRoot. It needs an equation and an initial guess to get started. It then iterates until the result no longer changes, and reports the final value. A consequence of such an approach is that only one solution can be found for a given intitial value, even if the equation has many solutions. As an example, here's a nasty equation: sin (x^2) = cos x. Note that both Solve and NSolve are utterly foiled:

Solve[Sin[(x^2)] == Cos[x]]

NSolve[Sin[(x^2)] == Cos[x]]

Now we attack the problem with FindRoot. First make a plot. This will give us a reasonable guess. Here we seek the zeros of the function sin (x^2) - cos x.

Plot[Sin[(x^2)] - Cos[x], {x, -2, 2}]

We use 1 as an initial value, hoping to find the solution near x = 1.

FindRoot[Sin[(x^2)] == Cos[x], {x, 1}]

Beautiful! And a bit easier than working Newton-Raphson by hand.

One can solve a system of equations with these commands as well:

sol = Solve[{2  x + 3  y == 6, x - y == 4}]

We can easily check the result if we recall how replacement rules work:

2  x + 3  y  /.  sol

x - y /. sol

The syntax for solving a system with FindRoot is as follows:

FindRoot[{y == 1/x, Sin[y] == Cos[x]}, {x, 1}, {y, 3}]

Recursive functions using  NestList.

Sometimes, like in the Newton-Raphson method, one likes to perform an operation repeatedly to the outcome of the previous operation. A simple example is to start with the number 1, and then multiply by 2, then multiply the answer again by 2, and keep on going like that. The numbers we get are (Subscript[x, 0]) = 1,  (Subscript[x, 1]) = 2,  (Subscript[x, 2]) = 4,  
      (Subscript[x, 3]) = 8, .... and so on. The relation of each number of its predecessor is of course

(Subscript[x, n + 1]) = 2  (Subscript[x, n])

There is a nice way to view this sequence of (Subscript[x, n])at once. First we introduce the doubling function

f[x_] := 2*x

To get the 10-th iteration of our doubling procedure we simply give the command

Nest[f, 1, 10]

But if you want to see all iterations up to the 10-th, you simply give the command

NestList[f, 1, 10]

Of course doubling a number every time doesn't give very unexpected answers. But if we use more sophisticated functions f, things may different. For example consider the recursion

(Subscript[x, n + 1]) = 2.5*(Subscript[x, n])(1 - (Subscript[x, n]))

with starting value (Subscript[x, 0]) = 0.5.

f[x_] := 2.5*x*(1 - x)

Exercise 1.
a) Print the first hundred iterations of this recurrence. What do you conclude about the long term behaviour of the sequence you found?
b) Change the 2.5 into 3.2 and repeat the experiment.

Plotting curves in the plane with ImplicitPlot.

We have used the Plot command to plot the graphs of functions of one variable. But how can one plot curves that are not the graphs of functions, for example a closed curve such as an ellipse? After loading the Mathematica package Graphics`ImplicitPlot` (which should have occured automatically at the beginning of this session), the command ImplicitPlot is available for this purpose.  For example, here is a plot of the ellipse with equation 4  (x^2) + (y^2) = 4.

ImplicitPlot[4  (x^2) + (y^2) == 4, {x, -2, 2}]

This command takes as its first argument an equation in two variables. Note that equations use the double equal sign (a special syntax is needed since the single equal sign is used for making asignments, a fundamentally different operation).

A handy option for both Plot and ImplicitPlot is PlotStyle, which controls the appearance of the curve. The setting PlotStyle->Hue[1] will produce a red plot, while PlotStyle->Hue[.7] will produce a blue plot.

Exercise 2.
Make an implicit plot of the following curves. Choose the x-range (and y-range if necessary) in such a way that you get a good idea of what the curve looks like.
a) (x^2) - (y^2) = 1
b) 2  (x^2) - 5 x  y  - (y^2) - 3  x + y = 0
c) (x^2) - 4 x y  +  4 (y^2)  - x + y + 5  =  
        0
d) (x^3) + (y^3) - 3 x y = 0 (folium of Descartes)
e) 4 (x^6)+4 (y^6)+12 
           (x^2) (y^4) + 12 (x^4) (y^2) -4  (x^4)-4 (y^4)-8  (x^2) (y^2) +(y^2) = 0
(double cardioid)

Polynomial equations in one variable

Mathematica certainly knows the so-called abc-formula,

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

It also knows the Cardano formulas for cubic equations,

Solve[x^3 + p x^2 + q x + r == 0, x]

And it knows the Ferrari formulae for fourth degree equations as well,

Solve[x^4 + p x^3 + q x^2 + r x + s == 0, x]

Although this is quite impressive, one may ask what the use of this is. Let us continue one more step and try to solve a fifth degree equation

Solve[(x^5) - 2 (x^4) - 6 (x^3) + 13 (x^2) - 4 x - 1 == 0, 
      x]

In this case the programmers did not provide anything but a tautological answer that the solution x is a root of (x^5) - 2 (x^4) - 6 (x^3) + 13 (x^2) - 4 x - 1 = 
        0. Actually, the formulas to solve quadratic, cubic and quartic equations had the only merit that they reduce the solution of these equations to the problem of square root and cube root extraction. The reason that Mathematica does not give us an answer for fifth or higher degree equations is that they cannot be solved generally in terms of repeated root extractions. 

So one has to reconsider what is meant by a solution of an equation. This is done in a branch of mathematics called algebra. However, in many applications people are only interested in numerical approximations of solutions. And so will we. We certainly have numerical solutions in our example,

sol = 
      NSolve[(x^5) - 2 (x^4) - 6 (x^3) + 13 (x^2) - 4 x - 1 == 0, 
        x]

Here is a plot

Plot[
      (x^5) - 2 (x^4) - 6 (x^3) + 13 (x^2) - 4 x - 1, {x, -2.7, 
        3}]

It should be said that we were lucky to get 5 solutions. There might have been less, as is shown by the following example.

sol = NSolve[(x^5) + 3 x - 1 == 0, x]

Here is a plot

Plot[(x^5) + 3 x - 1, {x, -2, 2}]

Again we found 5 solutions, but 4 of them (containing the symbol i) are so-called complex numbers. We shall deal with complex numbers next time. For the moment we shall be interested only in real solutions and note that in the present example there is only one real solution.

Exercise 3.
a) Can you think of a 4th degree equation which has no real zeros?
b) And what about a 5th degree equation?

Interpolation

Interpolation is a procedure which is often applied to give approximate answers to function evaluations or solutions of equations. Of both cases we shall give an example. Interpolation is based on the idea that for very small x-intervals the graph of a function almost looks like a straight line. Here is an example

Plot[Sin[x], {x, 1, 1.1}]

In ancient days, when people had to compile handcomputed tables for the sine,cosine and log function, interpolation was often used.  
Suppose for example that, after diligent computation (using Taylor series for example) someone has found that log(1.8) = 0.587787 and log(1.9) = 0.641854. We now like to evaluate log(1.87). Based on the assumption that the log-function for 1.8≤x≤1.9 looks practically like a straight line, the value of log(1.87) is approximatedly a weighted average of the values log(1.8) and log(1.9).

Exercise 4.
Use interpolation to find an approximate value of log(1.87). Compare this with the actual value as computed by Mathematica.

Interpolation can also be used to find approxiamtions of zeros of functions quickly. Consider the function

f[x_] := (x^5) - 3 (x^3) - 5 x - 1

It certainly has zeros as can be seen from the following plot.

Plot[f[x], {x, -3, 3}]

Let us concentrate on the zero in the neighbourhood of 2.

Exercise 5.
By making plots in smaller and smaller intervals near x=2, determine the value of the zero of f near 2 up to three decimal places.

Let us now proceed a bit more clever. After doing the above exercise you could have reached the point where you plotted f in a very small interval, like [2.068, 2.07]

Plot[f[x], {x, 2.068, 2.07}]

Notice that the graph in this small region almost looks like a straight line. We can now pretend it is a straight line and determine the point of intersection with the x-axis.

Exercise 6.
a) Determine the equation of a straight line passing though the points 
(2.068, f[2.068]) and (2.07, f[2.07]).
b) Determine the x-value of the point of intersection of this straight line with the x-axis and call this point (Subscript[x, 0]) . 
c) Determine f[(Subscript[x, 0])] and compare this with the values f[2.068] and f[2.07].

Conclude from this exercise that with only the knowledge of the values of f at 2.068 and 2.07 we have constructed a point (Subscript[x, 0]), which is a much better approximation of the zero of f near 2. We could now repeat the procedure by choosing another point (Subscript[x, 1]),very close to (Subscript[x, 0]), and perform interpolation between (Subscript[x, 0])and (Subscript[x, 1]).

Exercise 7.
What point (Subscript[x, 1])would you choose?

Newton-Raphson Method.

Consider the function

f[x_] := -1 - 2  x + (x^2) + (x^3)

A plot shows that f has three real roots:

Plot[f[x], {x, -2, 2}]

We can implement the Newton-Raphson method by introducing the recursive function r(x) . Note that we must use the immediate assignment operator = to ensure that the derivative is evaluated immediately.

r[x_] = x - f[x]/(f')[x]

Here we carry out 9 iterations of Newton-Raphson, with a starting value of x = 1. We carry out calcuations to 40 significant digits.

NestList[r, N[1, 40], 9] // TableForm

Exercise 8.
Carry out 9 iterations of the Newton-Raphson method above for the starting values 0, -2, 2. Can you see from the plot of f(x) why the different starting values lead to the various end results?
Suppose now we start with 1/2. In the first few iterations you see some negative numbers, far away from the zeros. Can you explain this by consulting the plot?

Exercise 9.
Use Newton-Raphson to approximate (Sqrt[1.25]) to 40 significant digits. Compare the speed of convergence with the Taylor series you used in the previous Notebook.

Equations in several variables

Consider the set of equations

(x^2)  + (y^2) = 1 and y = 2 (x + 1)

and solve it.

Solve[{(x^2) + (y^2) == 1, y == 2 (x + 1)}, {x, y}]

Exercise 10.
a) Did you expect to get two solutions? In particular, check the correctness of the first solution by substituting it into the equations. 
b) Solve the above system by hand, by elimination of y from the equations.
(This is only to show you that no mystery is involved in the solution of these equations. From now on we do our computation in Mathematica again)

Hopefully you can see that the first solution is an obvious solution. Here is a geometrical picture of the equations

ImplicitPlot[{(x^2) + (y^2) == 1, y == 2 (x + 1)}, {x, -1.5, 
        2}, {y, -2, 2}]

Exercise 11.
We repeat the above experiment, but now with
(x^2) + (y^2) = 1 and y = (x + 1)/4.
a) Can you see an obvious solution?
b) Solve the equations and show that the second solution gives rise to the remarkable equality (8^2) + (15^2) = (17^2).

In the previous exercises we intersected the circle with lines having slopes 2 and 1/4. Let us now take an arbitrary slope, which we call t.

Solve[{(x^2) + (y^2) == 1, y == t (x + 1)}, {x, y}]

The second solution is a point which changes as t changes. It runs along the unit circle. 
We have obtained a parametrisation of the circle other than the usual
x = cos(t),  y = sin(t). 
Just for fun, here is an animation of the running point. Enter the command below (don't try to understand it ). You will see a series of pictures. Double click on one of them.

For[n = -7, 
      n < 7, n++, 
t = 
        Tan[Pi*n/15] // N; 
p1 = 
        ImplicitPlot[{(x^2) + (y^2) == 1, y == t (x + 1)}, {x, -1.5, 
            2}, {y, -2, 2}, 
          DisplayFunction -> Identity]; 
p2 = 
        Graphics[{PointSize[ .035], Hue[ .7], 
            Point[{(1 - t^2)/(1 + t^2), 
                2*t/(1 + t^2)}]}]; 
Show[p2, p1, 
        AspectRatio -> Automatic, Axes -> True, 
        PlotRange -> {{-1.5, 2}, {-2, 
              2}}]
]