Functions

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2006

Introduction

What is a function?

In the most general sense, a function is any conceptual object that accepts input and produces corresponding output based upon that input.  Of course, ordinary mathematical functions, such as sine or cosine, obviously satisfy this definition, as would vector-valued functions which produce output in the form of vectors rather than simple numbers.  More sophisticated functions might receive input in the form of a symbolic expression and return output in the form of a transformed expression or a plot.  Perhaps the output might consist of symbolic, numerical, and graphical parts bundled together in a list.  The output might consist of a list of possible solutions to a problem rather than one specific result.  The output could even be an entire Mathematica notebook!

Functions lie at the heart of Mathematica — any operation it performs can be represented by a function which, when presented with appropriate input, returns one of the many output objects it recognizes.  Therefore, most of our work with Mathematica will consist of constructing functions which perform desirable tasks.

In this notebook we introduce many of the most important concepts and techniques used to create Mathematica functions; however, the richness of the Mathematica language precludes an exhaustive treatment. Moreover, this notebook is rather long and many readers may find plodding through syntax to be exhausting in itself.  Therefore, I encourage you to skip around through the notebook, learning the basics now and returning for more detail later as needed.  Short and simple examples appear throughout the notebook and links to more extended and more interesting examples are provided.  Feel free to branch out to these examples frequently — although some of these examples employ features of Mathematica that may not have been covered yet, people do not learn in the linear fashion of a traditional textbook.  One tends to learn best only those tools that are needed to solve the problem at hand, and part of solving problems is discovering what tools are available that might be helpful.  Only after a tool or technique has been recognized as potentially useful to a current problem is the motivation to master it strong enough to promote effective learning.

Acknowledgments

Prof. MacDonald (deceased) introduced the concepts of explicit and implicit functions to distinguish between functions with explicit arguments and named expressions that are used as functions without explicit arguments.  He also coined the term skeletal functions to describe unnamed pure functions whose arguments are given as slots rather than named variables.  I retain this descriptive nomenclature here.

Initialization

ClearAll["Global`*"];
Off[General::spell,  General::spell1]

Implicit functions

An implicit function is simply a named expression which contains one or more symbols that are considered variables.  For example, below we define f to be a quadratic function of x, where a, b, and c are considered parameters, although here the distinction between variables and parameters is merely a matter of attitude.

Clear[x, a, b, c];
f  =  a (x^2)  +  b x  +  c;

As presently defined, f is an implicit function which simply returns its symbolic definition.

f

Numerical values can be obtained either by using replacement rules to specify the values of the variables and parameters

f  /. {a -> 1, b -> 2, c -> 3, x -> 4}

which does not alter the definition of the function

f

or by assigning values to the parameters

a = 1;  b = 2;  c = 3;
f

and variables

x = 4;
f

However, assignment of values precludes many of the operations we might wish to perform on such a function.  For example, we can no longer differentiate this function with respect to x

D[f, x]

because x is now a number instead of a variable.

?x

Although we can undo this unfortunate assignment

Clear[x];
D[f, x]

you will soon find that unintended or obsolete assignments can be a major source of frustration to new users of Mathematica.  For example, if we neglect to clear the variables used on the rhs and forget that a and/or b have already been assigned values, a new function

g  =  a Sin[x]  +  b Cos[x]

does not contain the parameters a and b as intended.  A proper definition would require two new symbols, such as d and e, but symbol creation consumes memory (both computer and human) and slows down evaluations.

Clear[a, b, c, x]

A good technique for assigning parameters to a function is to give a set of replacement rules a name that clearly identifies it purpose and its association with the relevant function.  For example, the coefficients of the quadratic function called f above are identified below as its parameters.  Note that we employ a string variable to avoid evaluation of the expression for f.

parameters["f"] = {a -> 1, b -> 2, 
          c -> 3};

f /. parameters["f"]

f /. x -> 2 /. parameters["f"]

This technique enhances the readability of notebooks and reduces the chance of errors due to spurious assignments.

Express f = (1/(1  +  (1/(1 + (z^2))))) as an implicit function of z.  Plot f for z -> sin(x).

Express f = (1/(1  +  (1/(1 + (z^2))))) as an implicit function of z.  Evaluate ∂(f/∂)x for z -> sin(x).

Explicit functions

Explicit functions receive their variables as arguments.  The definitions can be made using either immediate or delayed assignments, as for implicit functions, but we recommend delayed assignment to avoid conflicts with pre-existing symbols.  The arguments can be expressed as symbols, patterns, or a mixture of types.  The functions may be represented in terms of several rules, with the rule that applies to a particular set of arguments chosen by pattern matching.

Patterns as arguments

Often the most convenient method for specifying the arguments of an explicit function is to use patterns rather than explicit symbols.  Consider a function

Clear[g]

g[x_] = (x^2) Exp[-(x^2)];

whose argument is expressed in terms of the pattern x_.   The single underscore (Blank) is the symbol for a blank pattern that may assume any value, symbolic or numeric.  Attaching the blank pattern to the symbol x in the form of x_ indicates that the pattern will be referred to as x on the right-hand side (rhs) of the definition statement.  Note that this usage of the symbol x does not conflict with any global definitions that might already exist.  Thus,

g[x]

gives the expected result while

g[y]

now gives the corresponding expression with x replaced by y.  Similarly,

g[2]

now evaluates to an exact numerical expression and

g[2. ]

evaluates to an approximate number.  There is also nothing to stop us from using a string as the argument,

g["stuff"]

g["stuff"] // FullForm

even if the result is uninterpretable as a mathematical function.  It is also possible to use a function as the argument.

g[g[z]]

This technique can be very useful for recursive functions, such as iterated function systems, or for symbolic manipulation and will be explored at greater depth later.

Similarly, functions of several variables can be constructed using different symbols to identify each pattern.

g[x_, y_] =  (x^(2 + y)) Exp[-(x^2)];

{g[x, y], g[x, 0], g[x, 2], g[x, -4], g[2, 2], 
      g[x, "stuff"]}

Each of these forms evaluates in the expected manner.  However, why didn't we use a different symbol for this new function?  Will a second definition for the same symbol, namely g, conflict with the first?

g[z]

Evidently, there is no conflict between these definitions!  Harmonious coexistence arises from the method that Mathematica uses to interpret expressions.  Mathematica stores the rules associated with each symbol and then, when evaluating expressions, scans its rule tables for any rules associated with the symbols that appear in that expression, beginning with those defined by the user and proceeding to the built-in rules that apply to any symbol of the same type.  The appropriate rule is selected by pattern matching.  The rules associated with the symbol g can be inspected using the ? operator; sometimes ?? will give even more information.

?g

Here we find two rules associated with the symbol g.  If our expression matches the pattern g[x_] the rule for that pattern is applied, or if our expression matches g[x_,y_] the corresponding rule is applied.  There is no conflict because the patterns are unique.  On the other hand, no evaluation is performed when the expression involves an unrecognized pattern.

g[x, y, z]

Therefore,  this technique can be used to define generalized functions which actually consist of several different functional relationships that are assigned the same symbolic name by virtue of their conceptual relationship but which apply to different input patterns.  Of course, care should be exercised in the design of functions with multiple definitions to avoid overlapping or ambiguous patterns.

Explain the behavior of the following functions: g[x_]=π or h[_]=π.  Are they equivalent?

Explain the behavior of the following function:
	Clear[f];
	f[x_] = (x^2);
	f[x_,y_] = (x^y);
Evaluate f[2], f[3,2], f[x,(y/Log)[x]].  How do these several definitions coexist peacefully?

Explain the behavior of the following function:
	area["rectangle",b_,h_] = b h;
	area["triangle",b_,h_] = b h /2
Add similar definitions for parallelograms and trapezoids (you may need more arguments!) and verify that the new definitions do not conflict with the previous definitions.  Why is it better to place the symbolic arguments before the patterns?

Delayed assignment

Immediate assignment in the form lhs=rhs is appropriate when rhs provides an immediate result that can be delivered without further processing of the lhs.  However, if rhs is more complicated, such as the result of a plotting command, that cannot be evaluated directly, then we should use delayed assignment of the form lhs:=rhs.  Consider the following function intended to produce a graph of the derivative of a function of a single variable over a specified domain.

PlotDerivative[f_, {x_, xmin_, xmax_}] := 
      Plot[Evaluate[D[f[x], x]], {x, xmin, xmax}]

An example of the use of this function is given below.

PlotDerivative[AiryAi, {x, -10, 4}];

Attempts to define this function using immediate assignment fail because neither the derivative nor the graph can be constructed before a specified choice of f is available.

Clear[f]

BadPlotDerivative[f_, {x_, xmin_, xmax_}] = 
      Plot[Evaluate[D[f[x], x]], {x, xmin, xmax}]

Another advantage of delayed assignment is that the risk of conflict with pre-existing symbols is reduced.  Without clearing f, a different, but equally useless, result is obtained from the failed function above if an earlier definition of f[x] is present.  The definition using delayed assignment is immune to this problem.

Example: spherical Bessel function

For some unknown reason, Mathematica lacks built-in definitions for spherical Bessel functions.  These functions are defined as

(Subscript[𝒿, ℓ])(
        z) = (Sqrt[(π/2)]) (z^(-1/
              2)) (Subscript[J, ℓ + (1/2)])(z)

(Subscript[y, ℓ])(
        z) = (Sqrt[(π/2)]) (z^(-1/
              2)) (Subscript[Y, ℓ + (1/2)])(z)

where (Subscript[J, n]) is the regular and (Subscript[Y, n]) is the irregular cylinder function.  We can define our own spherical Bessel functions, without fancy notation, as follows.

sbessj[n_, z_]  :=  (Sqrt[(π/2)]) 
          (z^(-1/2)) BesselJ[n + (1/2), z];
sbessy[n_, z_]  :=  (Sqrt[(π/2)]) 
          (z^(-1/2)) BesselY[n + (1/2), z];

We can also define a function which plots these functions for given order over a specified range.  For detailed information about plotting commands, see plotting.nb.  For the present purposes it is sufficient to recognize that plotsbess is defined as an explicit function of three variables which specify the order of the Bessel function and the minimum and maximum arguments to be plotted.  This function uses the built-in Mathematica function Plot to produce a figure that displays both the regular and irregular spherical Bessel functions over the specified range.

plotsbess[n_, xmin_, xmax_] := 
      Plot[{sbessj[n, x], sbessy[n, x]}, {x, xmin, xmax}]

plotsbess[0, 0.001, 10];

This simple definition for sbessj works for most arguments, but fails for x = 0.

sbessj[0, 0]

We will see how to correct this type of problem shortly.

exercises

The Stirling approximation for Log[n!] is \(TraditionalForm
      \`Log[\(n!\)] ≈ \ 
        n\ Log[n]\  - \ n\  + \ \(1\/2\) Log[2  π\ n]\  + \ 
          1\/\(12\ n\) + \ ⋯\).  
a) Write functions logStirling1 using just the first two terms (the lowest reasonable approximation), logStirling2 including the next term, and logStirling3 including the next two terms.
 b) Also write functions which evaluate the percentage error in these approximations.  
 c) Make a LogPlot which displays the percentage errors for all three approximations as functions of n.  [Hint: you will probably need to load the package Graphics`Graphics` to use LogPlot.]

Write a new version of the Stirling function for which logStirling[n,1] gives the first approximation, logStirling[n,2] the second, and logStirling[n,3] the third.  Furthermore, logStirling[n] with only a single argument should give the most accurate formula, namely logStirling[n,3].

Arguments in List or Sequence form

Many explicit functions receive their arguments in the form of a sequence, such as f[a,b,c], while others receive a list of arguments, such as g[{a,b,c}].  However, we often find our arguments represented in an inappropriate form and need a mechanism for converting a List to a Sequence.  For example, consider the following function defined in terms of a sequence of three arguments.

Clear[f, a, b, c, d];
f[a_, b_, c_, d_] := d (a^((b^c)))

We might as a result of earlier work have some of the desired arguments in the form of a list, but when we present that list to our function

argumentList = {a, b, c};

f[argumentList, d]

no evaluation is performed because the argument sequence does not match the pattern employed to define the function.  A convenient solution to this problem is provided by the function Sequence, which stores its arguments in the form of a sequence which can be inserted into an argument list.  Applying Sequence to the argumentList replaces the head List with the head Sequence.

argumentSequence = Sequence @@ argumentList

f[argumentSequence, d]

In fact, a sequence of this form can be inserted into any argument list and will act as if the arguments were provided directly.

Clear[g];
g[stuff, argumentSequence, moreStuff]

Functions formed from command sequences

We will often need functions which require many steps to produce their results.  Perhaps the simplest method for assembling a multistep function is to define a function in terms of an expression which contains a sequence of commands.  The basic syntax would be: 
	f[args_] := (expr1; expr2; ⋯ ; output)
where the initial steps return their output to the next step rather than to the user, with only the final step returning a printed result because it does not end with a terminal semicolon.  For example, the following function finds the coefficient of (x^n) for (((a + b x))^m) using Mathematica functions Expand and Coefficient.

Clear[f, n, m]

f[n_, m_]  :=  (temp  =  Expand[(((a + b x))^m)];  
        Coefficient[temp, x, n])

f[3, 13]

Although this function achieves its goal, it has the undesirable effect of leaving behind detritus in the form of unwanted symbols which should have been temporary.

temp

The simple example suggested above can be written more efficiently using pure functions which do not require creation of intermediate symbols, but it is often difficult to formulate more interesting multistep functions in readable fashion without creating local variables.  Traditional programming languages solve this problem by providing constructions known as subprograms to isolate the local variables which are used for intermediate steps and then are discarded when no longer needed.  These subprograms also protect their local variables from corruption by external commands.  The importance of protecting local variables is illustrated in a later section.

Pure functions

A pure function is represented by the basic syntax Function[var,body] where var represents the argument and body the definition of the function.  The argument can be a single variable or a list of variables.  The function is then applied to an argument using the syntax Function[var,body][arg] to replace var by arg in body.

Pure functions with names

If we intend to use a function several times we should assign it a name and refer to it by name rather than repeat its definition each time it is used.

h  =  Function[x, a (x^2)  +  b x  +  c]

h[Sin[ϕ]^2]

h[h[y]]

Here h is simply a symbol that represents a pure function.

?h

The formal arguments used to define the body of the function are local variables that will not interfere with assignments made elsewhere for symbols with the same names.  In the example below, the assignment given to x does not affect the definition of a function that uses x as its local variable.

x = "something";
g = Function[x, (x^2)];

g[y]

However, parameters that appear in the body of the functions are treated as global symbols.  Therefore, assignments of those parameters will affect subsequent evaluation of the functions and must be made with care.

h[a]

a = 2;
h[z]

Clear[a]

It is also important to recognize that functions defined in this manner are not evaluated until arguments are presented.  Therefore, it is not necessary to use delayed assignment to prevent premature evaluation of the right-hand side.  Notice how in the cell below the definition using Function behaves properly while similar definitions using immediate assignment fail, the first by failing to protect the local variable and the second by premature evaluation.

Clear[f1, f2, g, y, z];
f1[x_] = Expand[x];
f2[y_] = Expand[y];
g = Function[x, Expand[x]];

f1[(((y + z))^2)]

f2[(((y + z))^2)]

g[(((y + z))^2)]

We call functions of this type pure because we define the behavior of the function, not its arguments.  In fact, we can manipulate a pure function by name without supplying any arguments at all.  For example, suppose that we require the derivative of the quadratic function defined above.  The derivative can be evaluated in an obvious fashion using the built-in function D

D[h[y], y]

but this method requires that h be given an explicit argument.  We prefer the pure form using the conventional prime notation

h'

to obtain the derivative as a pure function also.  A pure function can be evaluated with any suitable argument.

{h'[stuff], h'[(2/3)]}

The logistic map is defined by the function f[x_]:=μ x (1-x) where μ is the growth parameter.  Write a pure function that embodies this mapping.  Use NestList to produce an iterated function sequence {y,f[y],f[f[y]],f[f[f[y]]],...} nested 4 times.

Pure and explicit function definitions can be mixed.  Study the behavior of the following function:
	f[μ_]=Function[x,μ x (1-x)]
Evaluate the following expressions and interpret the results: 
a) f[a]  	b) f[a][x] 	c)f[2.5][0.8]	d) NestList[f[μ],y,4]

Once again define:
	f[μ_]=Function[x,μ x (1-x)]
Then evaluate the following expressions and interpret the results: 
a) f[a] ' 	b) f[a]'[z] 	c) f[3.5]'[0.25]

Pure functions of several variables

Pure functions with several variables are defined by using a list of variable names in place of var and the same variable names in body.  When the function is evaluated it must be given a corresponding set of arguments, which are substituted in order.

s  =  Function[{x, y, z}, x (y^(-3)) Erf[z]]

Clear[x];
s[x, (1/y), 4]

s[x, x, x]

Of course, you do need to supply enough input to obtain a complete evaluation.

s[1]

Notice that even though the variables appear in Function without blank patterns (underscore) attached, there is no conflict between symbols used outside of the function definition and those used inside, even when the same names are involved.

x = "stuff";  y = "more stuff";
h = Function[{x, y}, x + y]

h[x, y]

Clear[a, b];
h[a, b]

The positions of the arguments matter, not their names.

Clear[x, y, z];
Function[{x, y, z}, x (y^2) (z^3)][z, x, y]

Express as pure functions:	
a) f[x_, y_] := (x^y)	b) g[x_,y_]:={Cos[x],Sin[y]}

Anonymous functions

Why bother giving a name to a specialized function that we will probably use just once and then discard?  Omission of a name avoids storing unnecessary definitions and conserves computational resources; it also relieves us of the responsibility for conceiving a memorable name that does not conflict with other symbols.  Nameless functions are obviously called anonymous.

Clear[a, b, c]

Function[x, a (x^2)  +  b x  +  c]

Function[x, a (x^2)  +  b x  +  c][q]

Function[x, a (x^2)  +  b x  +  c][Tanh[y]]

Function[x, a (x^2)  +  b x  +  c]@Tanh[y]

Skeletal notation

Since the arguments of a pure function are only local variables, why bother naming them?  If a pure function is only to be used once, why bother naming it or using the word Function?  There is a very convenient notation that can be used to define pure functions without cluttering the world with useless names.  Thus, body& is equivalent to Function[var,body] with each variable replaced by a Slot.  If only one variable is needed, it can be represented by the unnumbered slot #; but if several are needed their slots are numbered sequentially as #1, #2, etc.  MacDonald coined the term skeletal notation to describe this method of expressing functional relationships in bare form without fleshy names.

Clear[a, b, c, x]

\(\(a\ #\^2\  + \ b  #\  + \ c\  &\)\ [x]\)

\(\(\(h\)\(=\)\(a\ #\^2\  + \ b  #\  + \ c\  &\)\(\ \)\)\)

h[stuff]

g  =  (#1  +  
                2  #2) Exp[-#(3^2)] &;

g[x, y, z]

Many built-in functions are pure — we can manipulate such functions without specifying arguments.  For example, here we define f1 to be the tangent function and f2 to be its derivative.

Clear[f1, f2];
f1 = Tan;
f2 = f1'

We find that f2 is returned as a pure function given with skeletal notation.  There is no need to name the slot.  If we want a particular value, we simply fill the slot with the desired argument.

{f1[30 Degree], f2[30  Degree]}

A pure function can be applied to a list which contains the required arguments using the Apply operator @@.  The syntax f@@expr replaces the head of expr with f; thus, if the expression is a list of the form expr=List[a,b,⋯], normally printed as {a,b,⋯},  the result is simply f[a,b,⋯].

args  =  {x, y, z};
g @@ args

FullForm[args]

If the list is long enough the slots will be filled from the beginning of the list,

g @@ {x, y, z, "more stuff"}

but if it is too short the expression cannot be evaluated fully, leaving an unfilled slot.

g @@ {x, y}

Translate the following expressions into skeletal form and give examples of their use.
a) f[x_] := (x^2) Cos[x]	b) g[x_]:=Expand[x]		c) h[x_,y_,z_]:=Min[x,y,z] Max[x,y,z]

Express the logistic map μ x (1-x) using a pattern argument for the parameter μ and a slot for the independent variable x.  Then evaluate Nest[f[α],y,4] where f is your function.

Write a skeletal function which selects the (n^th) element of a list.  Test by mapping your function onto a list of the form {{a,b,c},{1,2,3},{"eeny", "meeny", "miny", "moe"}}.

Write and test a skeletal function which produces polar coordinates {r,θ} from Cartesian coordinates {x,y}.

Example: simplification using skeletal pure functions

Often one develops a symbolic result in a series of steps, with the result of each step fed to a function which performs the next transformation.  Such calculations are conveniently performed using the postfix operator (//) to add each step of the sequence as the calculation is developed.

For example, suppose we wish to create a function which produces the coefficient of (x^3) that appears in the expansion of (((a  +  b x))^13).  A simple way to accomplish this task is to supply the symbolic expansion to a skeletal function that extracts the desired coefficient.  The result of the first step, the expansion, is the argument that fills the slot in the next step reserved by #.  The presence of & indicates an anonymous pure function.

Expand[(((a + b x))^13)] // 
      Coefficient[#, x, 3] &

Generalizing this procedure, we can write a function with extracts the coefficient of (x^n) from the expansion of (((a  +  b x))^m) for arbitrary n, m.  Notice that we still use delayed assignment to prevent premature evaluation of the part of the definition using Expand.

Clear[f];
f[n_, m_] :=  
      Expand[(((a + b x))^m)] // Coefficient[#, x, n] &

f[3, 13]

Notice that the internal representation employs the standard form (body)&[arg], although I prefer to define such functions using postfix form.

?f

An important virtue of this style of coding is that it avoids creation of unwanted symbols for intermediate results.  Therefore, we often use anonymous pure functions to perform multistep transformations in a single cell without creating superfluous definitions for intermediate results.

Example: interpolating functions

Many Mathematica functions, such as Interpolation or NDSolve, return results in the form of interpolating functions.  An InterpolatingFunction is used to represent an approximate function whose values are obtained by interpolation within a table.  Thus, InterpolatingFunction[domain,table] is a pure function such that InterpolatingFunction[domain,table][x] computes the function at argument x defined within a specified domain using values contained in a table.  In standard output only the domain is printed, with the table indicated schematically by  <> .  Note that if the argument is outside the domain, a warning message is issued and an extrapolated value is returned; however, extrapolation tends to go bad very quickly outside the domain of the interpolating function and should be avoided.

The motion of a pendulum is described by the differential equation θ''[
          t]  ==  -Sin[θ[t]].  The following function evaluates that motion for a range of times assuming initial conditions θ[0] == θ0 and θ'[0] ==  ω0.  Notice that both the differential equation and the initial conditions must be formulated as equations, using == instead of =.  It is also necessary to specify the time range.  We use delayed assignment and give arguments in the form of blank patterns.

swing[θ0_, ω0_, tmin_, tmax_]  := 
      NDSolve[{θ''[
              t] ==  -Sin[θ[t]], θ[
              0] == θ0, θ'[
              0] == ω0},  θ,  {t, tmin, 
          tmax}]

A particular evaluation for specific times and initial conditions

swing[((3  π)/4), 0, 0, 10]

returns a solution in the form of a replacement rule which replaces θ by an interpolating function.  We can now produce an explicit function of t by using this replacement rule.

swing1[
        t]  =  θ[
          t]  /. swing[((3  π)/4), 0, 0, 
            10][[1]]

A phase diagram of angular velocity versus angle can then be obtained using ParametricPlot.

ParametricPlot[Evaluate[{swing1[t], D[swing1[t], t]}], {t, 0, 10}, 
        AspectRatio -> Automatic];

Notice that the phase trajectory is flattened, like an American football, whereas the phase trajectory for simple harmonic motion would be circular.  The pendulum approaches the simple harmonic oscillator as the amplitude of the swing approaches zero.

Produce a figure that displays phase trajectories for initial angles of {(Subscript[θ, 0]), (π/10), 
        9  (π/10), (π/10)}, where we use Mathematica iterator notation.  [Hint: Table can be used to form a list of pairs of functions of the form {{(Subscript[θ, 1])[t],  (Subscript[θ, 1]^′)[
            t]}, {(Subscript[θ, 2])[t],  (Subscript[θ, 2]^′)[
              t]} ⋯} , where each pair represents the phase trajectory for a different initial condition, and ParametricPlot can plot this entire family of trajectories simultaneously.]

Recursive functions

A function is described as recursive if it is defined in terms of itself.  Perhaps the simplest example is the factorial function defined by the recursion relation factorial[n]  =  n factorial[n - 1] with initial condition factorial[0] = 1.

factorial[0] =  1;
factorial[n_] :=  n factorial[n - 1]

Unfortunately, this function as written behaves poorly if given an argument which is not a nonnegative integer.  We will learn how to restrict such functions shortly, but at the present time we prefer to focus on the operation of valid recursive functions without worrying about such problems.

Table[factorial[n], {n, 0, 10}]

The simple function given above is somewhat inefficient because each time it is evaluated it must assemble an entire sequence of values.  This can be a serious inefficiency for functions based upon a complicated recursion relation.

Trace[factorial[10]]

A more efficient technique, called dynamic programming, is to use a definition which remembers its values.

factorial2[0] =  1;
factorial2[n_] := factorial2[n] = n factorial2[n - 1]

The first call for a dynamic recursion relation may actually slow the evaluation because rules must be established for the values it will remember later, but subsequent evaluations will be simpler and faster.

factorial[100] // Timing

factorial2[100] // Timing

factorial2[100] // Timing

Trace[factorial2[10]]

However, dynamic programing costs memory to store the rules.  Compare the rules associated with factorial versus factorial2.

?factorial

?factorial2

Therefore, it is not always obvious that dynamic programming is preferred — which resource needs to be conserved, time or memory?  Furthermore, there are limits to recursion by either method.

factorial[1000] // Timing

Although such limits can sometimes be increased, usually it is better to rethink the problem!

ClearAll[factorial2]

Interest is compounded monthly at rate mpr, monthly percentage rate.  Write a recursive function for the monthly balance given a unit initial investment.

Write a recursive version of the double factorial function, which is defined for positive integers as dfact[n]=n(n - 2)(Subscript[⋯n, min]) where (Subscript[n, min]) = 1 for odd or (Subscript[n, min]) = 2 for even integers.  What happens if you give this function a negative or noninteger argument?  The next section deals with restrictions on arguments.

Fibonacci numbers satisfy the recursion relation (Subscript[F, n]) =  (Subscript[F, n - 1]) +  (Subscript[F, n - 2]) with (Subscript[F, 1]) = (Subscript[F, 2]) = 1.  Write a recursive function for (Subscript[F, n]) using dynamic programming.  Note that you can make your definition for (Subscript[F, n])_.   [Do not worry about illegal arguments at this time.]

Conditions on arguments

Type checking

A function definition can be restricted to arguments of a particular type using the syntax x_type where type may be Integer, Real, Complex, List, or Symbol.

Clear[f];
f[x_Integer] :=  (x^2);
f[x_Real]  :=  (x^3);
f[x_Complex]  :=  Abs[x]^2;
f[x_List]  :=  (((Plus @@ x))^2)

This peculiar function gives different results for arguments that have the same numerical values but different types.

{f[2],  f[2. ]}

f[1 + I]

Here Plus@@x sums a list, which is then squared to produce the desired result.

f[{1, 2, 3, 4}]

However, f[a] remains unevaluated because a is a symbol of unknown type.

f[a]

We can add a definition for symbols also, if we like.

f[x_Symbol] :=  "whatever you want"

f[a]

This still doesn't give a result for f[(a^2)] because (a^2) does not match the pattern _Symbol.

f[(a^2)]

When in doubt, use the function MatchQ to determine whether an expression matches a pattern.

{MatchQ[a, _Symbol], MatchQ[(a^2), _Symbol]}

?f

Write a recursive factorial function that only responds to integer input.

Definitions with conditions

Conditions can be applied to the rules that define an explicit function.  This technique can be used with multiple definitions to produce functions whose limiting cases evaluate properly or to define discontinuous or piecewise functions.  The basic syntax uses the Condition operator (/;) to apply conditions to one or more pattern arguments 
		f[⋯, y_ /; test, z_, ⋯ ] := rhs 	
or to a definition
		f[x_] := rhs /; test
where the right-hand side is returned only if all conditions evaluate to True.  Conditions can be used with SetDelayed (:=) or ReplaceDelayed (:>) but generally should not be used with Set (=) or Replace (->) because evaluation occurs immediately, after which the condition is no longer present in the definition of the function.  Also note that all names appearing in a condition must appear in the pattern to which it is applied; otherwise, use conditional statements If, Which, or Switch as discussed in a later section.  Conditions should be applied to the smallest element containing all necessary symbols and can be applied either on the lhs or the rhs.  Thus,

Clear[f];
f[x_, y_]  :=  x  -  y  /; x > y;
f[x_, y_]  :=  y  -  x  /;  x ≤ y

successfully computes |x - y| for numerical arguments

{f[4, 2],  f[0.2, 3], f[x, y]}

but

Clear[g];
g[x_ /; x > y, y_]  :=  x  -  y ;
g[x_ /;  x ≤ y, 
        y_]  :=  y - x

{g[4, 2],  g[0.2, 3], g[x, y]}

does not because y is not present in the pattern x_.  Nor does it help to apply the condition after y_ (try it!).  On the other hand,

Clear[f];
f[x_  /;  x > 0]  :=  x;
f[x_  /;  x ≤ 0]  :=  -x

{f[-6],  f[0],  f[π]}

Clear[g];
g[x_ ]  :=  x    /;  x > 0;
g[x_ ]  :=  -x   /;  x ≤ 0

{g[-6],  g[0],  g[π]}

both compute |x| because the condition is formed properly on either side of the assignment.

When possible, it is more efficient to apply the condition to the argument rather than to the rhs because the latter method requires evaluation of each rhs for which a rule exists whereas the former method can be accomplished with fewer evaluations.

Conditions can be combined using logical operators to form a single logical condition that evaluates to either true or false.  Composite conditions should be enclosed in parentheses to ensure proper parsing.  The following example illustrates some composite conditions.  Note that the function remains unevaluated if the truth of condition cannot be determined or if a suitable definition is absent.

Clear[f];
f[x_] := "I am a positive even integer"  /; (x > 0 && 
            Mod[x, 2] == 0)
f[x_] := "I am a negative odd integer"  /; (x < 0 && 
            Mod[x, 2] == 1)
f[x_] := "I am not an integer"  /; (Mod[x, 2] ≠ 0 && 
            Mod[x, 2] ≠ 1)

f /@ {q, -2, 8, -7, π}

Example: step function

The unit step function is defined to be +1 for positive arguments or 0 for negative arguments.  By convention, we choose the intermediate value of (1/2) for the origin.  This function can be created using multiple definitions that test the argument.

Clear[step];
step[x_  /; x > 0]  := 1;
step[x_  /; x == 0]  := (1/2);
step[x_  /;  x < 0]  :=  0

{step[-(1/2)],  step[0],  step[12],  step[x],  
      step[1 - x] /. x -> 0.4}

All of these cases give the expected output except step[x], which remains unevaluated because the condition is ambiguous, being neither definitely True nor definitely False.

Note that the condition can be applied to either the left or the right sides of the assignment.

Clear[step];
step[x_ ]  :=  1    /; x > 0;
step[x_ ]  :=  (1/2)  /; x == 0;
step[x_ ]  :=  0    /;  x < 0

{step[-(1/2)],  step[0],  step[12],  step[x],  
      step[1 - x] /. x -> 0.4}

Example: spherical Bessel function at origin

Recall our previous definition for the regular spherical Bessel function.

Clear[sbessj];
sbessj[n_, z_]  :=  (Sqrt[(π/2)]) 
          (z^(-1/2)) BesselJ[n + (1/2), z];

Unfortunately, this simple definition cannot handle the origin.

{sbessj[0, 0], sbessj[1, 0], sbessj[n, 0]}

Since we should know the limits for these functions as x -> 0, or can deduce them using Limit, where we specify that the limit is to be evaluated for decreasing x.

Limit[sbessj[0, x], x -> 0, Direction -> -1]

Limit[sbessj[1, x], x -> 0, Direction -> -1]

We can include these properties in the definition of the function itself by applying conditions to the arguments and including several definitions covering various special cases.

Clear[sbessj];
sbessj[n_, 
          z_  /; z ≠ 
              0]  :=  (Sqrt[(π/2)]) (z^(-1/2)) 
          BesselJ[n + (1/2), z];
sbessj[0, 0]  :=  1;
sbessj[n_ /; n > 0, 0 ]  :=  0

Plot[
        Evaluate[Table[sbessj[n, x], {n, 0, 3}]], {x, 0, 10}];

This is certainly an improvement, but it is still not as complete as we would like; for example, it does not provide explicit derivatives.

D[sbessj[2, z], z]

An improved version will be developed later.

Exercises

Write and plot a function whose values are 1 - |x| for x ≤ 1 or 0 for x ≥ 1 without using Abs.

The double factorial function is defined for positive integers as dfact[n]=n(n - 2)(Subscript[⋯n, min]) where (Subscript[n, min]) = 1 for odd or (Subscript[n, min]) = 2 for even integers.  Write a recursive version of this function using conditions that restrict the argument to positive integers.  Be sure to test the function.  [Hint: you can both specify the argument type and apply a condition to the argument.]

Write a function which returns +1 for positive even integers, -1 for positive odd integers, or 0 for anything else using conditions (do not use If).  Be sure to test your function.  [Hint: the functions OddQ and EvenQ are helpful.]

Suppose that the earth is modelled by a sphere of uniform density with mass M and radius R.  Then, according to Newton's law of gravity, the gravitational force per unit mass, g(r), that acts on a body would be G M/(r^2) for r ≥ R or G M r/(R^3) for r ≤ R, where G is the gravitational constant.  It is simplest to employ dimensionless units in which {G->1,M->1,R->1}; for any particular case, the scale can be adjusted later by multiplying g(r) by the surface gravity G M/(R^2).  Write a function for g(r) and plot it in dimensionless form.

Pattern tests

Some conditions are best applied using predicates, which are built-in functions which ask questions that can be answered as either True or False.  Most of these predicates have function names which end in Q for question.  A list of the available predicates is easily obtained, plus a few interlopers such as LegendreQ (the Legendre function of the second kind is usually designated (Subscript[Q, n])).

?*Q

Among these some of the most important are:

?IntegerQ

?EvenQ

?OddQ

?NumberQ

?NumericQ

?VectorQ

?MatrixQ

Unfortunately, there are also several exceptions to the naming convention for predicates.

?Positive

?Negative

Predicates are used to test arguments using the syntax x_?predicate.  Note that in this context the question mark is interpreted as a PatternTest, which should be distinguished from the information escape character, with the same appearance, used to obtain interactive help.  The information escape cannot be embedded in composite expressions, whereas the pattern test can be used in expressions liberally.

For example, the following function

Clear[f];
f[x_?EvenQ] :=  (x^2);
f[x_?OddQ] :=  -(x^3)

is defined for integers

Table[f[n], {n, -5, 5}]

but is undefined for real numbers.

f[1.5]

A pattern test can be used with a type-specific pattern.

Clear[f];
f[x_Integer?
          Positive] := "I am a positive integer"
f[x_Integer?
          Negative] := "I am a negative integer"
f[x_Real?
          Positive] := "I am a positive real number"
f[x_Real?
          Negative] := "I am a negative real number"
f[x_Symbol] := "I don't know what I am!"

f /@ {q, 14, -27, π, N[π]}

More generally, x_?test is used to test a pattern using a pure function that returns either True or False when presented with an argument.  Several tests can be combined using logical operators.  For example, the following function is defined only for arguments which are positive prime numbers.

Clear[f, x];
f[x_?(PrimeQ[#]  &&  Positive[#] &)]  :=  (1/x)

f /@ {-4, -3, x, 16, 17, "stuff"}

The construction x_?(⋯) indicates that the pattern x_ is tested using the function contained within the parentheses.  The function is composed of two logical functions connected by a logical operator, in this case And or &&, to make a single logical function.  Our composite logical function receives a single argument in its slot, #, and is identified as a pure function by the & operator.  The parentheses are needed with skeletal notation to identify the beginning and end of the pure function.  The same effect could have been achieved by defining the test using Function[var,body] instead.

Rewrite the function above using conditionals instead of a pattern test.

Write a recursive factorial function which only responds to arguments which are nonnegative integers.  Why are these restrictions necessary?

Write a recursive double factorial function for nonnegative integers using a pattern test.

Example: a function which tests orthogonality

For example, suppose that we need a function that tests for the orthogonality of a matrix.  A matrix is orthogonal if when multiplied by its transpose an identity matrix is obtained.  Since orthogonality is meaningless for other objects, we use the predicate MatrixQ to restrict the definition to matrices.  We verify that the matrix is square using Dimensions and compare the product with its transpose to the identity matrix of the same length.

\(\(\(Clear[testOrthogonality];\)\(
\)
    \)\)
\(\(testOrthogonality[matrix_?MatrixQ] := 
        MatchQ[Dimensions[matrix], {x_, x_}]\  && \ 
          matrix\  . \ Transpose[matrix] == 
            IdentityMatrix[Length[matrix]];\)
\)
testOrthogonality[x_?NumericQ]  :=  (x^2) == 1

In the following cases our function successfully identifies a manifestly orthogonal matrix, rejects a matrix that is not square, and does not venture a guess about a symbol of uncertain character.

testOrthogonality[{{(1/(Sqrt[2])), (1/(Sqrt[2]))}, {-(1/(Sqrt[2])), 
          (1/(Sqrt[2]))}}]

testOrthogonality[{{(1/(Sqrt[2])), (1/(Sqrt[2])), 0}, {-(1/(Sqrt[2])), (1/(Sqrt[2])), 
          1}}]

testOrthogonality[x]

We also included a definition to handle expressions which evaluate to numbers, using NumericQ rather than NumberQ.  (Convince yourself that NumberQ is not up to the task!)

testOrthogonality /@ {1, 2, -1,  I,  (1/(Sqrt[2])) (1 + I),  
        Cos[0], Cos[0.25]}

A slightly more difficult case:

testOrthogonality@(({{(1/3), ((2 (Sqrt[2]))/3), 0}, {(2/3), -(1/(3 (Sqrt[2]))), (1/(Sqrt[2]))}, {-(2/3), (1/(3 (Sqrt[2]))), (1/(Sqrt[2]))}}))

Unfortunately, this version does not work well with matrices that contain approximate real numbers — we can use Equal with exact numbers but for approximate numbers the most we can expect is that the matrix times its transpose produces an identity matrix within some small but finite tolerance.  Check that if we use N to convert the above matrix to real numbers that the simple orthogonality test fails.  Then write an improved version that handles approximate numbers correctly.  [Hint: test the discrepancy with the largest absolute value.]

Arguments matching a specialized pattern

Suppose that we need a function that acts upon data only of a specialized pattern, such as a list of pairs.  The basic syntax is x:pattern where x is the name used on the rhs of the definition to refer to the pattern and pattern is one of the many types of patterns recognized by Mathematica.  For example, the function

Clear[f];
f[pairs : {{_?NumericQ, _} .. }]  :=  
      Sort[pairs,  #2[[1]] ≥ #1[[1]] &]

sorts a list of pairs in order of increasing first elements.  Sort[list,test] sorts a list using a logical function test to test the ordering of pairs of neighboring elements.  Our test employs an anonymous function with skeletal notation.  We used the part specification [[1]] to extract the first elements of each pair in the list and the logical operator ≥ to compare successive pairs that fill the two slots #1 and #2 in that anonymous function.  Sort then arranges the elements of the list to ensure that the test applied to neighboring elements of the output list is True for all neighbors.  The pattern contains two sets of curly brackets to describe a list of lists, with .. indicating one or more repetitions of the inner pattern.  The inner pattern describes a pair whose first element in numerical (according to the pattern test) while the second element is arbitrary (a blank pattern).

Thus, this function successfully orders a list in which all first elements of each pair are numerical expressions

Clear[x, y];
f[{{1, 2}, {-1, x}, {-4, y}, {12, 3}, {Cos[0.2], 
          Sin[0.2]}}]

but does nothing when one or more of the first elements of a pair cannot be tested properly

Clear[a, b];
f[{{1, 2}, {-1, x}, {a, b}, {-4, y}, {12, 3}}]

because we used a pattern {{_?NumericQ,_}..} that required each first element to be numerical.  Of course, we could add other definitions to handle other patterns as needed.

Mathematica provides a very flexible language for describing an almost limitless variety of patterns.  For more information about using patterns, refer to expressions.nb.  Click here for a table of pattern templates.

Write a function which arranges a list of pairs in order of decreasing second elements.

Write a function which arranges a list of sublists in order of decreasing second elements, requiring only that the sublists have two or more elements with the second being numeric.

Overlapping definitions

Earlier we warned about possible ambiguities that might be encountered if a function has several definitions which might overlap or conflict with each other.  Under those circumstances, which rules take precedence?  The general rule is to look for matching patterns beginning with the most specific definition and proceeding to more general definitions if needed.  If two definitions are judged equally specific, the first is used.  However, as noted by Tam, Mathematica's notions of generality are often inscrutable.

Consider, for example, the following functions discussed by Tam.  These functions appear to differ only in the order their definitions are given.

Clear[f];
f[x_Integer]  :=  (x^2);
f[x_Integer?Positive]  :=  1/x;
f[x_Integer  /;  x > 50]  :=  -x

Clear[g];
g[x_Integer?Positive]  :=  1/x;
g[x_Integer  /;  x > 50]  :=  -x;
g[x_Integer]  :=  (x^2)

If the argument has the value 60, is a "positive integer" more specific than "integer > 50"?  Surely both of those tests are more specific than "any integer", yet f uses the definition for any integer while g uses the definition for positive integer without questioning its size.

{f[60], g[60]}

Apparently Mathematica did not distinguish specificity but used only order.  However, we do not consider this observation to be robust and strongly recommend that function definitions should only involve compatible or nonoverlapping conditions.

Protection of local variables

Need for scoping

Protection of local variables is often crucial to the proper operation of a function.  The following example is discussed by Bahder.  Consider a function defined in terms of an integral.

Clear[h, x, y];
h[x_]  :=  Integrate[Sin[x + y], {y, 0, π}]

This definition works under many circumstances,

h[x]

h[2]

but fails when the argument conflicts with the dummy variable of integration

h[y]

As revealed by using Trace to follow the evaluation of the expression h[y],

Trace[ h[y] ]

the problem arises when y replaces the pattern x_ too early in the evaluation sequence and is confused with the integration variable.  One solution is to use Module to protect the dummy variable.  The Module command provides a method for encapsulating a sequence of commands together with a declaration of local variables.  The basic syntax is Module[{vars},body] where {vars} is a list of local variables and body is the body of the function consisting of a sequence of commands separated by semicolons.  The final result is the output of the last statement, provided that it is not terminated by a semicolon.

Clear[h, x, y];
h[x_]  :=  
      Module[{y},  Integrate[Sin[x + y], {y, 0, π}] ]

h[y]

This method works because the local variable is given a unique tag each time that Module is evaluated so that its local name will not conflict with any symbols that exist outside of that Module.  The basic concept of restricting the definition of symbols to limited parts of a code is known as scoping.  The scope of a variable is the range of code wherein it is defined.  Symbols with the same names can be used for different purposes if their scopes do not overlap.

Define a function f[n_]:=(∑_(k = 0)^n) (Subscript[a, k]).  Show that a scoping command, such as Module, is needed to protect the summation index.  Evaluate and interpret: a) f[n], b) f[k], c) f[10]-f[3].

Method for scoping

Every time Mathematica evaluates a module, it creates new symbols by attaching a unique tag to the symbol for each local variable, such that the local variable var is represented by the symbol var$nnn where nnn is the current value of the system variable $ModuleNumber.  $ModuleNumber is incremented every time that a Module is evaluated to ensure that the local variables are unique and will not conflict with any global symbols.

What is the current value of $ModuleNumber?  Evaluate a Module statement from the previous section and check the new value of $ModuleNumber.

We can peek at the operation of Module using Trace.  Upon evaluation of the following statement

Trace[ h[y] ]

you should see an output in which the {y} used to declare y as a local variable is replaced by {y$} and later the evaluation of the integral features both y and y$nnn where nnn is the current value of $ModuleNumber.  Thus, there is no confusion between the global symbol y used as an argument and the local symbol y$nnn used as a dummy variable of integration.

Sometimes one sees the symbols for local variables in an error message.  Alternatively, you can also access these symbols if you know the appropriate tag number.  For example, the following module resets the seed employed by the random number generator using the number of seconds between 1900 and the present, rounded to an integer.  Using Trace we can learn the symbol assigned to the local variable and then determine its value even though Module attempts to hide local variables.

Evaluate the following statement and then determine the seed used.

Trace[ 
      Module[{seconds}, seconds = Round[ AbsoluteTime[] ]; 
        SeedRandom[seconds]; Random[]] ]

Modular programming

Well-written programs tend to be modular, with each module receiving a limited amount of input, performing a specific task, and returning clearly defined results.  The internal programming of each module should be complete enough to cover all special cases and contingencies, thereby producing accurate results for valid input or complaining to the user about invalid input.  The modules should also contain enough documentation to be readable by an informed user without requiring human compilation of computer code.  Mathematica provides several constructions which serve these purposes.

Scoping commands

Several methods are available and the most important differences between these methods concern the treatment of local variables.  Scoping commands may be nested, with scoping rules applied from the inside out.

Module

The Module command provides a method for encapsulating a sequence of commands together with declaration of local variables.  The basic syntax is Module[{vars},body] where {vars} is a list of local variables and body is the body of the function consisting of a sequence of commands separated by semicolons.  The final result is the output of the last statement, provided that it is not terminated by a semicolon.  Alternatively, Return[expr] can be used to return an expression.

Local variables can be given initial values using the syntax Module[{x = (Subscript[x, 0]),⋯},body] where any variable in the local list can be initialized or not.  The initial values can be numerical, symbolic, or the result of evaluating a function.  However, for best results the initial values should be based only upon constants or explicit arguments — avoid using global symbols in the initialization section of Module.

The following example takes an electrostatic potential ϕ and displays its electric field and equipotential lines in two dimensions.  The electric field is evaluated by forming the gradient of the potential in the initialization part of the module.  Two plots are constructed as local variables but are not displayed until the final statement, where Show combines them into a single figure.

Needs["Graphics`PlotField`"];

plot[ϕ_] := 
      Module[{Efield = -{D[ϕ, x], D[ϕ, y]}, plot1, plot2}, 
        plot1 = ContourPlot[ϕ, {x, -2, 2}, {y, -2, 2}, 
            ContourShading -> False, DisplayFunction -> Identity]; 
        plot2 = PlotVectorField[Efield, {x, -2, 2}, {y, -2, 2}, 
            MaxArrowLength -> 1, ScaleFactor -> None, 
            DisplayFunction -> Identity]; 
        Show[plot1, plot2, 
          DisplayFunction -> $DisplayFunction];]

We test the function using a simple dipole field.  To avoid confusing clutter, the strong fields near the point charges are not displayed.  Fancier versions of this function might include additional variables for the plotting ranges and options for the plotting commands that could be used to improve the display, but such embellishments are beyond the scope of the present discussion.

ϕ = 
        (1/(Sqrt[(x^2) + (((y - 0.5))^2)])) - 
          (1/(Sqrt[(x^2) + (((y + 0.5))^2)]));
plot[ϕ];

Notice that the local variable Efield is not defined outside the module.

Efield

Write a function f[x,y] using Module that returns (∑_(k = 0)^n) ((x^k)/(k!)) where n is the largest integer less than or equal to y.  Of course, your function should only respond to positive y.  Verify that the function works properly and that the local variables are protected properly.

With

With is another scoping command designed to assign values to local constants.  The basic syntax is With[{x=x0,y=y0,⋯},expr] where the expression expr is evaluated using local constants x, y, ... assigned local values (Subscript[x, 0]), (Subscript[y, 0]), ... .  The principal distinction between With and Module is that With employs local constants that are evaluated only once, whereas Module employs local variables whose values may change many times as the body of the module is evaluated.  Thus, code written with With is sometimes more transparent to the reader than equivalent code based upon Module.  The local constants do not escape the scope of With and may have the same names as external global symbols.  The expression evaluated by With may itself be a Module.

a  =  "stuff";  b = "more stuff";  
    c = "yet more stuff";

ff[x_]  :=  
      With[{a = 2, b = 3, c = 4},  a (x^2) + b x  +  c]

ff[x]

{a, b, c}

Note that you cannot reassign a local constant within With.

gg[x_]  :=  
      With[{a = 2, b = 3, c = 4},  a = 12; b = -12; c = "junk"; 
        a (x^2) + b x  +  c]

gg[y]

Note that With protects local constants by appending $ to the corresponding symbols, a method similar to that used by Module.  Verify by using Trace[ff[x]].

Plot the function (((Sin[x  +  y]/(x + y)))^2) using With to assign the value y = π.  Verify that y can be assigned a value that is not changed by your assignment within With.

Block

Another type of scoping command is Block.  The principal difference between Block and Module is that Block treats the values assigned to symbols as local but treats their names as global, whereas Module treats the names of local variables as local.  The basic syntax is either Block[{x,y,...}, expr] or Block[{x = (Subscript[x, 0]),y = (Subscript[y, 0]),...},expr] where {x,y,...} are assigned local values during the evaluation of the expression expr, and may change during that evaluation, but their names exist outside the scope of Block.  If initial values {(Subscript[x, 0]),(Subscript[y, 0]),...} are supplied, these are used as local starting values.  Note that commands, such as Table and Do, implicitly use Block for handling iterators.

Clear[f, x]

f[x_] := 
      Block[{a = 2, b = 3, c = 4},  a (x^2) + b x  +  c]

{a, b, c}

f[x]

{a, b, c}

Notice that Block has not changed the values of the global variables.  Similarly, Table can use an iterator whose name has been assigned a global value without interference with that value.

Table[(a^2), {a, 0, 4}]

a

However, Block does not always prevent confusion between local and global symbols.  Compare the following versions of a function which adds reciprocal powers of its argument.

f[x_] := Module[{i}, Sum[(1/(x^i)), {i, 0, 5}]]

g[x_] := Block[{i}, Sum[(1/(x^i)), {i, 0, 5}]]

{f[2], f[x], f[i]}

{g[2], g[x], g[i]}

The version using Block fails because the local variable conflicts with the argument of the function.  Normally one should use Module instead of Block, but there are circumstances where Block can be useful.  For example, if one needs to temporarily modify a system parameter within a small section of code, Block is helpful but Module does not work because it appends extra stuff to the name of the local variables.  Bahder suggests the following example.

Clear[f];
f[0] = 1;
f[n_] := n f[n - 1];

blockFactorial[n_] := 
      Block[{$RecursionLimit = n + 3},  f[n]]

moduleFactorial[n_] := 
      Module[{$RecursionLimit = n + 3},  f[n]]

blockFactorial[280]

moduleFactorial[280]

Nevertheless, Module or With is preferred under most circumstances.

Conditional functions

If, Which, and Switch

A conditional function is a command whose actions depend upon conditions.  We have already explored how conditions can be applied to definitions or to rules using the condition operator (/;).  More sophisticated conditional branching can be implemented using the If, Which, or Switch functions.

The basic syntax for the If function is If[condition,t,f], where the expression t is returned if the condition is true or the expression f is returned if the condition is false.  However, the entire command is returned unevaluated if the condition is indefinite.

Clear[step];
step[x_]  :=  If[x > 0, 1, 0]

{step[-1], step[0], step[2], step[x],  
      step[x - 1] /. x -> 2}

Although this version of the step function is more compact than the earlier version using three conditional definitions, the value for step[0] is not what we prefer; a better version based upon Which is given below.  Also note that the result for an argument which yields an ambiguous condition remains unevaluated.  An alternative form which always returns something is If[condition,t,f,u], where the expression u is return if the status of the condition is uncertain.

Clear[f];
f[x_]  :=  If[x > 0, 1, 0, 0]

{f[-1], f[0], f[2], f[x],  f[x - 1] /. x -> 2}

We chose a different name because this variant isn't really a well-defined step function — the interpretation of unknown variables remains problematical.

A more flexible conditional function is the Which command with basic syntax Which[test1,expr1,test2,expr2,⋯].  Each test is evaluated in turn and the expression corresponding to the first true case is returned.  With this command we can construct a step function that returns the desired value at the transition point but does not produce unreliable results for uncertain input.  A default value can be given simply by choosing True as the final test in the sequence.  However, if a test is encountered that cannot be evaluated as either True or False, the statement is returned unevaluated.

Clear[step];
step[x_]  :=  
      Which[x > 0, 1, x == 0,  (1/2),  x < 0,  0]

{step[-1], step[0], step[2], step[x],  
      step[x - 1] /. x -> 2}

Finally, an even more sophisticated conditional is Switch which combines pattern testing with conditional evaluation.  The basic syntax is Switch[expr,form1,value1,form2,value2,⋯].  The input expression expr is compared with each pattern form until a match is found and the corresponding value value is returned.  Note that a default value can be given simply by ending with a blank pattern.  The following example comes from the Mathematica Book.

\(switchfunction[x_] := 
      Switch[x, 
_Integer, FactorInteger[x], 
_Rational, N[x], 
_Real, 
        Rationalize[x], 
_Complex, Im[x], 
_Symbol, 
        1 + x\^2, 
_, "The Head is " <> 
          ToString[Head[x]] <> "."]\)

switchfunction /@ {12, (1/3), 9.1, 
        2 + 7 I, π, "xxx yyy"}

Write a function whose values are 1 - |x| for x ≤ 1 or 0 for x ≥ 1.  Make several variations using conditioned arguments, If, or Which.

Use If or Which to write a function that gives the gravitational acceleration for a uniform sphere.  (Refer to the earlier version of this problem for details.)

Finite-difference methods for differential equations represent the solution on a uniform grid (Subscript[x, i + 1]) = (Subscript[x, i]) + dx with spacing dx as a vector f[(Subscript[x, i])] for i = 1, n and approximate the second derivative by 
	\(\(f\^″\)[
              x\_i] = \(f[x\_\(i + 1\)]\  - \ 2  f[x\_i] + f[x\_\(i - 1\)]\)\/dx\^2\) 
Thus, the second derivative can be constructed by matrix multiplication, (f^″) = A . f, where A is an operator acting on the solution vector.  Use Table and Switch to write a function A[n_,dx_] that constructs the desired matrix.

Use Switch to write a function that in a single statement defines the area of simple geometric figures, including triangle, rectangle, parallelogram, trapezoid, etc.  Compare with the earlier version using multiple definitions.

Boole and Piecewise: new in Version 5.1

Two new conditional functions were introduced in Version 5.1.  Boole[expr] returns 1 if the logical expr evaluates to 1 if true or 0 if false.  This simple function is equivalent to If[expr,1,0] but is handy nonetheless.  A more sophisticated conditional function is Piecewise[{{val1,cond1},{val2,cond2}⋯},val] which returns the first value for which the corresponding condition is true.  If none of the conditions evaluate to True, the final value is returned; if omitted, the default final value is 0.  If any of the conditions are indeterminate, the result is only partially evaluated.  Piecewise can be nested and can be used in integrals, differential equations, and in most applications suitable for conventional functions.

For example, we can now define spherical Bessel functions by providing an explicit formula for n = 0 that interprets the removable singularity at the origin.

Clear[sbessj];
sbessj[0, z_] = 
        Piecewise[{{1, z == 0}}, 
          Sin[z]/z];
sbessj[n_, z_] = (Sqrt[(π/2)]) (z^(-1/2)) 
          BesselJ[n + (1/2), z];

This version behaves well numerically.

Plot[
        Evaluate[Table[sbessj[n, x], {n, 0, 3}]], {x, -10, 
          10}];

Plot[
        Evaluate[Table[D[sbessj[n, x], x], {n, 0, 3}]], {x, -10, 
          10}];

It also has the advantage that symbolic derivatives are evaluated properly.  Notice that Piecewise takes care of the limit z -> 0 for n = 0 correctly.

Table[Simplify[D[sbessj[n, z], z]], {n, 0, 2}] // TableForm

The gravitational field of a uniform sphere with mass M and radius R is
	r < 
          R  ⟹  g  =  ((G M r)/(R^2))
	r > 
          R  ⟹  g =  ((G M)/(r^2))
where G is the gravitational constant.  The corresponding potential is then
	ϕ  =  -(∫_∞^r) g[
                s] ds
Use Piecewise to define the field and then evaluate the potential.   It may help to use Assuming[{R>0,r>0}, ⋯ ] to inform Integrate that it may assume that radii are positive real numbers.  Then plot ϕ[r]/ϕ[R] versus r/R.

Loop commands

Basic syntax

The Do command can be used to execute a set of instructions a fixed number of times, whereas the While or For commands execute a set of instructions while a condition remains true.

The basic syntax for Do is Do[body,iterator] where body is a sequence of commands separated by semicolons (except for the last) and iterator is a standard Mathematica iteration code.  If several iterators are included, separated by commas, nested loops are evaluated from the innermost iterator to the outermost.  Although one can often write functions more efficiently using Nest or related functions, for numerical calculations using intermediate results it is often useful to embed a Do loop within a module.

The basic syntax for While is While[condition,body] where condition is tested while body is iterated until the condition becomes False.  You must make sure that the condition will eventually fail, perhaps by including an iteration counter in body and a test on that counter in condition.

Finally, the basic syntax for For is For[start,test,incr,body] where start is evaluated first and usually initializes a counter (and perhaps other variables), then test and body are evaluated until test fails.  The iteration variable is incremented by incr.  A trivial example is given below.  Notice that 1) start is evaluated first,  2) since the test evaluates to True, 3) body is evaluated next, then 4) incr changes something and 5) test is evaluated again to determine whether body should be reevaluated.

For[(x = 20; y = 1), x > 0, 
      x = x - 2, (y = y*x; Print[y])]

An advantage of Do is that by prescribing the number of iterations in advance one can be confident that the loop will eventually terminate.  On the other hand, loops which terminate when a condition is satisfied are more flexible.  For example, when summing a series one can continue until convergence is obtained.  However, the user must make sure that the loop does eventually terminate.  Therefore, one usually includes a maximum number of iterations as part of the termination condition.

Example: summation using Do, While, or For

For example, below we use Do, While, and For to evaluate an approximation to the following sum with symbolic result

exactsum = Sum[(n^(-2)), {n, 1, ∞}]

and numerical result

N[exactsum]

First, we use a Do loop to sum 100 terms.  Notice that we initialize the sum in the declaration of local variables for Module and return a numerical result as the final statement of the module.  The Do loop by itself accumulates the sum, but does not print the final result.  By placing the loop within a module, we both protect the local variable and ensure that the final result is displayed.  It is good practice to always encapsulate loop commands within a scoping command.

Module[{sum = 0}, Do[sum = sum + (1/(x^2)), {x, 1, 100}]; 
      N[sum]]

Alternatively, While expresses that range in terms of a condition upon the loop counter.

Module[{sum = 0, x = 1}, 
      While[x ≤ 100, sum = sum + (1/(x^2)); x = x + 1]; 
      N[sum]]

Finally, For includes the initialization but is otherwise similar to While (although the components differ in order).

For[{sum = 0, x = 1}, x ≤ 100, x = x + 1, 
      sum = sum + (1/(x^2))]; N[sum]

a) Write and test factorial functions using Do, For, and While loops.  
b) Evaluate 257! using these functions.  Note that the recursion limit encountered for the recursive factorial function does not affect loops.  
c) Compare the speeds of these three versions of the function.

Use a For loop to determine the minimum nmax needed to reduce the difference between (∑_(n = 1)^nmax) (n^(-2)) and the limiting value ((π^2)/6)to less than 0.001.  Check using Sum.  [Important: the sum approaches the limit from below and the counter is incremented before leaving the loop.]

Fancy incrementation operators

Mathematica provides many of the incrementation operators that appeal to C programmers, but perhaps to few others.  For example, we can increment x within a For loop either using the assignment statement x=x+1 or the incrementation operator ++x. Although one never needs these somewhat cryptic constructions, they are listed below for completeness.

({{Incrementation Operators}, {({{incrementation, new value, return value}, {x++, x -> x + 1, old value}, {++x, x -> x + 1, new value}, {x--, x -> x - 1, old value}, {--x, x -> x + 1, new value}, {x += dx, x -> x + dx, new value}, {x -= dx, x -> x - dx, new value}, {x *= c, x -> x*c, new value}, {x /= c, x -> x/c, new value}})}})

Note that in some cases the old value of the loop counter is printed even though the counter is updated; in other words, in some cases the counter is printed before updating.  Thus,

x = 1;  x++

prints the value before incrementing, and then stores the new value

x

whereas

x = 1;  ++x

x

increments before printing.  These two forms usually behave the same way within loops that do not print their counters, but one should be aware of the difference between them.

As an example, the For loop at the end of the preceding section is rewritten below using incrementation operators in place of assignment statements.  The results are the same, of course, but we have saved a few characters at the expense of readability.

For[{sum = 0, x = 1}, x ≤ 100, ++x, sum += (1/(x^2))]; 
    N[sum]

Documenting functions

Several types of documentation are available.  The most obvious are the text cells that we have been using to present this course!  This type of documentation describes the notebook and its contents in general, but not the specifics of a function.  A function or module should be contained within a single input cell.  Comments which explain the steps within that cell can be included using the special notation (* comment text *), where everything between (* and *) is treated as a comment and is not evaluated.

Another important type of documentation is a usage message.  When we type

?BesselJ

we receive a usage message that provides a brief explanation of the function and its arguments.  This is especially important when the function definitions are not immediately accessible to the user; perhaps they are buried in some package.  In any case, it is good practice to provide usage messages for your functions if only to remind yourself what they do a few weeks later.  It is also advisable to provide usage messages for any options that you give your functions.  The basic syntax for a usage message is f::usage = "text".  Recall, for example, our version of the spherical Bessel function.

\(\(\(Clear[sbessj];\)\(
\)
    \)\)
sbessj::usage  =  "sbessj[n,z] computes the spherical Bessel \
function of the first kind (based on BesselJ) and is regular at the origin \
for integer order n.";
\(\(\(sbessy::usage\  = \ "sbessy[n,z] computes the spherical Bessel function of the second kind (based on BesselY) and is irregular at the origin for integer order n.";\)\(
\)
    \)\)
sbessj[0, z_] = 
        Piecewise[{{1, z == 0}}, 
          Sin[z]/z];
sbessj[n_, z_] = (Sqrt[(π/2)]) (z^(-1/2)) 
          BesselJ[n + (1/2), z];
sbessy[n_, 
          z_]  :=  (Sqrt[(π/2)]) (z^(-1/2)) BesselY[n + (1/2), 
            z];

?sbess*

Usage messages are a special case of a more general messaging function that can be used to report errors or warnings.  To create a message associated with a symbol, assign a string to symbol::tag where tag identifies the particular message.  A message is printed when the function Message[symbol::tag,parameters] is evaluated.  A message string can include one or more parameter slots, numbered `1`, `2`, etc., which are filled sequentially.  Note the back-quotes in the parameter numbering.

Our previous versions of factorial functions remained unevaluated if given improper input.  In the following version we include messages that inform the user of the nature of improper input.

factorial::usage = "factorial[n] returns n! for nonnegative \
integers and error messages otherwise";
factorial::noninteger = "`1` is not an integer";
\(\(\(factorial::negative = "negative argument `1` not allowed";\)\(
\)
    \)			\)
\(factorial[n_] := \ 
      Which[
	Not[IntegerQ[n]], \ Message[factorial::noninteger, n], 
	
        n < 0, \ Message[factorial::negative, n], 
	
        True, \ \(n!\)\ ]\)

?factorial

factorial[-2]

factorial[π]

factorial[10]

Messages are usually placed in packages that contain well-developed functions which are intended to be used in many applications.  It is often important to ensure that such functions do not return unpredictable results for improper input and it is then very helpful to the user, who might even be yourself, to provide diagnostic information.

Compiled functions

Unlike traditional compiled languages which translate instructions directly into machine code that performs simple operations, Mathematica is an interpreted language which evaluates instructions by applying rules to symbols that match patterns.  Consequently, Mathematica code is sometimes rather slow compared with, say, Fortran or C code for the same function.  Fortunately, there is a mechanism in Mathematica for recovering some of the economy of compiled languages by using Compile to produce machine code that operates upon numbers rather than symbols.  Thus, Compile[{(Subscript[x, 1]), (Subscript[x, 2]), ⋯},body] produces a pure function of type CompiledFunction[{(Subscript[x, 1]), (Subscript[x, 2]), ⋯},-CompiledCode- ] assuming that all arguments are real numbers.

Clear[f, x];
f =  Compile[{x},  (x^2) Exp[-(x^2)] ]

Timing[Table[f[x], {x, 0. , 10. , 0.0001}];]

Clear[g, x];
g[x_]  :=  (x^2) Exp[-(x^2)];
Timing[Table[g[x], {x, 0. , 10. , 0.0001}];]

Thus, for this simple function compilation improves the speed by about a factor of 3.  The ratio depends upon the function and the machine, but compilation is often quite helpful for complicated functions that are expected to be used only with numeric data of particular types.

When the arguments are not real numbers, it is necessary to specify their types explicitly using the syntax Compile[{{(Subscript[x, 1]), (Subscript[t, 1])},  {(Subscript[x, 2]), (Subscript[t, 2])}, ⋯},body] where (Subscript[t, i]) indicates the type and may be _Real,  _Integer,  _Complex, or True|False.  Refer to Compile for other variations that apply to more complicated arguments.

Example: Mandelbrot set

An elementary example of procedural programming is provided by the classic Mandelbrot set, which was often the subject of the first program written by a proud new PC owner back in the early 1980's.  Here we develop such a program using many of the concepts presented in this notebook, adding refinements as we go.

The algorithm consists of the iteration

z → (z^2) + c

where c = cx + I cy is a complex number with real part cx and imaginary part cy and where the starting value is (Subscript[z, 0]) = 0.  Either z will find itself trapped near the origin or, after some number of iterations, it will escape to infinity.  The Mandelbrot set consists of all points c for which the iteration beginning with (Subscript[z, 0]) = 0 produces a sequence that is trapped.  Pretty pictures are created by applying a coloring or altitude algorithm to display the rate at which points outside the Mandelbrot set flee the scene.

Therefore, we require a function of (cx,cy).  Since it will be necessary to terminate the iteration after a finite number of steps for points which either escape very slowly or simply orbit the origin, it is convenient to specify the maximum number of iterations as a third argument.  The value returned by the function will be either the number of iterations needed to obtain |z| > 2, which one can show is threshold for certain escape, or the maximum number of iterations permitted.  We use Module to isolate the local variables of no external importance and use While to perform the iteration while testing the magnitude of z.  We use Return to clearly identify the result of the function.  Those familiar with any standard procedural programing language, such as Fortran or C or Pascal, should immediately recognize this algorithm.

\(\(\(\(Clear[mandelbrotRate];\)

    mandelbrotRate[cx_, cy_, lim_]\  := \ 
      Module[{c, z, ct = 0}, c = cx\  + \ I\ cy; z\  = \ c; 
        While[\((Abs[z] < 2.0)\)\  && \ \((ct ≤ lim)\), 
          ct = ct + 1; z\  = \ z*z + \ c;]; 
        Return[ct]\ ]\)\(\ \)\)\)

The escape rate can now be displayed in the form of a DensityPlot using a coloring function that assigns each point a hue.  Below we evaluate the hue using the logarithm of the escape rate.

DensityPlot[
        Log[1 + mandelbrotRate[x, y, 50]], {x, -2, 0.6}, {y, -1.2, 
          1.2}, PlotPoints -> 100, Mesh -> False, 
        AspectRatio -> Automatic, 
        ColorFunction -> Hue];

The figure above has relatively coarse resolution.  If we start playing with the resolution and attempt to zoom in on particular regions or to alter the color scheme or plotting format, we quickly find that the calculations are maddeningly slow.  We can do a little better by using Compile to optimize the calculation.  However, because Compile assumes that the arguments are real numbers, we need to specify that lim is an integer.  Also note that we now use the ++ operator to increment ct; that looks more elegant to users of C but is equivalent to the previous version.

Clear[mandelbrotRate];
mandelbrotRate  = 
        Compile[{cx, cy, {lim,  _Integer}},  
          Module[{c, z, ct = 0}, c = cx  +  I cy; z  =  c; 
            While[(Abs[z] < 2.0)  &&  (ct ≤ 
                    lim), ++ct; z  =  z*z +  c;]; 
            Return[ct]] ];

It is also useful to save calculations in a data array so that we can vary the display format without repeating those calculations.  We now obtain higher resolution in similar or shorter time.  Note that the data table is arranged to produce a plot with the same orientation.

MandelbrotData  =  
        Table[mandelbrotRate[cx, cy, 50], {cy, -1.2, 1.2, 
            0.01}, {cx, -2, 0.6, 0.01}];

The display below uses ListDensityPlot, where the option MeshRange assigns numerical values to array elements so that we can read coordinates from the plot.  Click on the figure, then hold down the Ctrl key while moving the cursor — you should see cursor coordinates in the lower-left corner of the window.  We also used a skeletal function for the ColorFunction option; note here that the parentheses are needed to define the extent of that function.

ListDensityPlot[MandelbrotData, Frame -> False, 
        Mesh -> False, 
        MeshRange -> {{-2, 0.6}, {-1.2, 1.2}}, 
        AspectRatio -> Automatic, 
        ColorFunction -> (Hue[Log[1 + #]] &)];

Having an algorithm that works, we should include proper documentation before moving on to other projects.  You may find the internal comments anal-compulsive, but comments are essential to understanding more complicated functions.

\(\(\(Clear[mandelbrotRate];\)\(
\)
    \)\)
\(\(\(mandelbrotRate::usage\  = \ \*""mandelbrotRate[cx,cy,lim] computes the escape rate for c = cx + I cy under the quadratic map z->\!\(z\^2\)+c beginning with z=0.  The maximum number of iterations is given by lim, which must be an integer."";\)\(
\)
    \)\)
\(\(mandelbrotRate\  = \ 
        Compile[{cx, cy, {lim, \ _Integer}}, \ 
	
          Module[{c, z, ct = 0}, 
 (*\ 
              define\ point\ in\ complex\ plane\ *) 
		c = 
              cx\  + \ I\ cy; 
 (*\ initialize\ map\ *) 
		z\  = \ 
              c; 
 (*\ 
              iterate\ quadratic\ map\ until\ z\ exceeds\ 2\ or\ iteration\ limit\ is\ reached\ *) 			
		While[\((Abs[z] < 
                    2.0)\)\  && \ \((ct ≤ 
                    lim)\), 
			\(++ct\); 
			z\  = \ 
                z*z + \ c;]; 
 (*\ 
              return\ iteration\ counter\ *) 							
		Return[
              ct]]\ ];\)\)

?mandelbrotRate

Produce an appealing display of seahorse valley, which is the Mandelbrot set in the vicinity of c  ≈  -0.745  +  0.107  I.  Explore this area with ever increasing resolution and have some fun.

Optional arguments

We have already encountered many built-in functions in which one or more arguments may be omitted.  We have also created functions that have definitions that apply for various numbers of arguments.  It is often useful to specify default values for missing arguments.  Mathematica distinguishes between two types of optional arguments.  Positional arguments are interpreted according to the order they appear in the sequence of arguments, while options are replacement rules which may be placed after the positional arguments in arbitrary order.  Positional arguments may be omitted, and replaced by pre-defined defaults, beginning from the end and moving forward.  Thus, in the design of the function you must carefully consider which arguments the user is most likely to omit.  Options, on the other hand, may be given in any order or number.

Optional arguments with default values

The basic syntax for defining a function with an optional positional argument is 
	f[x_,⋯,a_:def,⋯]:= rhs
where x_,⋯ represents the set of mandatory arguments and a_:def represents an optional argument that is assigned a default value def if omitted.  Optional positional arguments should be placed after mandatory arguments and before optional replacement rules.  If there are several optional positional arguments, arguments are dropped from the end of the list.  In the example below, the optional exponent is given a default value of unity.

Clear[f];
f[x_, n_:  1]  :=  (x^n)

{f[π],  f[π, 2], f[f[π, π], π], 
      f[π, f[π, π]]}

Supplying options with defaults to your functions

The pattern x___ using a triple underscore indicates a sequence of zero or more objects that collectively will be named x; this pattern can be used to indicate that an arbitrary number of options may be present.  It is usually a good idea to specify explicitly that the options must occur in the form of replacement rules by using the more precise pattern specification opts___Rule.

For example, suppose that we wish to create a function which plots some quantity of interest.  Rather than deciding all of the plotting options in the definition of that function, we can give the user the flexibility to choose those options without revising the function.  Recall our earlier definition of a function that plots spherical Bessel functions.

plotsbess[n_, xmin_, xmax_] := 
      Plot[{sbessj[n, x], sbessy[n, x]}, {x, xmin, xmax}]

If sbessj and sbessy are not currently defined, go back to sbess and activate them before proceeding and use the back button to return here.  As it stands, this version does not display sbessj very well for n > 0.

plotsbess[5, 0, 20];

We can improve this function by allowing the user to select options, such as the vertical scale.

Clear[plotsbess];
plotsbess[n_, xmin_, xmax_, opts___Rule] := 
      Plot[{sbessj[n, x], sbessy[n, x]}, {x, xmin, xmax}, opts]

With one option we can control the scale

plotsbess[5, 0, 20, 
        PlotRange -> {Automatic, {-1, 0.5}}];

and with several more we can add labels and other embellishments.

\(\(plotsbess[5, 0, 20, PlotRange -> {Automatic, {\(-0.2\), 0.2}}, 
        GridLines -> Automatic, Frame -> True, 
        FrameLabel -> {"x", \*""(Subscript[𝒿, 5]),(Subscript[𝓎, 5])""}, 
        PlotLabel -> "Isn't this better?"];\)\)

Default values for the options to be associated with the symbol f are established using the Options[f] command.  Note that we clear any existing options and attributes as well as the function definition using ClearAll instead of simply Clear, which would have left the options intact.  We also use With to assign default options to a local constant included in the Plot function.  To apply this technique to a built-in function that expects a sequence and will not accept a list, simply apply Sequence to the list of options; namely defaults=Sequence@@optionList.

\(\(\(ClearAll[plotsbess];\)\(
\)
    \)\)
\(\(Options[
          plotsbess]\  = \ {PlotLabel -> "Spherical Bessel Functions", \ Frame -> True, \ 
          FrameLabel -> {"x", \*""(Subscript[𝒿, ℓ]),(Subscript[𝓎, ℓ])""}, \ 
          GridLines -> Automatic};\)\ 
\)
plotsbess[n_, xmin_, xmax_, opts___Rule] :=  
      With[{defaults = Options[plotsbess]}, 
        Plot[{sbessj[n, x], sbessy[n, x]}, {x, xmin, xmax}, opts,  
          defaults]]

Notice that we give plotsbess the user options opts first, followed by the default options; then, if the function is called with an option in the argument list for which a default has been defined, the argument takes precedence over the default.  We now obtain a well-formatted plot without supplying as many arguments.

\(\(plotsbess[5, 0, 20, PlotRange -> {Automatic, {\(-0.2\), 0.2}}, 
        FrameLabel -> {"x", \*""(Subscript[𝒿, 5]),(Subscript[𝓎, 5])""}];\)\)

Finally, if we want the default options to use arguments of the function when it is evaluated, it is better to define our own version of the options list using delayed assignment and pattern arguments.  In the example below the Bessel function order is included in the default plot and frame labels by using an options list that is evaluated with the plotting function.  You can, of course, override the defaults by including appropriate replacement rules in the argument list.

\(\(\(ClearAll[plotsbess];\)\(
\)
    \)\)
\(\(options[plotsbess, 
          n_] := \ {PlotLabel -> "Spherical Bessel Functions for ℓ=" <> ToString[n], \ Frame -> True, \ 
          FrameLabel -> {"x", {Subscript[𝒿, n], 
                Subscript[𝓎, n]}}, \ 
          GridLines -> Automatic};\)\ 
\)
plotsbess[n_, xmin_, xmax_, opts___Rule] :=  
      With[{defaults = options[plotsbess, n]}, 
        Plot[{sbessj[n, x], sbessy[n, x]}, {x, xmin, xmax}, opts,  
          defaults]]

plotsbess[5, 0, 20, 
        PlotRange -> {Automatic, {-0.2, 0.2}}];

This technique of specifying default options while permitting the user (including yourself, of course!) to choose the options most appropriate to a particular application will give your functions considerable flexibility and versatility.  We strongly urge you to use this technique liberally.

Modifying default options for built-in functions

In the preceding example we used With to supply the default options because we chose to associate the options with our function, plotsbess, rather than with the built-in function Plot.  If we had changed the options for Plot instead, then it would not have been necessary to include the default options; however, that procedure is risky because it can affect the operation of the built-in function later in undesired ways.  Nevertheless, we illustrate the technique here for the sake of completeness even though we do not recommend this procedure.

The function

SetOptions[symbol,option1->value1,option2->value2,⋯]

sets a sequence of default options associated with a particular symbol to specified values.  Because SetOptions expects a sequence, we apply Sequence to the list of options we want for plotsbess.  Note that we save the default options for Plot first so that we can reset them before leaving this section.

saveOptions[Plot] = Options[Plot];
\(\(\(SetOptions[Plot, Sequence @@ Options[plotsbess]];\)\(
\)
    \)\)
plotsbess[n_, xmin_, xmax_, opts___Rule] :=  
      Plot[{sbessj[n, x], sbessy[n, x]}, {x, xmin, xmax}, opts]

We now obtain a formatted plot without supplying as many options

plotsbess[5, 0, 20, 
        PlotRange -> {Automatic, {-0.2, 0.2}}];

but can still override the defaults

\(\(plotsbess[5, 0, 20, PlotRange -> {Automatic, {\(-0.2\), 0.2}}, 
        FrameLabel -> {"x", \*""(Subscript[𝒿, 5]),(Subscript[𝓎, 5])""}, 
        PlotLabel -> "override plot label"];\)\)

However, we must remember to reset the options for Plot — otherwise, the options appropriate to plotsbess would contaminate all other plots made during a session in which the options for Plot were modified.

SetOptions[Plot, Sequence @@ saveOptions[Plot]];

More examples

Interesting examples require too much code and space for inclusion within this notebook, but several can be accessed using the hyperlinks provided below.  These examples illustrate many of the programming techniques introduced in this notebook.

cicada.nb — explores the fascinating phenomenon of synchronized emergence of insect populations.

logistic map.nb — explores the competition between order and chaos in discrete models of nonlinear dynamics.

cobweb.nb — produces cobwebbing diagrams which illustrate the behavior of iterated function systems, such as the logistic map.  Some techniques of linear stability analysis are introduced.

\(holdstuff = NestList[HoldForm[2\^#] &, 2, 4]\)

stuff = holdstuff //. HoldForm[x_] -> x

ilg /@ stuff

(∫_a^b) 
          f dx  ≈  ((b - a)/(3  n))[
        (Subscript[f, 0]) + 4  (Subscript[f, 1]) + 2  (Subscript[f, 2]) + 4  (Subscript[f, 3]) +  ⋯  + 
          2  (Subscript[f, n - 2]) + 4  (Subscript[f, n - 1]) + (Subscript[f, n])]

\(\((δ f)\)\^2\  = \ \ ∑\+i\((\ \(∂f\/∂x\_i\) δ x\_i)\)\^2\)