Departments of Chemical Engineering and Mathematics
Mathematica Session 4

Numerical computing in Mathematica

1. Solving Equations Numerically

Last session, you tried using the Solve and NSolve commands to solve the equation

(e^x) + (e^(-x)) = 2  x + 4.

This was, of course, entirely unsuccessful. And yet this equation has two perfectly good solutions, as we can see by the fact that the graphs of y = (e^x) + (e^(-x)) and y = 2  x + 4 intersect in two points:

Plot[{(E^x) + (E^(-x)), 2  x + 4}, {x, -2, 3}]

The problem is that Solve and NSolve work only for polynomial equations and closely related types; when there's an exponential, trigonometric or logarithmic function around they often can't cope. In part that's because equations involving those kinds of functions often simply can't be solved using symbolic methods.

In this section, we introduce the idea (already perhaps familiar to some of you) of  "root finding" for equations of this type by numerical iteration. We'll be mostly working with equations of the form f(x) = 0, which means the above equation should really be written

(e^x) + (e^(-x)) - 2  x - 4 = 0.

Simple Iteration

The basic idea behind most numerical equation solving techniques is this: start with an approximate solution, and use it to find a better approximate solution. We make use of the fact that certain functions (not all of them, as you'll see) have the property of converging when iterated.

Iterating a function means evaluating it repeatedly, "feeding back" each piece of output into the function.

Converging is when the results of this process tend to a limit. As a dead simple example of convergent iteration, try

Cos[1.0]

then evaluate

Cos[%]

repeatedly and see what happens. The value you're moving towards is a root of the equation x = cos x.

This method can be used, then, to find solutions of equations of the form x = g(x). So how could we use it to solve something of the form f(x) = 0? A simple idea would be to set g(x) = x + f(x); then the equations x = g(x) and f(x) = 0 would be equivalent. Let's try this for our equation.

Type

f[x_] := (E^x) + (E^(-x)) - 2  x - 4

g[x_] := x + f[x]

Now choose an initial "guess"; -0.7 might be a good one, since from our graph we know this lies close to one of the solutions we're trying to find. Type

g[-0.7]

then evaluate

g[%]

repeatedly.  What happens?

What happens if you use

g[x_] := x + (1/2) f[x]

instead? This, of course, has the same root; does it behave in the same way?

Exercises and explorations 1a

1.	The above can be streamlined quite a bit by using the Mathematica command NestList, which is exactly right for this sort of thing. Try

g[x_] := x + f[x]

NestList[g, -0.7, 10]

and

g[x_] := x + (1/2) f[x]

NestList[g, -0.7, 10]

Note that you can get a quick visualisation of the iterations by feeding the NestList output straight into ListPlot:

ListPlot[%, PlotJoined -> True]

Experiment with different start points for these iterations. Can you get either to converge to the other solution of f(x) = 0, which is, from our graph, at around x = 2.0?

2.	There's a standard diagram for representing iterations of this type, sometimes called a "cobweb" or a "staircase", of which the following is an example:

The straight line is y = x, the curved one is y = g(x) and the solutions correspond to the points of intersection. If we start at a point on y = x and follow the directions "move vertically to the curve, then move horizontally to the line" over and over again, we get a path that corresponds to our iteration.

The following is some "quick and dirty" code for generating an animation that shows a point moving along this path.

\(iterationFigure[g_, x1_, n_, {x_, xmin_, xmax_}] := 
	
      \((xn = Nest[g, x1, n - 1]; 
		
        Plot[{x, g[x]}, {x, xmin, xmax}, 
          PlotStyle -> {RGBColor[0, 0.5, 0], RGBColor[0, 0, 1]}, 
			
          Epilog -> {RGBColor[1, 0, 0], PointSize[0.02], Point[{xn, xn}]}, 

          			PlotLabel -> xn]; 
	\ \ 
        Plot[{x, g[x]}, {x, xmin, xmax}, 
          PlotStyle -> {RGBColor[0, 0.5, 0], RGBColor[0, 0, 1]}, 
          Epilog -> {RGBColor[1, 0, 0], PointSize[0.02], Point[{xn, g[xn]}]}, 
          
			PlotLabel -> xn])\)\)

Type it in, then try

g[x_] := x + (1/2) f[x]

Table[iterationFigure[g, 0.3, n, {x, -2, 2.2}], {n, 1, 10}]

then view the sequence of graphs as an animation.

(If you're dead keen, we have some more sophisticated code that fills in the path as it goes along. It is on the course web page -- ask one of the tutors if you would like to get hold of it.)

Use iterationFigure to investigate the iteration based on g(x) = x + α f(x) for different values of α.  You could try, for example,

g[x_] := x + (1/3) f[x]

Table[iterationFigure[g, 0.3, n, {x, -2, 2.2}], {n, 1, 10}]

Try to find α's that give you really quick convergence (a) to the negative root and (b) to the positive root. [Hint: don't forget that α can be a fraction, and don't forget that it can be negative.]

What's distinctive about the graph of g(x) near the root in the cases that converge quickly? That don't converge at all?

Newton's Method

Iterations of the above type are all very well, but half the time they don't work, and even when they do work they're normally pretty slow. Newton's Method (sometimes known as the Newton-Raphson Method) goes one better; it works most of the time, and it's nearly always very very fast.

The basic idea is the same: given an equation of the form f(x) = 0, do an iteration using a function g(x) = x + α f(x). The only difference is that α is no longer a number, but a function of x. Theoretical considerations show that the value -1/f ' (x) is best, giving g(x) = x - f(x)/f ' (x).

Exercises and explorations 1b

1.	Type

g[x_] := x - f[x]/(f')[x]

and run

NestList[g, 0.3, 10]

Is convergence indeed fast? Can you choose a start point that gives the other root?

Here, Newton's method very quickly converges to more than the 5 decimal place accuracy that Mathematica normally displays. To see more digits, pass the NestList output to the N command, for example:

N[%, 15]

2.	What do you see if you animate

Table[iterationFigure[g, 0.3, n, {x, -2, 2.2}], {n, 1, 10}]

Try focusing in on a smaller interval around the root:

Table[iterationFigure[g, 0.3, n, {x, -0.8, -0.7}], {n, 1, 10}]

Repeat, this time using a different start point and focusing on the other root. Why do you think Newton's method tends to converge so quickly?

3.	Try out Newton's method on some other equations, such as x = cos x, (x^2) = 3, or whatever (remember to put each in the form f(x) = 0 first). Does it perform well? Try it on the equation 4  (x^3) - 8  (x^2) - 3  x + 9 = 0.  Why do you think its performance is less impressive here? (NB this is quite a hard question!)

4.	Mathematica, you'll be pleased to learn, has a really useful built-in command, FindRoot, that does Newton's Method (actually, a more sophisticated version) automatically. Try

FindRoot[(E^x) + (E^(-x)) == 2  x + 4, {x, 0.3}]

FindRoot[(E^x) + (E^(-x)) == 2  x + 4, {x, 1.3}]

FindRoot[x == Cos[x], {x, 1.0}]

FindRoot[4 (x^3) - 8 (x^2) - 3 x + 9 == 0, {x, 1.0}]

and so on. Incidentally, if you are wondering what sequence of iterates FindRoot is using, you can use the following trick of inserting a Print command as part of the" function" being evaluated (NB the ";" here) , for example:

FindRoot[Print[x]; (E^x) + (E^(-x)) == 2  x + 4, {x, 0.3}]

2. Calculating Friction Factors

As an example engineering application of the root finding techniques you've just been looking at, we are going to look at an equation that appears in fluid mechanics, which gives the value of the friction factor for flow in a pipe. Here it is (cf. Irving Shames, Mechanics of Fluids, 3rd edition, p. 364):

(1/(Sqrt[f]))=1.14-2.0log((e/D)+(9.35/((Subscript[Re, D])(Sqrt[f]))))

The root of this equation gives the friction factor, f, that we're after, and it depends on the Reynolds number, (Subscript[Re, D]), and the relative roughness, e/D, where D is the pipe diameter and e is the roughness (a property of the material from which the pipe is made). One other thing to note about this equation is that "log" here stands for logarithm to the base 10, which is expressed in Mathematica as Log[10, ...].

We'll translate this equation almost directly into Mathematica notation, except that we'll change the  variable names as follows: a single variable "RR" in place of e/D, "rey" for the Reynolds number, and "ff"  for the friction factor, to avoid possible confusion with the f[x] that was used earlier:

(1/(Sqrt[ff])) == 1.14 - 2.0  Log[10, RR + (9.35/(rey (Sqrt[ ff])))]

Exercises and explorations 2a

1.	Take RR = 0.0001, and make plots showing the graphs of 1/(Sqrt[ff]) and 1.14 - 2.0 log(...) on the same axes for values of rey = (10^4), (10^5) and (10^6). You'll need a plot range for ff up to about 0.05. Make a note of roughly where the graphs intersect in each case.

2.	Use FindRoot  to find the roots of the friction factor equation for  RR = 0.0001 and  rey = (10^4), (10^5) and (10^6). (Answers: 0.0310578, 0.0185208, 0.0134385.)

3.	(Optional) Use Newton's method to find the roots of the friction factor equation for the same values of RR and rey.

Graphing the friction factor

When you need a particular value of the friction factor for doing a calculation, using Mathematica's FindRoot (or a similar method in a different piece of software, or on a calculator) is a pretty optimal technique. But it's also very useful to have a graph of how the friction factor varies with the Reynolds number and the relative roughness. Indeed, you've probably already seen a graph looking something like this (cf. Shames, p.365):

Here, e/D (= RR)  is regarded as a parameter for the friction factor equation, that is, a number that we keep fixed whilst working out the friction factor, f, as a function of (Subscript[Re, D]), but that we want to vary to get different curves for f.

Logarithmic Plots

One of the main things to consider when making a graph of the friction factor is that the range of values for (Subscript[Re, D]) that are of interest span several orders of magnitude, between about (10^4) and (10^7). This is barely possible to represent well using an ordinary linear scale, and so in the graph above the horizontal scale has been made logarithmic (and the vertical scale also, for reasons that we'll touch upon in the exercise questions below). On a linear scale, the distance between two points represents a constant difference in value, so the distance between, say  (Subscript[Re, D])  =  10000 and (Subscript[Re, D])  =  20000 would be the same as that between  (Subscript[Re, D])  =  20000 and (Subscript[Re, D])  =  30000. On a logarithmic scale, a given distance represents a constant ratio, so in the graph above, the distance between  (Subscript[Re, D])  =  10000 and (Subscript[Re, D])  =  100 000 is the same as that between (Subscript[Re, D])  =  100 000 and (Subscript[Re, D])  =  1 000 000.

Mathematica has commands to handle all the standard types of logarithmic plots, but you won't find these if you simply launch Mathematica. Because they're for a specialised purpose (in the opinion of the designers of the system!), they are stored in an external file -- a package -- that you have to read in. There are nearly a hundred of these packages that come with Mathematica, and you can find out about them by browsing the Help System under "Add-ons". The one that contains what we're after is called Graphics`Graphics` (note: those are back-quotes `, not '), and you can load it by typing

Needs["Graphics`Graphics`"]

Exercises and explorations 2b

For small values of (Subscript[Re, D]), below about 2300, the friction factor can be calculated by a very simple relation, that is independent of the roughness:

f=(64/(Subscript[Re, D])).

Make plots of this relation using the regular Plot and three different kinds of logarithmic plots:

Plot[64/rey, {rey, 100,  2500}]LogPlot[64/rey, {rey, 100,  2500}]LogLinearPlot[64/rey, {rey, 100,  2500}]LogLogPlot[64/rey, {rey, 100,  2500}]

Work out what each of these plotting commands is doing -- you may find the option GridLines->Automatic useful. On one of these plots, the relation appears as a straight line, can you explain why?

Constructing the graphs

To build graphs like the one shown above, we'll use FindRoot to solve the friction factor equation, for example:

\(FindRoot[
	
      Evaluate[\((1\/\@ff == 1.14 - 2.0  Log[10, RR + 9.35\/\(rey \@ ff\)])
            \) /. {RR -> 0.0001, rey -> 10\^4}], {ff, 0.02}]\)

Of course, we need to do this repeatedly, for different choices of RR, and for (Subscript[Re, D]) varying between (10^4)and (10^7). As in previous sessions, we will call up the Table command, and construct a list of {rey,ff} coordinate pairs. For example:

\(frictionFactors = 
      Table[{rey, 
			
          ff /. FindRoot[
              Evaluate[
                \((1\/\@ff == 1.14 - 2.0  Log[10, RR + 9.35\/\(rey \@ ff\)])
                    \) /. {RR -> 0.0001}], {ff, 0.02}]}, {rey, 10\^4, 10\^5, 
          0.5*10\^4}]\)

(Since FindRoot produces a rule as its output, we need to do the "ff /. ..." , like you saw last session, to extract the ff coordinate.) We will plot lists like this using a special log-log version of ListPlot:

LogLogListPlot[frictionFactors, PlotJoined -> True, 
      GridLines -> Automatic]

So far, so good. But a problem may have occurred to you: we want (Subscript[Re, D]) to range over three orders of magnitude, but Table generates a linear sequence of sampling values. You may want to think about this for a bit before moving on to our solution.

If we want the coordinate points to be linearly spaced in (Subscript[Re, D]) when plotted logarithmically, then (since exponentiation is the inverse of taking a logarithm) we can construct our coordinate pairs using an exponential sequence, like this:

Table[(10^a), {a, 4, 7, 0.1}]

And in full, we put (10^a) in place of rey:

\(frictionFactors = 
      Table[{10\^a, 
			
          ff /. FindRoot[
              Evaluate[
                \((1\/\@ff == 
                      1.14 - 2.0  Log[10, RR + 9.35\/\(\(10\^a\) \@ff\)])
                    \) /. {RR -> 0.0001}], {ff, 0.02}]}, {a, 4, 7, 0.1}]\)

Exercises and explorations 2c

Construct a graph showing friction factor curves for (Subscript[Re, D]) ranging from (10^4) to (10^7), and for the relative roughness varying from 0.0001 to about 0.01. 

	There are many "cosmetic" options that you can apply -- consider, for example: Frame, FrameLabel, FrameTicks and StyleForm (to change the fonts). Labelling each curve with its RR value is unfortunately not so straightforward, but it can be done with FrameTicks.