Chapter 1
Introduction

For convenience we turn off the spell checker by typing

Off[General::spell1]

1.1 Preliminary Remarks

Many vibration problems in engineering are nonlinear in nature. The usual linear analysis may be inadequate for many applications. An essential difference in the study of nonlinear systems is that general solutions cannot be obtained by superposition, as in the case of linear systems. Moreover, the nonlinearity brings many new phenomena, which do not occur in linear systems. To study nonlinear systems, one has to learn new mathematical techniques, which have been developed in many branches of applied mathematics, physics, and engineering. In the past several years, a number of powerful computer software packages have been developed that allow one to perform complicated symbolic manipulations. One of the best known is Mathematica; it is created by Wolfram Research, Inc.. In this book, we use Mathematica's symbolic programming techniques to implement various perturbation methods for studying the dynamics of weakly nonlinear systems. Instead of being burdened by the tedious algebra required to obtain the solutions, Mathematica enables us to focus our attention on understanding the techniques and the physics, thereby free our time for creative thinking.

Symbolic manipulators are powerful tools. They can be used to perform analyses in either an interactive or an automatic way. Learning Mathematica in the process of studying perturbation methods should enhance the concepts of Mathematica programming techniques as well as various perturbation techniques.

In Chapter 2, we describe how various perturbation methods, such as the Lindstedt-Poincaré technique, the method of multiple scales, and the method of averaging, can be used to obtain uniform asymptotic expansions of the solutions of the undamped and unforced Duffing equation. In Chapter 3, we apply the Lindstedt-Poincaré technique, the method of multiple scales, the method of averaging, the generalized method of averaging, the Krylov-Bogoliubov-Mitropolsky technique, and the method of normal forms to obtain uniform expansions of the solutions of an undamped and unforced single-degree-of-freedom oscillator with quadratic and cubic nonlinearities.

In Chapter 4, we use the method of multiple scales, the generalized method of averaging, and the method of normal forms to obtain uniform expansions of the solutions of a forced single-degree-of-freedom oscillator with quadratic and cubic nonlinearities. In Chapter 5, we apply the method of multiple scales, the method of normal forms, and the generalized method of averaging to study higher-order approximations for systems with internal resonances. In Chapter 6, we apply the method of multiple scales to obtain uniform expansions of the solutions of externally and parametrically excited linearly coupled or uncoupled multiple-degree-of-freedom systems.

In Chapter 7, we focus our attention on continuous systems with cubic nonlinearities. In Section 7.1, we discuss the concepts of adjoint and solvability conditions. Then, we use the method of multiple scales to construct uniformly valid asymptotic expansions of the solutions of hinged-clamped beams in Section 7.2 and cantilever beams in Section 7.3. In Chapter 8, we use the method of multiple scales to directly attack three continuous systems with quadratic and cubic nonlinearities. We treat buckled beams in Section 8.1, circular cylindrical shells in Section 8.2, and near-square plates in Section 8.3.

In Chapter 9, we show how one can use the method of multiple scales to construct higher-order approximations of the solutions of continuous systems having two-to-one internal resonances. We discuss two-mode interactions in buckled beams in Section 9.1 and four-mode interactions in suspended cables in Section 9.2.

1.2 Symbolic Computation

The Mathematica programming language allows not only for interactive manipulations of a wide range of objects, but it also encourages the automation of complicated computations. Mathematica actually works by using systems of rewrite rules, and the user can freely create and use his or her own systems of rules. This distinguishes it from traditional programming languages, which normally have no such features. We have applied rule-based and functional programming techniques in Mathematica throughout this book. Hence, we briefly cover the basic ingredients of rule-based programming in Sections 1.2.1--1.2.5 and some functional programming features in Sections 1.2.6--1.2.8.

For further information on Mathematica, consult the standard Mathematica reference: The Mathematica Book, Stephen Wolfram, Third Edition, Wolfram Media and Cambridge University Press, 1996. Full descriptions of the functions in the standard Mathematica packages are available in Mathematica 3.0 Standard Add-On Packages manual, Wolfram Media and Cambridge University Press, 1996. Alternatively, simply click on the Help menu to view the Help Browser.

In the Text cells of this book, we use bold-magenta to indicate user-defined functions or names and bold-black to indicate Mathematica's built-in functions.

1.2.1 Global Rules

Global rules are applied whenever the appropriate left-hand side is encountered or matched. The advantage of using global rules is that one does not need to apply explicitly a rule to get things done. However, it is almost impossible to prevent such global rules from being applied to a certain expression. There are two types of user-defined global rewrite rules: those using = and those using :=. The distinction between them lies in when the right-hand side is evaluated. Moreover, each of these rules has two forms, depending on where the rule is stored, resulting in four rules indicated by =,^=,:=, and ^:=.

1.2.1.1    =  rules

We can think of  = , an assignment statement in analogy with traditional procedural programming languages, as a (global) rewrite rule. The characteristic property of using the rule = is that the right-hand side is evaluated immediately when the assignment is made, and all subsequent matched left-hand sides are replaced by the evaluated right-hand side. For example, let us consider the rule

a = 2;

The output is suppressed by ending the input with a semicolon. From now on, whenever a is encountered in any expression, Mathematica will replace it by its value 2.

expr1 = (a^2) + 3 - a + Tan[a]

Let us try to evaluate expr1 for another value,

a = 3;

expr1

The new value of a is not recognized in expr1.

Let us try to clear the definition of a by using Clear or

a =.

Then, let us evaluate expr1,

expr1

Again, expr1 is not changed accordingly.

To be able to evaluate expr1 for different values of a, we need to either define expr1 before the assignment of a or use ":=", instead of "=". The precise meaning of := is taken up in the next section. Here, we illustrate the former choice.

Clear[a]; expr1 = (a^2) + 3 - a + Tan[a]

a = 2; expr1

a = 3; expr1

In traditional programming languages, the left-hand side of an assignment statement is required to be a simple identifier (i.e., a symbol). In computer algebra, the left-hand side can be arbitrarily complicated. For example,

look[here + 9] = Expand[(1 + x)^2 + y]

Note that the output of an = expression is the evaluated form of the right-hand side. Mathematica regards the left-hand side as a pattern. Whenever Mathematica finds something that matches this pattern, it replaces the something by the evaluated right-hand side.

look[there + look[here + 5 + 2^2]]

In this evaluation, the pattern look [here + 5 + (2^2)] simplifies to look [here + 9] , which is replaced by 1 + 2 x + (x^2) + y. The resulting expression does not match any pattern involving look and so it is left in unevaluated form. This rule is stored with look.

?look

Again, we see the evaluated form on the right-hand side.

There are some problems associated with left-hand sides that are not symbols. For example, suppose we try to make the following rule:

x + y = z

We get an error message stating that Plus is protected, which means new rules cannot be added for Plus. Every time Mathematica encounters Plus, it searches through the rules for Plus to see if anything applies. If we add a new rule for Plus, then that rule would have to be examined at every subsequent addition. When a rule of the form x + y = z is given, Mathematica interprets it as a rule of the form Plus[x, y]=z. Rules have to be stored somewhere and the default rule is that for the head of the left-hand side. One way of adding  a new rule for Plus is to unprotect Plus, make the rule, and then reprotect it.

\(Unprotect[Plus]; 
x + y = z; 
Protect[Plus]; \)

Now whenever Mathematica sees x + y, it rewrites it as z.

x + s + y + t

Alternatively, we can use UpValues to associate the rule with the unprotected argument of the left-hand side. (Definitions that attach a value to the head of the left-hand side are called DownValues of the head.) For example,

s + t ^= u

Note the caret ^ before the  =  sign. This rule is associated with the symbol s or t.

?s

t + s + t

A given symbol can have both up and down values. Let's give s a down value in addition to the up value it already has.

s[x_] := Exp[I a x];

Then looking at s shows both kinds of values:

?s

Finally, we can access the up and down values individually.

{UpValues[s], DownValues[s]}

Now let us try naively define the sine function using an  =  rule according to

m[x] = Sin[x];

It works properly for the symbol x but not for anything else.

{m[x], m[y], m[1]}

This is where the special symbol _ comes in. The form x_ means a pattern named x. We can show the internal representation of any expression, such as x_ , by either using

FullForm[x_]

or using

x_ // FullForm

An underscore _ in a pattern matches anything, so it is a dummy variable or a "wild card". If it appears on the left-hand side of " = " rule with a name, like x, then the left-hand side is rewritten as the right-hand side with x replaced by the desired variable. Let us use this to redefine m; that is,

m[x_] = Sin[x];

Now m works for any argument:

{m[x], m[y], m[1], m[s - t]}

Thus, one can use  =  rules to define functions.

1.2.1.2    :=  rules

Rules using  :=   do not evaluate the right-hand side immediately but instead leave it unevaluated until the function is actually used. They can be used with simple left-hand sides or with left-hand sides containing patterns. As an example, let us define and compare two rules using  =  or  := .

\(Clear[x]; 
lhs1 = Expand[\((x + 1)\)^2]; 

    lhs2 := Expand[\((x + 1)\)^2]; \)

If these are evaluated, they give the same result.

{lhs1, lhs2}

If we now give a value to x, then lhs1 and lhs2 will use this value in different ways.

x = p + q;

{lhs1, lhs2}

If the left-hand side of a  :=  rule contains a pattern, then on a subsequent occurrence of the left-hand side with actual arguments, the formal (dummy) arguments (or names of patterns) on the right-hand side are replaced by the actual arguments and then the right-hand side is evaluated. Thus, each time the left-hand side of such a rule matches something, it is replaced by a new evaluation of the right-hand side. To see the difference, we again set up the two rules

f[k_] = Expand[k^2]

g[k_] := Expand[k^2]

We note that the right-hand side of the rule for f is evaluated immediately when it is defined, whereas the right-hand side of g is kept unevaluated. Now, try out these two definitions on the same value.

{f[p + q], g[p + q]}

Because the right-hand side of the rule for f was evaluated immediately when it was defined, there is nothing to expand, and k is replaced with p+q without expanding the result. On the other hand, the right-hand side of the rule for g was not expanded when defined and it retains the whole expression Expand[k^2]. When the two functions are subsequently used, f[p + q] is just replaced by (((p + q))^2), while g[p + q] is replaced by Expand[(p + q)^2] , which evaluates to (p^2) + 2  p q + (q^2). The internal representation of such a definition has the following form:

FullForm[g[k_] := p]

To see the actual head of an expression before it is evaluated, we wrap the expression with Hold and then request its FullForm

FullForm[Hold[g[k_] := p]]

Thus, the symbol  :=  is the infix form of SetDelayed. This makes dramatically clear the distinction between the evaluation when the rule is given and the evaluation when the rule is used.

When the head of the left-hand side is protected, we can use UpValues to associate the rule with the unprotected argument of the left-hand side. For example, when we try to differentiate a user-defined definite integral with constant bounds, we obtain

Clear[int, f, x]

D[int[f[x], {x, 0, 1}], x]

which is not what we intended. We can either unprotect D and add a rule or use  ^:=  to assign the definition with int; that is,

D[int[x_, a_], y_] ^:= int[D[x, y], a]

Now we check to see if Mathematica gives us the answer we want.

D[int[f[x], {x, 0, 1}], x]

As a second example, we define Cos[n π]=(((-1))^n) and Sin[n π]=0 for integer n:

\(Unprotect[Cos, Sin]; 

    Cos[n_\ π] := \((\(-1\))\)\^n /; IntegerQ[n]\)Sin[n_ π] := 0 /; IntegerQ[n]Protect[Cos, Sin];

If we declare m to be an integer (here UpValues is used),

IntegerQ[m] ^= True;

then

{Cos[m π], Sin[m π]}

1.2.2 Local Rules

Local rewrite rules (lhs -> rhs  or  lhs :> rhs) are useful for making substitutions without making the definitions permanent. These rules are applied to an expression using the operation  /.  (means ReplaceAll) or  //.  (means ReplaceRepeated in Mathematica).

1.2.2.1    ->  rules

We can use Solve to obtain the solutions of some algebraic equations. For example,

\(eq1 = u\^2 - 3\ u + 2 == 0; 
sol1 = Solve[eq1, u]\)

We can check the result by substituting sol1 into eq1 and obtain

eq1 /. sol1

The output of Solve is a list of lists of local rules. Thus, u -> 1 is a local rule, which is the analog of the global rule u = 1. The rule is applied to an expression by using /. , so the expression eq1/.sol1 means "use the rewrite rules u -> 1 and u -> 2 in eq1". The result of this is {0 == 0, 0 == 0}, which is then evaluated as {True, True}. The usual form of the right-hand side of /. is a list of local rules for some of the symbols that appear on the left-hand side.

u v w /. {u -> 3, w -> 5}

If the right-hand side is a list of lists of local rules, then the result is a list of modified expressions, one for each substitution in the list.

u v w /. {{u -> 3, w -> 5}, {v -> 2, w -> 3}}

Local replacement rules also work on expressions which have a Hold wrapped around them. For example, without a Hold, (∂_x) Cos[x] immediately evaluates to -Sin[x].

(∂_x) Cos[x]

However, the rewrite rule goes under the Hold to make the replacement.

Hold[(∂_x) Cos[x]] /. Cos -> Sin

Local rules with " -> " share with " = " rules the property that their right-hand sides are evaluated at the time they are defined.

1.2.2.2    :>  rules

The local analog of a " := " rule is a " :> " rule; i.e., a local rule that evaluates its right-hand side only when it is used. For example,

expr1 = Exp[I ((Subscript[ω, 2]) - (Subscript[ω, 1])) (Subscript[T, 0])];

Clear[a]

expr1 /. Exp[a_] -> Exp[Expand[a]]

{expr1 /. Exp[a_] :> Exp[Expand[a]], expr1 // ExpandAll}

where the immediate rule (lhs -> rhs) fails to do the expansion. The left-hand sides of the local rules involve patterns rather than just symbols. The difference is that the immediate rule for Exp[a_] replaces it by the evaluation of Exp[Expand[a]] which equals Exp[a]. Hence, when this is used with a equal to I((Subscript[ω, 2]) - (Subscript[ω, 1])) (Subscript[T, 0]), we get the result \(TraditionalForm
      \`E\^\(I\ T\_0\ \((\(-ω\_1\) + ω\_2)\)\)\). On the other hand, the delayed rule (lhs :> rhs) for Exp[a_] replaces it by the unevaluated Exp[Expand[a]] which, when used with a equal to I((Subscript[ω, 2]) - (Subscript[ω, 1])) (Subscript[T, 0]), gives \(TraditionalForm
      \`E\^\(\(-I\)\ T\_0\ ω\_1 + I\ T\_0\ ω\_2\)\). Alternatively, for this simple example, we can use ExpandAll to obtain the same result.

We can check how these expressions are represented internally.

FullForm[Hold[p /. q -> r]]

FullForm[Hold[p /. q :> r]]

Thus, /. is the infix form of ReplaceAll, the arrow  ->  is the infix form of Rule, and the arrow  :>  is the infix form of RuleDelayed, corresponding to Set and SetDelayed for  =  and  := .

1.2.2.3    /.  and  //.

There is another form of  /.  given by  //.  which applies a local rule repeatedly until the expression no longer changes. Internally,  //.  is represented by

p //. q -> r // FullForm // HoldForm

An example of the difference between  /.  and  //.  follows. This example uses a list of rules rather than just a single rule. When a list of rules is applied to a single expression, then the rule for each symbol is tried from the left until a match is found. In the following example, the right-hand side of the  /.  expression consists of a list of two rules for the same symbol, fac. This list is searched from the left until a pattern is found that matches the left-hand side of the  /.  expression. In the first case using  /. , as soon as a match is found, the evaluation is finished. In the second case using  //. , the rules are tried repeatedly from the left on the output of the previous evaluation until no matches are found.

fac[5] /. {fac[1] -> 1, fac[i_] -> i fac[i - 1]}

fac[5] //. {fac[1] -> 1, fac[i_] -> i fac[i - 1]}

In the first case, the left-hand side of the rule fac[1] -> 1 does not match anything in fac[5], but fac[i_] -> i fac[i - 1] does with i_ equal to 5, so the output is 5 fac[4]. In the second case, the left-hand side of the rule fac[i_] -> i fac[i - 1] continues to match a part of the existing expression until 120 fac[1] is obtained. Then the left-hand side of the rule fac[1] -> 1 matches, leading to 120*1,  which simplifies to 120 where neither rule matches, so the output is 120.

If such rules are given globally, then the order in which they are given does not matter since Mathematica will put the more specific rule, fac[1] = 1, first. However, in a list of local rules, applied with  //. , we are completely responsible for the ordering. Thus, the following gives the wrong answer:

fac[5] //. {fac[i_] -> i fac[i - 1], fac[1] -> 1}

If several local rules are given for different symbols, then these rules are applied simultaneously. For example,

{u, v, w} /. {u -> v, v -> w, w -> s}

If the substitutions are carried out sequentially, then the results are quite different.

{u, v, w} /. {u -> v} /. {v -> w} /. {w -> s}

In particular, this means that variables can be interchanged without introducing an intermediate temporary variable.

{u, v} /. {u -> v, v -> u}

1.2.3 Pattern Matching

Mathematica is a language based on pattern matching. Patterns refer to the structure of expressions. The symbol _ by itself, without any symbol on the left, can be used to describe a pattern. (Recall that the FullForm of _ is Blank[] ). For instance, the expression _^_ matches anything of the form x^y, where x and y are any (dummy) expressions.

Clear[a, x]

\(g1[_, _^_] := u; 
{g1[a, a^a], g1[a, b^c]}\)

A pattern of the form x_ is matched by any expression, which is then used for evaluating the right-hand side.

\(g2[x_^y_, z_] := u[x, y, v[z]]; 
{g2[a^a, a], g2[a^b, c]}\)

\(g3[x_^y_, x_] := u[y, v[x]]; 
{g3[a^a, a], g3[a^b, c], g3[a^b, a]}\)

A pattern of the form x_h is matched by any expression whose head is h.

\(g4[x_^y_Integer, z_] := u[x, y, v[z]]; 
    
{g4[a^a, a], g4[a^3, arg2]}\)

If we give a rule for an expression involving two separate arguments with underscores, then we are constructing a function of two variables. Such a function only works if it is given exactly two arguments.

\(g5[x_, y_] := x + y; 
{g5[a], g5[a, b], g5[a, b, c]}\)

Besides rule-schemes using a single underscore_, there are rule-schemes using double and triple underscores. A double underscore, __, is matched by any number of  arguments, excluding zero arguments, separated by commas, while a triple underscore, ___, is matched by any number of arguments, including zero arguments. The form x__ means a sequence of one or more arguments, named x, and x__Head means a sequence of one or more arguments, named x, with the head of each argument is Head. Similarly, the form x___ means zero or more arguments, named x, and x___Head means zero or more arguments, named x, with the head of each argument is Head.

\(g6[x__] := Length[{x}]; 
{g6[], g6[a], g6[a, b], g6[a, b, c]}\)

\(g7[x___] := Length[{x}]; 
{g7[], g7[a], g7[a, b], g7[a, b, c]}\)

Default values and double and triple underscores are important techniques, which give optional number of arguments to functions; that is, a variable number of arguments can be given to such functions. Alternatively, there is another rule in which a specific argument can be optional.

\(g8[x_^y_, z_] := u[x, y, z]; 
{g8[a^b, c], g8[a, c]}\)

Because a = a^1, one might think that g8[a, c] should match the pattern g8[a^1, c], which would mean that it should be rewritten as u[a, 1, c], but of course Mathematica cannot guess that this is what we intend. However, there is a provision to take care of such default values. To assign a default value v to a pattern, one can write _ : v. So, the effect we wanted to achieve is given by the form

\(g9[x_^\((y_:  1)\), z_] := u[x, y, z]; 
{g9[a^b, c], g9[a, c]}\)

In this case, the default value 1 for the exponent is the natural and obvious choice, and Mathematica knows this. It has standard built-in default values for a number of such positions. The notation _.  tells Mathematica to use the built-in default value. Note the almost invisible period after the underscore. Thus, the effect we wanted at the beginning is given by a tiny modification of the original form.

\(g10[x_^y_. , z_] := u[x, y, z]; 
{g10[a^b, c], g10[a, c]}\)

There is still another way that optional arguments occur in Mathematica. Some functions, such as Plot, can take named optional arguments such as AspectRatio -> 1.   By incorporating such functions into our own definitions, we can use these named optional arguments too. Consider an example of a plotting function.

plotwithSin[func_, var_, opts___] := 
      Plot[{Sin[var], func[var]}, {var, 0, 2  Pi}, opts, 
        PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1]}];

plotwithSin[Cos[2  #] &, x];

plotwithSin[Cos[2  #] &, x, AspectRatio -> 1];

In an expression of the form g[x_^n_Integer, z_], the patterns are named with the symbols x, n,  and z, but there is no name for the whole compound pattern x_^n_Integer. There is a way to give names to such compound patterns so that they can be referred to directly on the right-hand side. The syntax consists of a name followed by a colon followed by the compound pattern.

\(g11[expr : x_^y_Real, z_] := z\ expr; 
    
{g11[\((a + b - 2)\)^0.5, 2], g11[\((a + b - 2)\)^3, 2]}\)

1.2.4 Using Patterns in Rules

Patterns play an important role in both global and local rules.

1.2.4.1 Patterns in Global Rules

For example, it is very easy to give rewrite rules for differentiating polynomials of one variable.

\(diff[x_^n_. , x_] := n\ x^\((n - 1)\); 

    diff[a_ + b_, x_] := diff[a, x] + diff[b, x]; \)

Notice the default value for n in the first rule. Try it out on some typical functions.

{diff[y, y], diff[z^2.1, z], diff[u^2 + v^3, u]}

But notice that diff does not know what to do with a constant times x, or just a constant for that matter, and we have no obvious way as yet to teach it what to do.

{diff[3, x], diff[3 x, x]}

We could try the following: first, give a rule for products.

diff[a_ b_, x] := a diff[b, x] + b diff[a, x]

Using this for a x gives

diff[a x, x]

The program does not know that a is supposed to be a constant, so we have to tell it that explicitly, with a last rule.

diff[a, x] = 0;

Then it gives the "correct" answer.

diff[a x, x]

However, this is not very satisfactory. We would like some general way to say that a is not a function of x. Section 1.2.5 below will continue this discussion.

1.2.4.2 Patterns in Local Rules

We will encounter this type of application quite often in the following chapters. In the following example, we show how to change the appearance of a matrix. Using Array to make a matrix with indexed entries, we then use a local rule to display the indices as subscripts.

Array[x, {3, 5}] /. x[i__] -> (Subscript[x, i]) // MatrixForm

where the two consecutive underscores next to i on the left-hand side of the rule represent one or more arguments.

We can also use pattern matching to replace the coefficients of various terms any way we want. For example,

Clear[f]

eq1=u''[t]+u[t] f'[(Subscript[x, 0])]+(1/2) u[t]^2 f''[(Subscript[x, 0])]+(1/6) u[t]^3 Derivative[3][f][(Subscript[x, 0])]==0;

eq2 = eq1 /. _.  u[t]^b_.  -> (Subscript[k, b]) u[t]^b

Sometimes we may obtain, in a bifurcation analysis, a data set, which includes a set of bifurcation points (Subscript[x, i]) and the corresponding control parameters cpSearch and cpFix, such as {{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[x, 4]),cpSearch},cpFix}. If we want to switch these two parameters and then use the same data format as our initial input to a program, we have at least two choices. First, we take apart the data list and then assemble them:

data1 = {{(Subscript[x, 1]), (Subscript[x, 2]), (Subscript[x, 3]), (Subscript[x, 4]), cp1}, cp2};

{Append[Drop[data1[[1]], -1], data1[[2]]], 
      data1[[1, -1]]}

Second, we simply apply a local rule using patterns to the data list as a whole:

data1 /. {{a__, b_}, c_} -> {{a, c}, b}

1.2.5 Restricting Pattern Matching with Predicates

So far, all of the rules we have considered have been "context free" rewrite rules. Whenever the pattern is matched, the rewriting is carried out. There can be a restriction on the head of the matching expression included in the pattern. However, there are also conditional rewrite rules, which are only applied when some condition is satisfied. First, we have to discuss predicates (functions that return the value True or False) in Mathematica because all of the conditions will always be expressed in terms of them. All built-in predicates that are defined for all expressions end with Q. It is easy to display all of them.

?*Q

Specifically requesting information on our function gives the usage message.

?OddQ

Predicates are used to control pattern matching. In general, they are applied using  /;  , which is the infix form of Condition.

FullForm[Hold[x /; y]]

However, the position of the predicate in an expression allows it to be used in different ways. If the predicate is placed at the end of a global rule definition, it will restrict the application of the rule. For instance, we can use predicates to extend our definition of differentiation in Section 1.2.4.1 to deal with arbitrary polynomials in a very natural way by adding a single conditional rule.

Clear[diff]

\(diff[x_^n_. , x_] := n\ x^\((n - 1)\); 

    diff[a_ + b_, x_] := diff[a, x] + diff[b, x]; 

    diff[a_\ b_, x_] := a\ diff[b, x] + b\ diff[a, x]; 

    diff[a_, x_] := 0 /; FreeQ[a, x]\)

Now constants and products are handled properly

diff[a[x] + b x^2 + c, x]

Using rule /; Predicate restricts the rule to those situations in which the predicate evaluates to True;  that is, to those expressions belonging to the type given by the predicate. An unrestricted rule is the same as a conditional rule where the predicate always equals True.

If a predicate is immediately applied after a pattern, it will restrict the pattern matching rather than the rule application. For instance,

fac1[1] = 1; fac1[n_] := n fac1[n - 1]

This will work perfectly well if positive integers are given as arguments.

fac1[5]

However, if another argument is given, then it will fail badly.

{fac1[-3], fac1[endless]};

A very large output is suppressed. What happens, of course, in these cases is that the value 1 is never encountered as an argument, so the function keeps calling itself recursively until the built-in recursion limit is reached. This bad behavior can be corrected by using a conditional rule.

fac2[1] = 1; fac2[n_] := n fac2[n - 1] /; n > 1;

{fac2[5], fac2[-2], fac2[something]}

One can also use the form _?Predicate, which restricts the pattern to something for which the predicate evaluates to True.

\(fac3[1] = 1; fac3[n_?Positive] := n\ fac3[n - 1]; 
fac4[1] = 1; 
    fac4[n_ /; Positive[n]] := n\ fac4[n - 1]; \)

{fac3[5], fac3[-2], fac4[5], fac4[-2]}

In n_?Positive, Positive is a pure function, while in the form using  /;  , the condition is the value of the predicate for the name of the pattern. In either case, Positive or Positive[n] is a positive test in the sense that the pattern is matched and the rule applied only if the test succeeds. Now try fac4 on a real number and see what happens.

fac4[3.6]

The rule is applied until the negative value -0.4 is reached, where the condition fails so fac4[-0.4] is returned in unevaluated form. Of course, the real problem is that we only intend the function to apply to integers. But this additional restriction can easily be added.

fac5[1] = 1; fac5[n_Integer?Positive] := n fac5[n - 1];

{fac5[5], fac5[3.6], fac5[-2]}

The predicate that appears after  /;  or ? can also be a user defined expression.

test[x_Integer?((# > 2) &)] := x + 1

{test[1], test[2], test[3], test[4], test[5]}

The form ?Predicate can only be used after single slots, but the form  /; Predicate can be used after any pattern, simple or compound. For instance,

\(f1[a_, b_] /; EvenQ[a + b] := a^b; 

    f1[a_, b_] /; OddQ[a + b] := a^\((\(-b\))\); \)

{f1[1, 3], f1[2, 3], f1[3, 3]}

Predicates also play an important role in manipulating lists.

Select[Range[-5, 5], Positive]

Cases[Range[-5, 5], _?Positive]

Cases[Range[-5, 5], a_ /; a > 0]

There is also another form of Cases in which an operation is applied to the entries that are selected.

Cases[Range[-5, 5], (a_ /; a > 0) :> Sqrt[a]]

The following example was the 1992 Mathematica programming competition question. The problem is to write a function called maxima that starts with a list of numbers and constructs the sublist of the numbers bigger than all previous ones from the given list. For instance, maxima[{3, 6, 4, 2, 8, 7, 9}] should return {3, 6, 8, 9}. The winning entry used a pattern with a condition in a local rule.

maxima[list_List] := 
      list //. {a___, x_, y_, b___} /; y <= x -> {a, x, b}

maxima[{3, 6, 4, 2, 8, 7, 9}]

1.2.6 Pure Functions

Given an expression expr involving a variable x, we can think of expr as describing a function with x being the argument. If we need to refer to this function by name, we can either define a rule for f in the form f[x_] :=  expr or use the object Function[x, expr]. Either f or Function[x, expr] is a name for this function, and we can use the two interchangeably. To apply them to an argument, we write f[arg] in the familiar way or in the so-called pure function form Function[x, expr][arg].

f[x_] := 1 + x^2g = Function[x, 1 + x^2];

{f[2], g[2]}

The name of the formal argument in a pure function does not matter. Function[x, x^2] is the same function as Function[y, y^2]. This fact is easy to see if you apply such a function to an argument, the result is the same in both cases.

{Function[x, x^2][a], Function[y, y^2][a], Sin[a]}

Note that Sin and all built-in functions, with no arguments, are pure functions. Because the names of the variables in a pure function do not matter, Mathematica denotes these variables with symbols. The symbols #1, #2, ... are used for the first, second, ... variable. The internal form of #i is Slot[i]. If we use these forms, then we do not give the arguments of Function[ ] in order to declare the names of the variables. Thus, instead of Function[x, x^2] we simply write Function[#1^2]. A convenient abbreviation for #1 is #, so we can simplify our example a bit more and write Function[#^2]. Finally, there is a postfix operator & for Function. That is, body& is the same as Function[body]. Using this operator, we arrive at the shortest form #^2 & that our example can take. We have used this form quite frequently in this book. Let us see the internal form of this short form of a pure function.

FullForm[#^2 &]

The operator & has a very low priority, just above assignment. Therefore, x =  body & is understood as x =  Function[body], but x ->  body & is interpreted as Function[x ->  body]. If we want the right-hand side of the rule to be the pure function, we use x ->  (body &). Beware of x ->  (body) &. Another case where the low priority of pure functions requires the use of parentheses is in predicates for patterns: x_?body& is not the same as x_?(body&). The latter is usually correct.

1 + # + #^2 & [3]

The whole expression to the left of & is part of the body of the pure function. No parentheses are necessary. In other situations, parentheses are necessary.

h[a] /. h -> (#^2 &)

Parentheses around the pure function are needed here; otherwise, the whole rule would be considered a part of the body of the function.

Pure functions written in the form with # and & are very concise and useful, especially, in functional operations, such as Map, Apply, Thread, and iteration functions, such as FoldList. For example,

eq1 = (∂_(t, t))u[x, t]/a + u[x, t] == F[x]/a

Clear[lhs2]

\(eq2 = lhs2 == rhs2; 
eq3 = lhs3 == rhs3; \)

We can multiply both sides of eq1 by a and use Thread and obtain

eq1 // Thread[a  #, Equal] & // ExpandAll

Moreover, we can add eq2 to eq3 by

eq2 + eq3 // Thread[#, Equal] &

1.2.7 Map and Apply

In Mathematica, two important commands that take functions as arguments are Map[ ] and Apply[ ]. The various versions of Map act on the arguments of an expression, while Apply acts only on its head.

1.2.7.1 Mapping Functions onto Expressions

The operation Map[f, list] or in infix form f /@ list applies the function  f  to each element of the list list. The second argument of Map[ ] need not be a list, however. Any expression of the form h[(Subscript[e, 1]), (Subscript[e, 2]),  ... , (Subscript[e, n])] will do. The result of the mapping is the expression h[f[(Subscript[e, 1])], f[(Subscript[e, 2])],  ... , f[(Subscript[e, n])]].

Clear[f]

{Map[f, a + b + c], #^2 & /@ {1, 2, 3}}

Map[f, expr, levelspec] has an optional third argument that specifies the levels at which to map. The default level is {1};  that is, at the elements of expr. In a matrix, the entries are at level 2 (there are two levels of lists). If we want to map a function at these entries, we can use Map[f, matrix, {2}].

expr1 = {{(Subscript[a, 1, 1]), (Subscript[a, 1, 2])}, {(Subscript[a, 2, 1]), (Subscript[a, 2, 2])}};

Map[entry, expr1, {2}] // MatrixForm

This is a very powerful facility and is one of our main tools in manipulating expressions. The purpose of the operations based on Map is to make it possible to treat lists as wholes. For instance, a really poor way to square the entries in a list is as follows:

\(list1 = {a, b, c}; 
Table[list1[\([i]\)]^2, {i, Length[list1]}]\)

The term list1[[i]] tears apart the original list by extracting its parts one at a time, ^2 squares each part, and then Table reassembles the parts into a new list. Map, however, allows us to treat mathematical structures as wholes so that we do not tear them apart and rebuild them again.

Map[ ] always maps the function at all elements of the given levels. The command MapAt[f, expr, poslist] allows you to apply a function at selective positions.

MapAt[G, a b + c d + e f, {2, 1}]

expr2 = a + b/a + c Exp[a + 1]

If we want to map a function f at all occurrences of a, we first find all the positions of a in our expression.

Position[expr2, a]

This list of positions is in the right form for MapAt and can be used directly.

MapAt[f, expr2, %]

This particular example could have been done more easily with a replacement rule.

expr2 /. a -> f[a]

In some applications, the function to be performed on an element (Subscript[e, i]) may depend on its position i. The operation MapIndexed[f, {(Subscript[e, 1]), (Subscript[e, 2]),  ... , (Subscript[e, n])}] behaves essentially like Map[ ], but it passes the position of each element as a second argument to f. The resulting expression is

MapIndexed[1 + #^2 &, {a, b, c, d}]

MapIndexed[1 + #^2 &, {{a, b}, {c, d}}, {2}]

Therefore, the function that is mapped must be a function of two arguments. The second argument is a position (note the list braces) that can be used to modify the operation performed on the first argument.

MapIndexed[#1^First[#2] &, {a, b, c, d, e}]

1.2.7.2 Apply

The rule Apply[h, expr] or h @@ expr replaces the head of expr by h.

a + b + c // FullForm

Replacing the head Plus by Times gives the product of the three terms in the sum above.

{Apply[Times, a + b + c], Times @@ (a + b + c)}

With a third argument, Apply[h, expr, levelspec] replaces heads in the parts of expr described by levelspec by h.

Apply[b, {{a[1, 1], a[1, 2]}, {a[2, 1], a[2, 2]}}, {2}]

This simple definition finds the average of a list of numbers

average[list1_List] := Apply[Plus, list1]/Length[list1]average[{a, b, c}]

Note that computing the average of a list in this way requires no do-loops or knowledge of the length of the list. Apply is frequently used if one first wants to prepare a number of ingredients and then apply some operation to them. The ingredients can be held in a list until they are ready and then the head of the list is changed to the appropriate operation by using Apply.

To generate more concise notations for the outputs created in the following chapters, we can introduce displayRule to meet our needs. For example,

output=(Subscript[u, 2]^((2, 0, 0)))[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])]+(ω^2) (Subscript[u, 2])[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])]==2 (Subscript[u, 1]^((1, 1, 0)))[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])]-δ (Subscript[u, 1])[(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])]^2;

\(\(displayRule = {
        \(\(Derivative[a__]\)[u\_i_]\)[__] :> 
          SequenceForm[
            Times\ @@\ MapIndexed[D\_\(#2[\([1]\)] - 1\)\%#1\ &, {a}], u\_i], 
        u\_i_[__] -> u\_i}; \)\)

output /. displayRule

1.2.8 Nest and Fold

Nest and its related operations NestList and FixedPoint apply a function to its argument many times. Nest[function, x, n] applies function to x and repeats the application n times

Nest[(((1 + #))^2)&, x, 3]

NestList[(((1 + #))^2)&, x, 3]

Here is an example producing a simple continued fraction.

Nest[1/(1 + #) &, a, 3]

An operation that is closely related to Nest is FixedPoint, which nests its operation until there is no change. For instance, everyone is familiar with what happens if the Cos key on a pocket calculator is pushed repeatedly. In principle, FixedPoint is what happens if it is pushed forever.

{Nest[Cos, 0.3, 5], Nest[Cos, 0.3, 10], FixedPoint[Cos, 0.3]}

Actually, FixedPoint stops after machine accuracy is achieved. We can use the option SameTest to get some control of FixedPoint.

FixedPoint[Cos, 0.3, SameTest -> (Abs[#1 - #2] < 10^-5 &)]

The second pair of functions, Fold and FoldList, do something similar to Nest and NestList, but for functions of two variables.

Fold[f, u, {x, y, z}]

FoldList[f, u, {x, y, z}]