Departments of Chemical Engineering and Mathematics
Mathematica Session 3

Symbolic Computation in Mathematica

1. Lists

In previous sessions, you've seen how the different kinds of brackets are used in Mathematica. You'll recall that square brackets are used to enclose the inputs to functions, as in

Sin[π]

or

D[(x^2), x]

whereas round brackets (parentheses) are used to group items algebraically, as in

3*(1 - 5)

or

(x/((x + 1) ((x^2) + 1)))

But you must have noticed that there's a third kind of bracket used in Mathematica, namely the curly bracket or brace. These crop up, for example, in graphics commands such as

Plot[{Sin[x],  Cos[x]},  {x, -π, π}]

Curly brackets are used when we have several pieces of data we'd like to keep in one place and treat as a single entity. A collection of data separated by commas and enclosed within curly brackets is called a list.

Defining your own lists is easy. You can, for example, type them in full, like this

oddList  =  {1, 3, 5, 7, 9, 11, 13, 15, 17}

Alternatively, if (as here) the list elements correspond to a rule of some kind, the command Table can be used, like this:

oddList  =  Table[2  n + 1, {n, 0, 8}]

We can pick out, say, the 5th element in oddList by typing

oddList[[5]]

(Note the double square brackets.) We can generate a list containing the elements of oddList in reverse order by typing

Reverse[oddList]

We can add new elements to lists by using Append or Prepend, as in:

Append[oddList, 19]

Prepend[oddList, -1]

Lists can be joined together:

Join[{-5, -3, -1}, oddList, {19, 21, 23}]

And note that we can assign the output of any command like this to a new variable:

newOddList = Join[{-5, -3, -1}, oddList, {19, 21, 23}]

Exercises and explorations 1

1.	Lists are fundamental to the way Mathematica works, and there are many commands to manipulate them. Work out what the following do:

Length[oddList]

Delete[oddList,  4]

Insert[oddList,  100,  3]

ReplacePart[oddList,  100,  3]

First[oddList]

Last[oddList]

2.	The command ListPlot is used to generate plots from lists of data. Try

ListPlot[{4, 9, 16, 25}]

ListPlot[{{2, 4}, {3, 9}, {4, 16}, {5, 25}}]

What's the difference between these two plots and why?

Test the effect of setting the PlotJoined option to True.

ListPlot[{4, 9, 16, 25}, PlotJoined -> True]

What other graphics options work with ListPlot, and what do they do?

3.	Type in the following list, which represents the results of an experiment:

resultsData  = 
       {{0.72, 0.35}, {1.42, 0.40}, {2.23, 0.47}, {2.91, 0.88}, {3.68, 
          0.61}, {4.40, 0.74}}

Produce a "scatter plot" (i.e. an x-y plot) of resultsData. Use list manipulation commands to produce a plot in which any obvious "rogue" points (points that clearly seem to disobey the rule obeyed by the others) have been discarded.

2. Equations, assignments and rules

Equations

You set up equations in Mathematica using the double equals sign, ==, as in

Solve[(x^2) + 4  x + 3 == 0, x]

Solve[-(x^2) == 4  x + 3, x]

Solve[{0 == 4  x + y, 3  x + ((11  y)/2) == 7}, {x, y}]

Assignments

As you have seen in previous sessions, the single equals sign has a different meaning. Roughly speaking, it's used for instructions rather than statements. The command

z = 3

doesn't mean "z equals 3" but "let z equal 3": the value 3 is assigned to the variable z.

The value doesn't have to be a number; it can be any Mathematica expression. Try

int1 = Integrate[3  (x^2), x]D[int1,  x]

or

equation  =  (x^2) - 7  x + 6 == 0Solve[equation,  x]

You can undo assignments using the Clear command:

Clear[z, int1,  equation]

It's a good habit to clear assignments to common variable names (x, y, z, a, b, etc.) when you have finished with them, so as to avoid any possible conflicts if you use them again later on in a Mathematica session.

Rules

Notice the form of the output from the previous Solve commands. You might expect the solutions of the equation (x^2) - 7  x + 6 = 0 to be expressed in Mathematica something like this:

But the actual form of the solution is different, in two ways. First, there seem to be too many curly brackets: we have a "list of lists", rather than a simple list. That's the case for technical reasons that needn't detain us here. More important, in this context, is the fact that the solutions don't appear as equations: instead of the double equals sign, we have the arrow symbol, "->"  (which is equivalent to  "->")  .

Expressions containing this arrow symbol are called rules. We have used these often in previous sessions, along with the symbol "slash-dot" (/.), to perform substitutions, as in (from Session 1)

P = R T/V

P /. V -> 0.1

To see why it might be useful for the solutions of equations to appear as rules, consider the following example:

y  =  2  (x^3) - 3  (x^2) - 12  x + 5

yprime  =  D[y, x]

solutions = Solve[yprime == 0, x]

stationaryPoints = {x, y} /. solutions

Clear[y, yprime, solutions, stationaryPoints]

Exercises and explorations 2

1.	Try the following uses of Solve:

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

Solve[{x + 2  y == 11, x*y == 12}, {x, y}]

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

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

For the last two, try using NSolve instead of Solve. What do you think NSolve does?

2.	Now try

Solve[(E^x) + (E^(-x)) == 2  x + 4, x]

Does using NSolve help here? (For something that does work on equations like this, see next session's notes.)

3. Differentiation in Mathematica

In this section, we draw together some of what you've seen about symbolic computation to explore one particular mathematical topic: differentiation. First, the basics. Recall that if you have an algebraic expression that you want to differentiate you can use Mathematica's D command, as in

D[Cos[(x^2)], x]

Where more than one variable is involved, D does partial derivatives:

D[x y, x]D[x (y^2) +  x Cos[y],  y]

If you've defined your own Mathematica function and you want to differentiate it, you have two options. Let's say, for example, that we define a function f as follows:

f[x_] := Cos[(x^2)]

To find its derivative we can type either

D[f[x], x]

or simply

f'[x]

Stationary Points on Graphs

The stationary points of a graph are the points at which the gradient is zero. Suppose we're setting out to find the stationary points on the graph of y = 2  (x^3) - 3  (x^2) - 12  x + 5. Let's start by simply differentiating this expression:

expr = 2  (x^3) - 3  (x^2) - 12  x + 5

D[expr, x]

This gives us a quadratic that is equal to zero at the stationary points, which gives us a way of finding the x-coordinates at those points: we just solve the equation "derivative = 0".

Solve[D[expr, x] == 0, x]

But what we actually need, of course, are the x- and y-coordinates. We can get those in one go by using the fact that the solutions of equations come in the form of rules:

{x, expr} /. Solve[D[expr, x] == 0, x]

We can show these on a plot by typing

plot1 = ListPlot[%, PlotStyle -> PointSize[0.03]]

(Recall that "%" refers to the last output.)

The following command gives us a graph of y = 2  (x^3) - 3  (x^2) - 12  x + 5:

plot2 = Plot[expr, {x, -3, 3}, PlotStyle -> RGBColor[1, 0, 0]]

Finally, to show the function and its stationary points on the same figure:

Show[plot1, plot2]

To be on the safe side, we should now type

Clear[expr]

Exercises and explorations 3a

1.	Use the above technique to find and display the stationary points of the degree-7 polynomial

\(TraditionalForm
      \`y = 2  x\^7 - 7  x\^6 + 21  x\^4 - 14  x\^3 - 21  x\^2 + 28  x\).

2.	Stationary points come in two types: extrema (that is, maxima and minima) and points of horizontal inflexion. By examining your graph, or otherwise, identify to which of these two categories each of the stationary points of the above polynomial belongs.

3.	Could this have been predicted without plotting any graphs? In general, how do you think you could tell, going only on the output from Solve and without plotting, whether a particular stationary point is an extremum or a point of horizontal inflexion? Test your idea, for example on the non-polynomial function

y=7(x^2)-42x-(84/x)+(21/(x^2))+(14/(x^3))-(7/(x^4)).

What's going on here mathematically?

Stationary Points on Surfaces

On a line graph, a stationary point occurs wherever the derivative is zero. The most commonly encountered type of stationary point is the extremum, that is the maxiumum or minimum; points of horizontal inflexion are in practice somewhat rarer.

On a surface, a stationary point occurs whenever both partial derivatives are zero. There are two commonly encountered types: the extremum and the saddle, along with several somewhat rarer types that we don't study in depth.

Exercises and explorations 3b

1.	Consider the function y=(((x^2)-(y^2))/((((1 + (x^2) + (y^2)))^2))). Generate a surface plot of this function by typing

expr = (((x^2) - (y^2))/(((1 + (x^2) + (y^2)))^2))

Plot3D[expr,  {x, -2, 2}, {y, -2, 2}]

2.	 Examine the picture on your screen. How many points can you see at which the surface appears to be locally horizontal and where, roughly, are they located? Which of these are maxima and which minima? Are there any others, and if so how would you describe it/them?

3.	 Calculate the positions of the stationary points by typing

Solve[{D[expr, x] == 0, D[expr, y] == 0}, {x, y}]

Were you right about their locations?

4.	 Generate a contour diagram of this surface by typing

ContourPlot[expr,  {x, -2, 2}, {y, -2, 2}, PlotPoints -> 60]

How do maxima, minima and saddles show up in such diagrams?

5.	 Use the above techniques to investigate the function y=((x((x^2) - (y^2)))/((((1 + (x^2) + (y^2)))^2))). How would you describe the stationary point at the origin? Does it fit into our classification? Why do you think this type of stationary point rarely arises in practical applications?

Taylor Series

Among the easiest mathematical functions of all to deal with are polynomials: functions of the form

\(TraditionalForm
      \`f(x) = a\_0 + \(a\_1\) x + \(a\_2\) x\^2 + \(a\_3\) x\^3 + \ 
          ⋯\  + \(a\_n\) x\^n\).

Unfortunately, most functions aren't polynomials. Fortunately, though, a great many (though not all) of them can be approximated, at least for certain values of x, by polynomials.

For example, consider the function y = (e^x). Close to x = 0, this can be approximated reasonably well by the straight line y = 1 + x, as we can see from the following graph:

Plot[{(E^x), 1 + x}, {x, -3, 3}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1]}]

(the exponential function appears in red, and the linear one in blue).

We can get a slightly better approximation, that holds over a slightly wider range, by using the quadratic expression 1 + x + (1/2) (x^2) in place of 1 + x, as follows

Plot[{(E^x), 1 + x + ((x^2)/2)}, {x, -3, 3}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1]}]

We can inprove things still further by adding a cubic term:

Plot[{(E^x), 1 + x + ((x^2)/2) + ((x^3)/6)}, {x, -3, 3}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1]}]

The successive polynomial approximations come from McLaurin's Theorem, which says that if all the derivatives of f are defined at x = 0 then, for some set of x-values including 0,

f(x)=f(0)+x f ' (0)+((x^2)/(2!))f ''(0)+⋯+((x^n)/(n!))(f^( (n)))(0)+⋯.

Exercises and explorations 3c

1.	Polynomial approximations are obtainable in Mathematica using the Series command:

Series[(E^x), {x, 0, 2}]

Series[(E^x), {x, 0, 3}]

and so on.

To remove the final "order" term and convert this truncated series into an ordinary polynomial, use Normal:

Normal[Series[(E^x), {x, 0, 2}]]

\(expr = Normal[Series[E\^x, {x, 0, 3}]]; 

    Plot[{E\^x, expr}, {x, \(-3\), 3}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1]}]\)

Generate graphs of the quartic and quintic (order 4 and 5) McLaurin approximations to (e^x). Does the increase in accuracy and range continue?

2.	 Repeat, but this time study arctan x instead (in Mathematica syntax, ArcTan[x]). What's similar and what's different about what you observe?

3.	 What do you expect to see if you type

Series[(E^x), {x, 1, 3}]

instead of

Series[(E^x), {x, 0, 3}]

Try it. Generate successive approximations using a 1 instead of a 0 in this middle position. Plot each of them on the same axes as (e^x). What do you see? Repeat for arctan x.What are the similarities and differences?