Departments of Chemical Engineering and Mathematics
Mathematica Session 5

Building your own functions

1. Defining Simple Functions

You have already seen some simple function definitions in previous sessions; for example this (stream) function in Session 2 for flow past a rotating cylinder:

ψ2[r_, θ_] := 
      v (r - ((a^2)/r)) Sin[θ] + ω (a^2) Log[r]

From the moment you evaluate this piece of code until the end of your session -- or until the system crashes, whichever is sooner! --  the function ψ2 can be used in combination with all Mathematica's built-in functions and commands.

You have seen that the combined symbol ":= "  ("colon equals") is used for defining functions, where you might have expected the simple  "=" sign, for assignment.  The difference between = and := in Mathematica is quite subtle: for the moment it's worth learning the rule of thumb that := is almost always best for defining functions, whereas = is almost always best for assigning values to variables, e.g.  z = 3.

The other peculiarity about function definitions is the use of the underscore symbol, _, on the left-hand side.  The following exercises demonstrates why underscores are essential.

Exercises and Explorations 1

Define a function quad by typing

quad[x_] := (x^2) - 4

Test your definition by typing stuff like

quad[3]

quad[x]

quad[x + 2]

quad[y]

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

Expand[quad[x + 2]]

D[quad[x], x]

Integrate[quad[t], t]

Now clear your definition by typing

Clear[quad]

and try the whole thing again, but this time without the underscore: that is, starting with

quad[x] := (x^2) - 4

What goes wrong?  How would you describe the role of the underscore symbol?

One final remark on function definitions. It is a good habit to always begin function names with a lower-case letter, even though Mathematica's built-in functions all begin with capitals.  Well, actually, precisely because they do.  The last thing you want is to define a new function and accidentally give it a name that's already been used by Mathematica, and if you stick to lower-case letters you avoid that.

2. Building New Functions

When we speak of a function in mathematics, we usually  mean a rule linking numbers to other numbers.  In Mathematica, though, and in computing generally, the term function is used to describe anything with a collection of inputs and outputs, whatever form those inputs and outputs take.

For example, Mathematica's Integrate facility is a function in this sense.  It takes input in the form of an algebraic expression and a variable name, as in:

Integrate[2  x Sin[(x^2)],  x]

and returns output in the form of another algebraic expression, namely -Cos[(x^2)].

Solve is also a function.  It takes input in the form of a set of equations and a set of variables to be solved for, and outputs the solutions in a special format as Mathematica "rules",

Solve[{x  +  2  y  ==  5,  x y  ==  2},  {x,  y}]

By now it may have occurred to you that, since nearly every use of Mathematica involves inputs and outputs, nearly all of what Mathematica does can be described in terms of functions. In graphics, too, the appearance of a graph -- e.g. the appearance of lines (thick/thin, solid/dashed) and colours -- is often adjusted by way of functions. Recall the RGBColor function from previous sessions, whose inputs are the (relative) magnitudes of red, green and blue, and whose output is, in effect, the resultant colour.

Since functions are so basic to Mathematica's workings, the idea of function-building can be carried across to virtually all contexts in which you want to get Mathematica to do something new, whether or not you're working with numbers.

Worked example: Critical contours

Here's an example. Consider the function \(TraditionalForm
      \`f(x, y) = x\^2 + y\^2 - 7\ x\ y + 10\ x\^3\ y + 10\ x\ y\^3\). Let's have a look at its contour plot:

cp1 = 
      ContourPlot[
        (x^2) + (y^2) - 7 x y + 10 (x^3) y + 10 x (y^3), {x, -1.2, 
          1.2}, {y, -1.2, 1.2}, ContourShading -> False, PlotPoints -> 60]

Not very informative. Let's try increasing the number of contours:

cp2 = 
      ContourPlot[
        (x^2) + (y^2) - 7 x y + 10 (x^3) y + 10 x (y^3), {x, -1.2, 
          1.2}, {y, -1.2, 1.2}, ContourShading -> False, 
        PlotPoints -> 60, Contours -> 20]

Better, but not much. The problem is that the really important contours -- the ones that give us the most useful information -- aren't automatically plotted by Mathematica. They're the so-called critical contours: the ones that pass through any saddle points the function may possess.

So, let's set ourselves the task of building a Mathematica function that finds these critical contours and plots them. As with all problem-solving exercises, it's a good idea if we work systematically, breaking down the problem into smaller sub-problems. In this case, such a plan of work might look like this:

1. Set up a set of equations for finding the stationary points of a 2-variable function.
2. Solve these equations.
3. Devise a way of deciding which of the solutions represent saddle points.
4. Use this decision method to select the saddle points only.
5. Find the values of the function at these saddle points.
6. Plot the contours corresponding to those values.

We can implement each of these stages as a Mathematica function.  Stage 1 to begin with:

statPtEqns[expr_, x_, y_] := {D[expr, x] == 0, D[expr, y] == 0}

Testing this on our example:

statPtEqns[(x^2) + (y^2) - 7 x y + 10 (x^3) y + 10 x (y^3), x, y]

Now stage 2, in which we actually solve the equations:

statPtSolutions[expr_, x_, y_] := NSolve[statPtEqns[expr, x, y]]

Testing again

statPtSolutions[(x^2) + (y^2) - 7 x y + 10 (x^3) y + 10 x (y^3), x, y]

Now for the hard bit, mathematically speaking. You have not met this in your maths course yet, but it can be shown that a stationary point of f is a saddle provided

(Subscript[f, xx])(Subscript[f, yy])-((Subscript[f, xy])^2)<0,

where (Subscript[f, xx]) is the second partial derivative of f with respect to x, etc. This gives us

saddleTest[expr_, x_, y_][solution_] := 
      TrueQ[(D[expr, {x, 2}] D[expr, {y, 2}] - D[expr, x, y]^2 /. solution)
             < 0]

This uses the TrueQ command, which is new to you. It outputs True if the condition inside is definitely true, and False otherwise, as the following tests suggest:

saddleTest[(x^2) + (y^2) - 7 x y + 10 (x^3) y + 10 x (y^3), x, y]
        [{x -> -0.827895, y -> 0.120788}]

saddleTest[(x^2) + (y^2) - 7 x y + 10 (x^3) y + 10 x (y^3), x, y]
        [{x -> -0.474342, y -> 0.474342}]

Stage 4, selecting the saddle points:

saddleSolutions[expr_, x_, y_] := 
      Select[statPtSolutions[expr, x, y], saddleTest[expr, x, y]]

and the test:

saddleSolutions[(x^2) + (y^2) - 7 x y + 10 (x^3) y + 10 x (y^3), x, y]

Stage 5, finding the associated contour values, might work like this:

criticalValues[expr_, x_, y_] := 
      expr /. saddleSolutions[expr, x, y]

Testing...

criticalValues[(x^2) + (y^2) - 7 x y + 10 (x^3) y + 10 x (y^3), x, y]

... reveals that we have some redundant information here. That's not fatal, just untidy, so we won't worry about it. The final stage is plotting the contours:

\(criticalContourPlot[expr_, {x_, xmin_, xmax_}, {y_, ymin_, ymax_}] := 
      ContourPlot[expr, {x, xmin, xmax}, {y, ymin, ymax}, 
        Contours -> criticalValues[expr, x, y], PlotPoints -> 60, 
		
        ContourStyle -> RGBColor[1, 0, 0], ContourShading -> False]\)

Testing:

ccp = 
      criticalContourPlot[
        (x^2) + (y^2) - 7 x y + 10 (x^3) y + 10 x (y^3), {x, -1.2, 
          1.2}, {y, -1.2, 1.2}]

Finally, we can show this plot on the same pair of axes as the "ordinary" contour plot we generated earlier:

Show[{cp2, ccp}]

Look carefully at this way of working: we've broken down our task into simple pieces, each dependent on the previous ones, and represented each by a simple Mathematica function, each of which calls one or more of the previous ones.

Exercises and Explorations 2

1.	Test the criticalContourPlot function on some other expressions. Try, for instance

criticalContourPlot[
      -1 + 2 (x^2) - (x^4) + (y^2), {x, -2, 2}, {y, -2, 2}]

criticalContourPlot[
      -1 + x + 2 (x^2) - (x^4) + (y^2), {x, -2, 2}, {y, -2, 2}]

criticalContourPlot[
      -1 + x - (x^4) + (y^2), {x, -2, 2}, {y, -2, 2}]

plus some of your own. What error messages do you get? Can you get the code to fail? How might you improve the performance of criticalContourPlot?

2.	Your main task for this section is to write and test a function that generates, for a given expression expr in one variable, a graph of expr showing each of its stationary points. You've already done this for particular expressions in Session 3: the task here is to make a function to do the job for any expression you feed it (any polynomial expression, anyway).

Plan your work carefully. Look back at the notes from Session 3 to remind yourself of the particular Mathematica features you're going to need to make use of. Break down the work into simple sub-problems, and try to devise Mathematica functions to solve each of these in turn. Test as you go along.

A prize for the most elegant solution we see? You never know.

3. Particle Paths for Fluid Flows (Optional)

In this final section, we want to offer a sophisticated example of function building, in the context of fluid mechanics.  In Session 2, we looked at the streamlines for some simple fluid flows, by plotting the contours for appropriate stream functions -- this being the simplest, for flow past a cylinder:

ψ(r,θ)=v (r-((a^2)/r))sin θ  .

To save typing later, let's set this up as a direct function of x and y:

ψ1[x_, y_] := 
      v ((Sqrt[(x^2) + (y^2)]) - ((a^2)/(Sqrt[(x^2) + (y^2)]))) 
        Sin[ArcTan[x, y]]

And then, for example,

streamlines1 = 
      ContourPlot[
        ψ1[x, y] /. {v -> 10, a -> 0.5}, {x, -2, 2}, {y, -2, 2}, 
        PlotPoints -> 30, ContourShading -> False, 
        Epilog -> {GrayLevel[0.5], Disk[{0, 0}, 0.5]}]

The velocity of a (small) particle moving with the flow is given, in terms of its x and y components, by the following partial derivatives of the stream function, ψ:

(Subscript[v, x])= ((∂ψ)/(∂y)),  (Subscript[v, y])=- ((∂ψ)/(∂x)).

Our plan is as follows:

1. Create a function that calculates these velocity components given any stream function, ψ.
2. Create a function that solves the simultaneous differential equations x'(t) = (Subscript[v, x]) (t) and  y'(t) = (Subscript[v, y])(t) numerically, so that we have a solution to describe a particle's x and y coordinates as a function of time.
3. Create functions to visualise a solution in different ways: as a static graph showing the path  followed by a particle over time, and as a  movie showing a particle moving in "real time" in the flow.

Stage 1: Recall that the D command gives us partial derivatives, for example

-D[ψ1[x, y], x]

This looks pretty gruesome, but Mathematica has no trouble with it. As an added complication, though, we have the fact that Mathematica's numerical solver needs its expressions to be in terms of x[t] and y[t] instead of just x and y. This can be quite easily achieved, as follows.

-D[ψ1[x, y], x] /. {x -> x[t], y -> y[t]}

We're now ready to write our function:

\(velocityEquations[streamFunction_] := {
		
        \(x'\)[t] == \((D[streamFunction, y] /. {x -> x[t], y -> y[t]})\), 

        		\(y'\)[t] == 
          \((\(-D[streamFunction, x]\) /. {x -> x[t], y -> y[t]})\)
}\)

Let's set up a particular stream function and test our new Mathematica function:

sf1 = ψ1[x, y] /. {v -> 10, a -> 0.5}

velocityEquations[sf1]

Stage 2: Mathematica's numerical differential equation solver is called NDSolve. We must give it the differential equations to be solved, a pair of initial conditions (i.e. the starting coordinates for x and y, such as -2.0 and 0.25), and the time range required for the solution (e.g. from t = 0 to t = 0.5):

NDSolve[
      Join[velocityEquations[sf1], {x[0] == -2.0, y[0] == 0.25}], {x[t], 
        y[t]}, {t, 0, 0.5}]

(Note the use of Join to combine two lists here -- cf. Session 3, Section 1.) The output of NDSolve is in the (rather baffling) form of interpolating functions. These are strange objects: they contain, in a hidden form, all the information we need about our solution. You'll meet one way of getting at that information below.

Annoyingly, there's an extra pair of curly brackets; you can get rid of these by using First:

First[
      NDSolve[Join[
          velocityEquations[sf1], {x[0] == -2.0, y[0] == 0.25}], {x[t], 
          y[t]}, {t, 0, 0.5}]]

Just as when we solve ordinary algebraic equations, our solutions arrive from NDSolve in the form of rules. To get things that we can operate on directly, therefore, we need to use the "/." operation:

{x[t], y[t]} /. 
      First[NDSolve[
          Join[velocityEquations[sf1], {x[0] == -2.0, y[0] == 0.25}], {
            x[t], y[t]}, {t, 0, 0.5}]]

We're now ready to write our solving function, where we generalise the above to any given initial point, {x0, y0}, and any given time range (t = tmin to t = tmax):

pathSolution[
        streamFunction_, {x0_, y0_}, {t_, tmin_, tmax_}] := {x[t], y[t]} /. 
        First[NDSolve[
            Join[velocityEquations[streamFunction], {x[0] == x0, 
                y[0] == y0}], {x[t], y[t]}, {t, tmin, tmax}]]

Let's test this function:

path1 = pathSolution[sf1, {-2.0, 0.25}, {t, 0, 0.5}]

Stage 3: You can plot the x-coordinate or y-coordinate against t using Plot. If you have assigned path1 as above, try:

Plot[path1[[1]], {t, 0, 0.5}]

Plot[path1[[2]], {t, 0, 0.5}]

(The [[...]] operation for picking out elements of a list was introduced in Session 3.) The following function will produce a plot of y against x -- i.e. the particle path in the (x, y) plane:

pathPlot[pathSoln_, {t_, tmin_, tmax_}] := 
      ParametricPlot[pathSoln, {t, tmin, tmax}, 
        PlotStyle -> RGBColor[1, 0, 0]]

(Note that tmin and tmax have a different meaning here from that in the definition of pathSolution: there, they specified the range of calculation; here, they specify the plotting range -- which cannot, of course, exceed the calculated range.) Let's test this on path1:

pathplot1 = pathPlot[path1, {t, 0, 0.5}]

You don't have to plot all of the path -- you can use the t range to select a portion of it:

pathPlot[path1, {t, 0.25, 0.5}]

You can superimpose path plots with streamline plots using Show, for example:

Show[streamlines1, pathplot1, AspectRatio -> Automatic]

Finally, here is a simple function to create particle path movies. As in previous sessions, we'll use Table to build up a sequence of plots:

pathMovie[pathSoln_, {t_, tmin_, tmax_}, frames_] := 
      Table[ListPlot[{pathSoln}, PlotRange -> {{-2, 2}, {-2, 2}}, 
          PlotStyle -> {RGBColor[1, 0, 0], PointSize[0.02]}], {t, tmin, tmax, 
          (tmax - tmin)/(frames - 1)}]

Test:

pathMovie[path1, {t, 0, 0.5}, 20]

pathMovie as defined here isn't ideal -- it would be better to have the particle motion overlaid on a streamline plot, and also to be able to view several particle paths at once. Both of these require Mathematica constructions a little beyond the present discussion; if you're keen you can find some extra code to achieve these kinds of results on the course web page.

Exercises and Explorations 3

1.	Using the functions introduced in this section, investigate particle paths for the stationary cylinder flow starting to the left of the cylinder and on or near the line y = 0 (take, for example,  x0 = -2.0, y0 = 0.0, +0.01, -0.01, +0.1, -0.1, etc). Use pathMovie to observe the changing velocity of a particle as it travels. What happens as particles go near the stagnation point at x = -0.5, y = 0?  What happens if you start a particle path at this point?

2.	Investigate the related, but more complex situation for the rotating cylinder case when ω > 2 v/a (cf. Session 2).
	Calculate particle paths starting at points such as A, B, C and D in the diagram below (for which ω = 70). What happens to the particles as they go near the saddle point, X, of the "critical contour"? What would happen if you started a particle path at X? (X is located at {0.0, -1.59307} for the ω = 70 case.)