Departments of Civil Engineering and Mathematics
CE 109: Computing for Engineering --- Mathematica Session 1

Mathematica: A System for Doing Mathematics by Computer

Engineers were some of the first professionals to recognise the benefits of using computers to perform calculations. Over the years, as hardware and software have become immensely more powerful and computers more and more widely available to individuals, the main thing that's changed is what we mean by "calculations". Once computers were reserved for processing large amounts of numerical data (what sometimes gets called "number-crunching"). While that's still one of the things they're good at, the story doesn't end there any more.

Mathematica is just one of the mathematical software tools available to the modern professional engineer: one of the more powerful, flexible and expensive ones, as it happens! It's described by its manufacturers as "A System for Doing Mathematics by Computer". In this introductory session, you'll explore what that means in practice.

The price you pay for all this power and flexibility is that Mathematica takes some learning. Unlike some word processing packages, for example, or even some of the more restricted mathematical software you sometimes see, its use isn't self-explanatory. It has a syntax all its own, and it can be very pedantic about errors in that syntax. Using it often feels more like programming a computer than using a package, and actually, that's no accident. Stick with it, though: it's worth it in the end.

1. Arithmetical calculations

To start with, it may be helpful to think of Mathematica as a kind of calculator. For example, try typing

2 + 2

To get Mathematica to perform this calculation, hold down the "shift" key and press "return" ("return" without "shift" just gives you a line break.) With any luck, you should get the output

We say that the expression 2 + 2 has been "evaluated".

This calculation may have taken a long time. That's not because Mathematica finds it difficult, but because the part of the program that actually does the calculating doesn't get loaded until the first time it's needed; it's this loading process that takes the time. All subsequent calculations should be much quicker, some very fast indeed. Here are a few

6 + 2

6 - 2

6*2

6 2

6/2

6^2

100!

1000!

Sin[Pi/3]

Sqrt[9]

Sqrt[50]

2^(1 + 4)

Notice the use of the symbols *, / and ^. Notice, too, the way Mathematica interprets the space between the 6 and the 2.

Exercises and explorations 1

1.	The last few examples all involve the use of brackets of some kind. But sometimes square brackets have been used, and sometimes round ones (parentheses). Does this matter? What are the rules governing what kind of brackets to use when?

2. 	Would it have mattered if we'd typed

sin[pi/3]

instead? What are Mathematica's rules about capital and lower-case letters?

3. 	Was the output from

Sqrt[50]

what you'd expected? What's happening here? What other examples can you find?

The Basic Input palette

If you look to the left-hand side of the screen , you may see a tall rectanglar window, the top portion of which looks like this:

If not, you can open this window (which is known as the "Basic Input" palette) by finding Palettes in the File menu and seclecting BasicInput.

This facility allows you to create Mathematica commands that look like typeset mathematics. Try, for example:

(6^2)

Sin[(π/3)]

(Sqrt[50])

Exercises and explorations 2

1.	What happens if you type 6, followed by "ctrl-6" (hold down the "Ctrl" key and hit 6), followed by 2? What does "ctrl-6" do? How about "ctrl-/"?

2.	 Try typing "esc, p, i, esc" (hit the escape key, then p, then i, then escape again). What happens? What if you leave out the 'i'? What other symbols can you create in this way?

Exact or approximate?

All the commands above are examples of exact arithmetic: calculations with whole numbers, surds, fractions and named constants such as π. This is something a calculator doesn't do: it works with numbers in the form of decimals, or in standard scientific form. Mathematica can also do this type of arithmetic. Try the following pairs of commands, taking note of any differences in the output:

(2/6)(2.0/6.0)

Cos[(π/6)]Cos[(3.1415927/6)]

Sin[2]Sin[2.0]

(Sqrt[50])(Sqrt[50.0])

Exercises and explorations 3

1.	Experiment with the Mathematica command N, as in

N[π]

N[π, 100]

N[(2/6)]

How would you describe what this command does? How could you see it being used?

2.	Investigate the Mathematica commands Floor, Ceiling and Round, as in

Floor[3.77]

Round[25]

2. Symbolic calculation

Where Mathematica starts to differ from most calculators is in its ability to perform symbolic calculations, for example:

Expand[(((x + y))^3)]

Factor[(x^3) + 3  (x^2)*y + 3  x*(y^2) + (y^3)]

Apart[(x/((1 + x) (1 + (x^2))))]

Together[%]

(Notice the use of the percent sign, %, to mean "last output").

D[(x^3), x]

∫3  (x^2) \[DifferentialD]x

Exercises and explorations 4

1.	What do Apart and Together do? How about Cancel and Simplify?

2.	Open up the Algebraic Manipulation palette (under Palettes in the File menu). This palette, unlike Basic Input, has the setting "Evaluate in Place". To find out what this means, first type, without evaluating,

∫((2 + 3  x + (x^2))/(2 + 2  x + (x^2))) 
        \[DifferentialD]x

Then select the fraction inside the integral, and click on the Apart[■] button. With the same area selected, click on Together[■]. Try using Expand[■] on the numerator, and so on.

3.	Using the Basic Input palette, type

(∑_(r = 1)^6)((1/(r^2)))

Express this sum as a decimal. Try summing from 1 to n, and from 1 to ∞ (you'll find the infinity symbol on the palette).

3. Graphics

Mathematica incorporates a wide range of two- and three-dimensional graphics functions. The simplest is Plot, which generates ordinary two-dimensional graphs, as in

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

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

(Notice here the use, for the first time, of curly brackets or braces, to enclose the elements of the x-range or to group collections of functions. More about this later in the course.)

You can take control of some of the characteristics of the plot by means of what are called option settings, for example:

Plot[Sin[x], {x, -π, π}, PlotRange -> {-1.5, 1.5}]

Plot[Sin[x], {x, -π, π}, 
      PlotRange -> {{-2 π, 2  π}, {-1.5, 1.5}}]

Plot[Sin[x], {x, -π, π}, AspectRatio -> 1]

Plot[Sin[x], {x, -π, π}, AspectRatio -> Automatic]

Plot[{Sin[x], Cos[x]}, {x, -π, π}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1]}]

More than one option can be set in this way, for example:

Plot[(E^(-x)), {x, 0, 3}, Background -> Hue[0.5],  
      PlotLabel -> "Decaying exponential function"]

For a complete list of options for Plot, together with their default settings (what Mathematica assumes the setting to be if you don't set the option explicitly), type

Options[Plot]

Exercises and explorations 5

1.	If you have not done so already, look up the entry for Plot in the Help system: choose Help from the Help menu, and type "Plot" in the search box. Find out how the AxesLabel option works, and produce a graph of tan(x) where the x-axis is labelled "x" and the y-axis "tan(x)".

2.	Produce a plot showing the relationship between the graphs of y = (x^2) and y = (Sqrt[x]). Your graphs should be in different colours and the scales should be the same on both axes.

4. Using Mathematica to look at beam deflections

In many engineering situations, piecewise-defined functions have an important role. The simplest way to define a function is to use one single rule that applies for all values of x, like this:

f(x) = (x^2) - 2  x.

A piecewise-defined function, on the other hand, is something like the following:

g(x)={({{-1, , x < 1, }, {(x^2) - 2  x, , x ≥ 1. }})

With a piecewise-defined function, then, there are two, three or several rules, each applying for different sets of values of x.

Mathematica allows you to define your own functions. In the case of ordinary, single-rule functions like f, this is what you type:

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

Note the use of the underscore symbol and the "colon-equals" sign. When you type this and "shift-enter", you'll get no output, unusually for Mathematica. That's because you're teaching Mathematica a new function rather than asking it to perform a calculation, From now until the end of your session--or until the machine crashes, whichever is sooner--you can use this new function. For example

f[0]

f[2]

Plot[f[x], {x, -3, 3}, PlotRange -> All]

Piecewise-defined functions (at least the simplest kind of them) can be defined in Mathematica using what's called an If-expression:

g[x_] := If[x < 1, -1, (x^2) - 2  x]

(read this as "if x is less than 1, then g(x) = -1, otherwise g(x)  =  (x^2) - 2  x").

These can be used in exactly the same way:

g[0]

g[2]

Plot[g[x], {x, -3, 3}, PlotRange -> All]

Note that piecewise-defined functions don't have to be continuous; the basic "step function" is a key example:

s(x)={({{0, , x < 0, }, {1, , x ≥ 0. }})

In Mathematica:

s[x_] := If[x < 0, 0, 1]

Plot[s[x], {x, -3, 3}, PlotStyle -> RGBColor[1, 0, 0]]

One application of piecewise-defined functions occurs in structural mechanics, as part of  "Macaulay's method" for solving beam deflection problems. The following is an example of a typical beam displacement formula (for the beam problem described below):

y=(1/EI)(((20 (W - 12500))/3)x+((25000 - W)/15)(x^3)-(250/3)(x^4)+((50000 + 7 W)/30)(((x - 10))^3)+(250/3)(((x - 10))^4)).

What's unusual about this, mathematically speaking, is that the terms (((x - 10))^3) and (((x - 10))^4) are actually ramp functions: that is, they're understood to take the value 0 for all values of x less than 10. Sometimes, though not always, this is indicated by the use of square brackets:

y=(1/EI)(((20 (W - 12500))/3)x+((25000 - W)/15)(x^3)-(250/3)(x^4)+((((50000 + 7 W)/30)[x - 10])^3)+(((250/3)[x - 10])^4)).

To represent a ramp function using Mathematica we could use an If-expression. The following is one way:

rampFunction[x_, a_, r_] := If[x < a, 0, (((x - a))^r)]

Plot[rampFunction[x, 10, 3], {x, 9, 11}, PlotRange -> All]

We'll apply the above Mathematica techniques to the following beam deflection problem: a beam of length 14 m and flexural stiffness 10 MN (m^2) is simply-supported at the points x = 0 and x = 10; a distributed load of intensity 2 kN/m is applied between the two supports, and a variable concentrated load W is applied downwards at the right-hand end of the beam:

(Note: positive deflections are taken upwards here, opposite to convention, since this simplifies the plotting later on.)
The bending moment analysis yields the following formula for y as a function of x and W (we first "clear" any previous values from the variable names x, y and W):

Clear[x, y, W]

flexstiff = (10^7)

y[x_] := 
      (1/flexstiff) 
        (((20 (W - 12500))/3) x + (((25000 - W))/15) (x^3) - 
            (250/3) (x^4) + (((50000 + 7  W))/30)*rampFunction[x, 10, 3] + 
            (250/3)*rampFunction[x, 10, 4])

You can make plots of y using Plot in the usual way, provided that you give W a value:

W = 3500; Plot[y[x], {x, 0, 14}]

(N.B. the use of ";" to separate the two commands.)

A function for y'(x) can be defined as follows (instead of typing this from scratch, you may want to make a copy of the y function definition and modify that):

yprime[x_] := 
      (1/flexstiff) 
        (((20 (W - 12500))/3) + (((25000 - W))/5) (x^2) - 
            (1000/3) (x^3) + 
            (((50000 + 7  W))/10)*rampFunction[x, 10, 2] + 
            (1000/3)*rampFunction[x, 10, 3])

Exercises and explorations 6

1.	Plot graphs of y(x) for values of W = 0, 2, 4, 6, 8 kN.

2.	For each of the values of W in question 1, obtain values for:
	  (i)  the maximum absolute displacement of the beam, and the point at which it occurs;
	  (ii)  the gradient of the beam at the point of support x = 10.
	Write up your answers in the table below.

Note: since you have expressions for both y(x) and y'(x) you could in principle determine the values in part (i) exactly by solving the equation y'(x) = 0, and so on. However, it is somewhat easier, and just as good if you only need 1 or 2 significant figures accuracy, to estimate values from graphs. Try the following techniques for help in reading graphs: (i) enlarge the graphic output: click on the graphic, move the mouse to the "resize" blob at the bottom right-hand corner, click (and hold) on the blob and drag the mouse to the right; (ii) use the option GridLines->Automatic in Plot; (iii) click on the plot graphic and hold down the <Ctrl> key -- the cursor's shape will change to a "+" -- as you move the cursor around, its coordinate position can be read off in the bottom left corner of the Mathematica window.

3.	The second derivative of y, y''(x), can be defined as the following Mathematica function:

ydprime[x_] := 
      (1/flexstiff) 
        (((2 (25000 - W))/5) x - 1000  (x^2) + 
            (((50000 + 7  W))/5)*rampFunction[x, 10, 1] + 
            1000*rampFunction[x, 10, 2])

Using ydprime, estimate the position of the (interior) point of contraflexure, where the curvature changes sign and y''(x) = 0.

({{W/kN, \(Maximum\ absolute
\ displacement/mm\), \(\(Point\ of\ maximum\)\ 
	displacement/m\), Gradient at x = 10, \(\(Point\ of\)
contraflexure/m\)}, {0
	,  ,  ,  ,  }, {2
	,  ,  ,  ,  }, {4
	,  ,  ,  ,  }, {6
	,  ,  ,  ,  }, {8
	,  ,  ,  ,  }})