Introduction to Mathematica

by

Kenneth Shiskowski
Eastern Michigan University

for

Principles of Linear Algebra
Using Mathematica

by

Kenneth Shiskowski and Karl Frinkle

Chapter 1: The Very Basics About Using Mathematica

Mathematica is an extremely powerful mathematical software package which incorporates text editing, mathematical computation and programming as well as 2D and 3D graphics capabilities.  You can literally write a mathematical textbook completely in Mathematica which includes all of the text and graphics in one smoothly flowing document.  If you have never or only slightly used Mathematica before, then it will take some effort to learn how it works - believe me that it is well worth the time expended for the ability to do mathematically almost anything you can dream of that a computer might be able to do for you.  In this Introduction to Mathematica, you will see but a fraction of its capabilities but hopefully enough to get you well on your way in doing calculus, linear algebra, differential equations and basic mathematical programming as well as computation and graphics.

This chapter will discuss the fundamentals of using Mathematica for the novice user.  If you are already familar with Mathematica you may wish to skip this chapter although we warn you to do so at your own risk.  The new user of Mathematica will find it quite difficult in the beginning, but with practice and patience you will master all of the basics and in time come to enjoy using Mathematica.

Mathematica files are called notebooks, and in a notebook you can place text along with input commands and their associated outputs which can be literally anything such as graphics, tables or lists, functions, etc.  You can group the material in a notebook into different types of cells which are indicated on the right side of the notebook by brackets.  At the top of the notebook you will see the tab Palettes and under it is the Writing Assistant which will allow you to create new cells and/or modify cells.  You can use Writing Assistant to change the font, color and size of the text in your cells and you can also do this using Format at the top of the notebook.  The word processing capabilities of Mathematica are very similar to that of Microsoft Word.

Each chapter of this notebook file is actually a section in Mathematica with the first cell of the section is subtitle cell which can be created by creating a horizontal ba between or just after a cell and then choosing Subtitle from the pull down menu at the very top left of the lower ruler at the top of your screen.  In order to create a subsection of this section (or chapter), do the same as just described but choose the Subsection from this menu. 

If you wish to delete a cell (use Ctrl X) or modify its entire contents in some way, then click on the cell tag or bracket on the right and then carry out the desired operation using Writing Assistant or the tabs at the top of the screen.  In text, in order to create a new paragraph in a cell, use Enter.  To do the same in an Input cell where commands are placed, use Enter as well.  If you wish to split a cell, then use Shift Enter with the cursor at the location of the split.  In using Writing Assistant or any pull down tab, if you click on the triangles on the left you will open or close one of the sections inside this tab.  Note that text paragraphs are not necessarily indented automatically so that you must indent them yourself manually if you want this to happen. 

If you wish to close a group of cells and only see the first cell of the group (which is probably the title cell of the group), then double click on the far right bracket for the group.  You will then see a cell bracket with an arrow to the right of the cell bracket of the title or first cell of the group.  If you double click on this arrow, then you can open all the cells of this group.  The copy Ctrl C and paste Ctrl V features of Mathematica are the same as those of Microsoft Word and other software.  If you wish to change the size, font or other feature of a collection of cells, select one of them by clicking on its bracket and then hold the Ctrl key down while you select the rest of the cells - now go to the Format Tab or other location and carry out your change.

It is strongly recommended that you save (use Ctrl S) your work constantly since like all software, Mathematica can clitch which could cause you to loose some or all of your material.  As well, you should have backup copies of all of your work on a separate computer or flash drive since from my own personal experience I know that unfortunate problems can occur.

If you are using Mathematica to do homework problems, it is strongly suggested that you place each problem in a single group of cells with the first cell the title of the problem.  This will make it much easier to organize your work both for yourself and the instructor who may read your material.  After you have finished working in a particular Mathematica notebook, it is also recommended that you delete all of your output from the file unless it will take too long to recompute it.  Most, if not all of the size of a Mathematica file will be due to graphics, especially 3D graphics, and such files can become very large and consequently take Mathematica quite a while to open or save, and at such times a clitch can happen.  Under the tab Cell, you have the command to Delete All Output which removes all output from the entire file - you might use this periodically while working in a notebook in order to shorten the file.  When you reopen a notebook where all output has been deleted you can reconstruct it all by going to the tab Evaluation and using the command Evaluate Cells.

If a Mathematica calculation is taking too long and/or you notice that there is an error in the input, then in order to terminate the calculation you should go to the Evaluation tab at the top of the screen and select Abort Evaluation.  This should immediately halt the calculation in its tracks unless Mathematica is stuck in some enormous loop and can not find its way out - then your only alternative is to turn your computer off, that is, gently pull the plug on it.

Never use capital letters to define a quantity in Mathematica as it might already be a built in command name which you can NOT override with something else.  As well, do not use the capital letters D and N for any kind of variable in Mathematica as they are also command names.  The commands Clear and ClearAll will undefine a quantity which you have named.  If you use the command Exit[ ] it should clear everything from memory which you have defined by quitting the Mathematica kernel which is the guts of Mathematica. 

Finally, if you are a novice or beginner at using Mathematica, then besides this introductory notebook there are many videos on YouTube which explain many of the basic features of Mathematica.  It is strongly suggested that you seek these out and hopefully will find a few useful ones for doing what you are interested in.

Chapter 2: Basic Arithmetic and Polynomial Algebra

Basic Arithmetic

In this chapter, we will start to use Mathematica to do some basic arithmetic and algebra computations.  In the arithmetic which is done first, we will add and multiply, factor positive integers into products of powers of primes, find the greatest common divisor GCD and the least common multiple LCM of two positive integers and more.  In the algebra, we will factor polynomials, divide one polynomial into another to get their quotient and remainder, solve for the roots of a polynomial and also solve equations for their unknowns, and more.

In order to create an Input cell where you can do your calculations, go to the tab Palettes and bring up Classroom Assistant.  Now click on the tab Create Input Cell with the cursor at the end of our prior work.  Now you will have a new Input cell as part of your current group of cells.  Both palettes Classroom Assistant and Basic Math Assistant have symbols such as π in them as well as the natural number e. 

Below is typed in the first Input cell the command 1 + 1.  If you hit Shift Enter with the cursor on this line, then Mathematica will carry out your command and produce 2 as an Output.  After you get the output, Mathematica automatically assumes you want another Input cell, and so typing right after an output will be in a new Input cell.  You can also go back an insert a new Input cell by creating a horizontal bar between two Input cells by clicking on the region between the two cells, and then using Alt 9 while Alt 7 creates a text cell.  You can insert a Text cell between Input cells by creating the horizontal line divider between cells by clicking on the space between the cells and then using Text Cell out of the tab Text Cells in Writing Assistant.  As well, if you type in a Mathematica command name such as FactorInteger, but now you have forgotten precisely how it works and need its Help file, then put the cursor in the command name and hit the F1 key.  

Note that Mathematica will also recognize a space in a product as multiplication, although for safety sake you might want to put in all of your multiplications directly.  Mathematica gives exact answers in a calculation if the inputs are all also exact values, but any value with a decimal point in it is treated by Mathematica as an approximation and it gives an approximate answer back.  The command N[V, m] will give the approximate value of the quantity V to m digits of accuracy.  Mathematica can answer most computational questions to an arbitrary number of digits of precision - look up the command WorkingPrecision to see how it can be done as an alternative to the use of the command N.  In the last example of this arithmetic section we also multiply three complex numbers using Product - a complex number in Mathematica is expressed as a + b I for a and b real numbers. 

In order to get a power of someting in Mathematica which is placing a superscript, use Ctrl 6 together to get an exponent location after you have already typed in the base.  If you want a fraction in the same standard way, then type the numerator followed by Ctrl / together and then the denominator in the location created.

1+1

Sum[k^2,{k,1,10}]

Sum[π^k,{k,1,3}]

N[Sum[π^k,{k,1,3}],10]

Sum[3.14159265^k,{k,1,3}]

Product[1/k,{k,1,5}]

FactorInteger[90]

(2^1) (3^2) (5^1)

GCD[210,90]

LCM[210,90]

QuotientRemainder[83,5]

5 16+3

Polynomial Algebra

Now we switch from arithmetic to algebra.  Our algebra will be mainly polynomial and similar to what we did in the arithmetic part above although we will also find the roots of polynomials as well.  After the first two computations of multiplying out two polynomials and then dividing the one of larger degree by the smaller degree one, we will name or define the two polynomials as Poly1 and Poly2, and then repeat the process to see that Mathematica understands what we want.  Also, we define the list called QR below which is the quotient first and remainder second in our division of Poly2 by Poly1 - a list is an ordered collection of objects which in Mathematica has { } around its elements with commas between the elements.  Then QR[[1]] is the quotient and QR[[2]] is the remainder in the division.  We will discuss lists and matrices in the next chapter.

Mathematica expresses its equations with a double equal sign == while it uses a single equal sign to make a definition or assignment of a quantity to a name such as in the use of Poly1 and Poly2 below.  Hence, 5 x + 3 y == 9 is an equation in Mathematica while Eqn1 = 5 x + 3 y == 9 assigns the name Eqn1 to this equation for the time you are using this notebook unless you decide to change it.  In Mathematica, you do not need to put in a multiplication sign as Mathematica usually places a space between objects which indicates multiplication. 

One last bit of useful information is that Mathematica uses the percent sign % to refer to the last computed output and %% to refer to the next to last computed output.  This can be helpful as we will see below.  We will use the % below when we want to change the roots of Poly1 to a set of rules which can then be substituted back into Poly1 to see that we get 0 (or very close to 0) back.  As well, it is probably better to assign a name to your quantities in order to be better able to use them later and know what you are specifically talking about, and we do this for the roots of Poly2.

Expand[(7 (x^3)+5 (x^2)-9 x+1) (-4 (x^6)-(x^5)+3 (x^4)-2 (x^3)+(x^2)-8 x+5)]

5-53 x+98 (x^2)-16 (x^3)-30 (x^4)-31 (x^5)+6 (x^6)+52 (x^7)-27 (x^8)-28 (x^9)

PolynomialQuotientRemainder[-4 x^6-x^5+3 x^4-2*x^3+x^2-8 x+5,7 x^3+5 x^2-9*x+1,x]

Poly1=7 (x^3)+5 (x^2)-9 x+1

Poly2=-4 (x^6)-(x^5)+3 (x^4)-2 (x^3)+(x^2)-8 x+5

Expand[Poly1 Poly2]

QR=PolynomialQuotientRemainder[Poly2,Poly1,x]

Expand[QR[[1]] Poly1+QR[[2]]]

NRoots[Poly1==0,x]

{ToRules[%]}

Poly1/.%

Poly2Roots=N[Roots[Poly2==0,x],10]

Rules={ToRules[Poly2Roots]}

Poly2/.Rules

Rules[[2]][[1]][[2]]

Poly2Roots[[2]][[2]]

Poly=Expand[-4 Product[x-Poly2Roots[[k]][[2]],{k,1,6}]]

Sum[Re[Coefficient[Poly,x,k]] x^k,{k,0,6}]

Product[k-(k+2) I, {k,0,3}]

Chapter 3: Lists, Sets, Strings, Sequences and Matrices

Now we look more closely at the different ways we can organize information in lists, sets, strings and matrices, and how these different structures work and can be manipulated.  In Mathematica, there are NO sets since Mathematica requires an ordering be placed on its data and so it deals with lists which can be treated as sets if you ignore the order of the elements in the list.  The list L of the elements a, b and c in this order is given as L = {a, b, c} in Mathematica.  Shortly, we will look at taking the union, intersection, complement and concatenation (or joining) of lists as well as taking out parts of a list.  The empty list is { } and L[[k]] is the kth element of list L.

A string S is a grouping or ordered collection of characters or symbols with quotes around them such as S = "Mary had a little lamb."  The string assigned to the name S is the sentence between the two double quotes.  Strings often show up as the title to Mathematica plots and similar structures.

A sequence or table in Mathematica is a list of objects created from a formula such as {1, 4, 9} can be created as the output from the input              Table[(k^2), {k,1,3}].

L1={a,b,1,3,c,5}

L2={7,2,5,2,5,1,c,b}

Union[L1,L2]

Intersection[L1,L2]

Complement[L1,L2]

Join[L1,L2]

Length[L1]

L1[[2]]

Squares=Table[k^2,{k,1,10}]

Sum[Squares[[k]], {k,1,10}]

Now we switch to matrices or two dimensional arrays of rows and columns of entries.  Mathematica considers a matrix M to be a list of lists where each element of M is one of the rows of the matrix M.  As such, M = {{1,2,3}, {4,5,6}, {7,8,9}} has {1,2,3} as its first row, {4,5,6} as its second row and finally {7,8,9} as its third row.  We will look at examples of adding, multiplying, inverting, transposing and finding determinants of matrices where multiplication of matrices is indicated by a dot or period between their names as given in the 5th line of input below.

M = {{1,2,3}, {4,5,6}, {7,8,9}}

MatrixForm[M]

K = {{-5,7,1}, {0,-4,6}, {2,-1,9}}

MatrixForm[K]

MatrixForm[M+K]

MatrixForm[M.K]

Det[M]

Det[K]

Det[M.K]

Inverse[M]

MatrixForm[Inverse[K]]

MatrixForm[Transpose[M]]

Chapter 4: Expressions Versus Functions

This chapter will look at the differences in Mathematica between expressions versus functions and how to manipulate and evaluate both.  You should think of an expression as the rule for some function f, that is, if f(x) = 5 x + 9, then 5 x + 9 is an expression which is the rule of the function f.  Mathematica treats a function very differently than it treats an expression - you should think of an expression as a string of symbols while a function is a string of symbols (its rule or expression) with a method of evaluating its expression at different values of the variable(s) in the expression.

Let's now look at the difference in Mathematica between the expression g = 5 x + 9 and the function f with rule f(x) = 5 x + 9, and how each of them can be evaluated at x = 1.  Note that normal function evaluation can be done so that f evaluted at x = 1 is f[1] while the expression g can be evaluted at x = 1 by the substitution command g /. x -> 1.  If you want to compose the function f with the built-in function sin[x], then f[sin[x]] will do it - this composition can also be done using f@f[x] or Nest[f, f[x], 1], and its result is an expression, not a function.  If you want f to be composed with itself more than once, then use Nest[f, f[x], k] to get f composed with itself k times.

g=5 x + 9
g/.x->1

f[x_]=5 x + 9

f[1]

f[x]/.x->1

g/.x->sin[x]

f[sin[x]]

h=f@f[x]

Simplify[h]

Nest[f,f[x],1]

Simplify[%]

Nest[f,f[x],1]/.x->2

Nest[f,f[x],5]

Simplify[%]

As the last topic of this chapter, let's do an example of a piecewise function and its graph.  A piecewise function is one whose rule is given in parts or pieces where each part is used only when certain conditions are satisfied.  Happily, Mathematica has a Piecewise command we can utilize.  Note that in defining the function f[x] below that && is the Mathematica notation for And in joining two statements.

f[x_]=Piecewise[{{Sin[x],x≤-2},{Cos[x],x>-2&&x≤3},{-x+5,x>3}}]

Plot[f[x], {x,-π,2π}, PlotStyle->{Red,Thick}]

Chapter 5: Plotting 2D and 3D Graphs Along with Animations

Section 5.1 Function, Parametric and Implicit Plotting in the xy-plane

We begin our plotting chapter with the basic 2D plotting of the graphs of functions y = f(x) and expressions as well as plotting parametric curves x = f(t), y = g(t).  We will plot single functions and parametric curves as well as several together and in combination in different colors.  Besides these two types of curves, we will implicitly plot equations in the two variables x and y such as the unit circle with equation (x^2)+(y^2)=1 or something more complicated.  The implicit plotting of equations can be done using the ContourPlot command and ContourStyle option to control color and thickness, etc of the plot.

Let's begin with simple function or expression plotting, and then move on to parametric curve plotting.  In the first example, remember to get base e for the exponential function from the palette Basic Math Assistant.  Also, in order to get a superscript or exponent, use Ctrl 6 together after the base is already in place.  In order to create a fraction, first type the numerator and then hit together the Ctrl and / keys to be able to place the denominator.  

The plot option of PlotStyle can control color and thickness for graphics as well as whether the graph is solid or dashed except when doing implicit plotting of equations using ContourPlot when PlotStyle switches to ContourStyle.

f[x_]=x Sin[(x^2) (E^x)]

Plot[f[x],{x,-1,2}, PlotStyle->Red, PlotRange->{-2,2}]

g=(E^(((-1)/2) Sin[3x]))

Plot[{f[x],g},{x,-1,2}, PlotStyle->{Directive[Red,Thick],Directive[Green,Thick]}, PlotRange->{-2,3}]

Plot[{f[x],g},{x,-1,2}, PlotStyle->{Directive[Red,Thickness[.005]],Directive[Green,Thickness[.01]]}]

ParametricPlot[{Sin[2 t],Sin[3 t]},{t,0,2 π},PlotStyle->Directive[Blue,Thick]]

ParametricPlot[{(E^(Sin[2 t]))Sin[t],(E^(Sin[2 t]))Cos[t]},{t,0,2 π},PlotStyle->Directive[Blue,Thick]]

ParametricPlot[{((E^(Sin[2 t]))-3 Cos[6 t])Sin[t],((E^(Sin[2 t]))-3 Cos[6 t])Cos[t]},{t,0,2 π},PlotStyle->Directive[Blue,Thick]]

ParametricPlot[{{Sin[2 t],Sin[3 t]},{Sin[t],Cos[t]}},{t,0,2 π},PlotStyle->{Directive[Blue,Thick],Directive[Red,Thick]}]

ContourPlot[{(x^2)+(y^2)==1, (x^4)+(y^4)==1},{x,-1,1}, {y,-1,1}, ContourStyle->{Directive[Red,Thick],Directive[Blue,Thick]}]

ContourPlot[{((((x-5))^2)/49)+((((y-10))^2)/144)==1, ((((x-5))^2)/49)-((((y-10))^2)/144)==1},{x,-10,25}, {y,-10,25}, ContourStyle->{Directive[Red,Thick],Directive[Blue,Thick]}]

Section 5.2 Animations in the xy-plane

Now we turn our attention to creating a movie or animation in the xy-plane whose frames consist of plots of the 3 different types from section 5.1.  Let's begin by plotting the function y = sin(x) from x = 0 to x = A where the animation parameter A goes from 0 to 4π.  It is followed by running two y = sin(x) animations at once based on the same animation parameter A.

Animate[Plot[Sin[x],{x,0,A}, PlotRange->{{0,4π},{-1,1}},PlotStyle->Directive[Blue,Thick]],{A,0.01,4π},AnimationRunning->False]

Animate[{Plot[Sin[x],{x,0,A}, PlotRange->{{0,4π},{-1,1}},PlotStyle->Directive[Blue,Thick]],Plot[Sin[x+A],{x,0,4π}, PlotStyle->Directive[Red,Thick]]},{A,0.01,4π},AnimationRunning->False]

Our next animation will work on parametric plots.  Let's use a parametric plot to produce an insect shape where the insect is gliding around the origin in our animation.

r=(E^(Sin[2t]))-4 Cos[12 t]+(1/2) Sin[24 t]

ParametricPlot[{r Sin[3t]+10,r Cos[2 t]},{t,0,48π},PlotRange->{{-25,25},{-25,25}},PlotStyle->Directive[Blue,Thick],PlotPoints->250 ]

Rot={{Cos[A], -Sin[A]},{Sin[A],Cos[A]}}

{{Cos[A],-Sin[A]},{Sin[A],Cos[A]}}

GlidingInsect=Simplify[Rot.{r Sin[3t]+10, r Cos[2 t]}]

Animate[ParametricPlot[GlidingInsect/.A->B,{t,0,48π},PlotRange->{{-25,25},{-25,25}},PlotStyle->Directive[Blue,Thick], PlotPoints->500],{B,0,2π},AnimationRunning->False, AnimationRate->0.4]

Lastly, we will do an animation involving implicit plotting of an ellipse where the center is moving along the circle with center at the origin and radius 10.  In this example, we use the Epilog option to put into each frame of our movie the circle of the ellipses' centers.

ellipses=((((x-10 Sin[A]))^2)/4.)+((((y-10 Cos[A]))^2)/25.)==1

Animate[ContourPlot[Evaluate[ellipses/.A->B],{x,-15,15},{y,-15,15},ContourStyle->Directive[Blue,Thick],PlotPoints->100, Epilog->{Red,Thick,Circle[{0,0},10]}],{B,0,2π},AnimationRate->.15,AnimationRunning->False]

Section 5.3 Function, Parametric and Implicit Plotting in Space

It is time to move function, parametric and implicit plotting from the xy-plane to space.  We begin by plotting functions z = f(x,y) which give surfaces in space followed by parametric curve and surface plotting.  A parametric curve or spacecurve is of the form x = f(t), y = g(t), z = h(t) with one independent variable t while a parametric surface looks like x = f(u,v), y = g(u,v), z = h(u,v) with two independent variables u and v.  Lastly, we do implicit plotting of equations in the 3 variables x, y and z.  The plot of a single equation in the variables x, y and z is a surface.

f[x_,y_]=Sin[x+y] Cos[x-y]+3

g=Sin[x y]+(E^(-((x^2)+(y^2))))

Plot3D[{f[x,y],g},{x,-3,3},{y,-2,2}, PlotStyle->{Red,Blue},AxesLabel->Automatic]

Plot3D[f[x,y]g,{x,-3,3},{y,-2,2}, PlotStyle->Red,Mesh->None,AxesLabel->Automatic]

Now we examine parametric curve plotting in space.  The spacecurve we will do first is a 3D rendering of the insect from section 5.2.

r=(E^(Sin[2t]))-4 Cos[12 t]+(1/2) Sin[24 t]

ParametricPlot3D[{r Sin[3t],r Cos[2 t],r},{t,0,48π},PlotStyle->Directive[Blue,Thick],PlotPoints->250 ,AxesLabel->Automatic]

ParametricPlot3D[{(E^(t/100)) Cos[t],(E^(t/100)) Sin[t],(t/10)},{t,0,48π},PlotStyle->Directive[Blue,Thick],PlotPoints->250,AxesLabel->Automatic ]

ParametricPlot3D[{Cos[2t],Sin[4t],Cos[6t]},{t,0,2π},PlotStyle->Directive[Blue,Thick],PlotPoints->250,AxesLabel->Automatic ]

Next, we look at parametric surface plotting.  The surface we do first is a torus or donut followed by three interlocking mutually perpendicular torii.

ParametricPlot3D[{(7+3Sin[u])Sin[v],(7+3Sin[u])Cos[v],2Cos[u]},{u,0,2π},{v,0,2π},PlotStyle->Blue,AxesLabel->Automatic ]

ParametricPlot3D[{{(7+3Sin[u])Sin[v],(7+3Sin[u])Cos[v],2Cos[u]},{2Cos[u],(7+3Sin[u])Sin[v],(7+3Sin[u])Cos[v]},{(7+3Sin[u])Sin[v],2Cos[u],(7+3Sin[u])Cos[v]}},{u,0,2π},{v,0,2π},PlotStyle->{Blue,Red,Green} ,AxesLabel->Automatic]

The final type of 3D plotting we do is implicit plotting of equations in the 3 space variables x, y and z which generate surfaces.  We start with a surface called a cylinder since it is a surface whose equation has only two of the three space variables in it - this equation is really a curve in the plane of its two variables which is a surface in the third direction.

ContourPlot3D[(x^2)+(z^2)==9,{x,-3,3},{y,-3,3},{z,-5,5},ContourStyle->Blue,BoxRatios->Automatic,AxesLabel->Automatic]

ContourPlot3D[{(x^2)+(z^2)==9,(x^2)+(y^2)==9},{x,-5,5},{y,-5,5},{z,-7,7},ContourStyle->{Blue,Red},BoxRatios->Automatic,AxesLabel->Automatic]

ContourPlot3D[{(x^2)+(y^2)+(z^2)==16,(x^2)+(y^2)==9},{x,-5,5},{y,-5,5},{z,-7,7},ContourStyle->{Blue,Red},BoxRatios->Automatic,AxesLabel->Automatic,Mesh->None]

ContourPlot3D[{(x^2)+(y^2)+(z^2)==16,(x^2)+(y^2)-(z^2)==9},{x,-5,5},{y,-5,5},{z,-7,7},ContourStyle->{Blue,Red},BoxRatios->Automatic,AxesLabel->Automatic,Mesh->None]

Section 5.4 Animations in Space and Using Multiple Animation Parameters

A word of warning in doing 3D animations.  This type of animation can use a great deal of memory, and so it is highly recommended to delete the output from memory before closing the file unless it took a great deal of time to produce the animation and you do not want to have to do it again.

We start with the insect example of the previous sections and we will get it to fly around and flap its wings.  After it, we animate the torus with one animation parameter followed by three animation parameters.

r=(E^(Sin[2t]))-4 Cos[12 t]+(1/2) Sin[24 t]

Animate[ParametricPlot3D[{r Sin[3t],r Cos[2 t],-A r},{t,0,48π},PlotStyle->Directive[Blue,Thick],PlotPoints->250 ,AxesLabel->{x,y,z},PlotRange->{{-6,6},{-6,6},{-6,7}}],{A,-1,0},AnimationRunning->False]

Animate[ParametricPlot3D[{(7+A Sin[u])Sin[v],(7+A Sin[u])Cos[v],Cos[u]},{u,0,2π},{v,0,2π},PlotStyle->Blue,AxesLabel->{x,y,z},PlotRange->{{-10,10},{-10,10},{-2,2}} ,PlotPoints->100],{A,1,3},AnimationRunning->False]

Animate[ParametricPlot3D[{(C+A Sin[u])Sin[v],(C+A Sin[u])Cos[v],B Cos[u]},{u,0,2π},{v,0,2π},PlotStyle->Blue,AxesLabel->{x,y,z},PlotRange->{{-11,11},{-11,11},{-5,5}} ,PlotPoints->100],{A,1,3},{B,1,4},{C,5,8},AnimationRunning->False]

Chapter 6: Solving Single as well as Systems of Equations

Mathematica can solve single equations as well as simultaneous systems of equations, both linear and non-linear.  It can also find approximate or exact solutions although for exact solutions the equation or system must be capable of being solved for exact solutions by some known method.  Also, Mathematica can find both real and complex solutions.  The solution given by the command Solve will give a list of replacements for the variables solved for in your solution.  In order to get the actual values of the solutions as a list and only for a specific variable solved for like z, you must use the input z /. soln where soln is the result of Solve and z is one of the variables solved for in Solve.  The command Solve seeks exact solutions while NSolve always gives approximate solutions.

Single Equations

We begin by solving some single polynomial equations for real and complex solutions, both approximate and exact.  Solve will give all real and complex solutions to a single polynomial equation, that is, it will find all of the real and complex roots of any complex coefficient polynomial.

P=5 (x^2)-7 x -4

soln=Solve[P==0,x]

approxsolns=N[soln,25]

P/.%

xsolns=x/.approxsolns

Expand[5 (x-xsolns[[1]]) (x-xsolns[[2]])]

Q=5 (x^3)-7 (x^2) +I x-4 I

soln=Solve[Q==0,x]

approxsolns=N[soln,15]

Q/.approxsolns

xsolns=x/.approxsolns

Expand[5 (x-xsolns[[1]]) (x-xsolns[[2]]) (x-xsolns[[3]])]

R=-9 I (x^7) -3(x^5)+I (x^4)-4 I (x^3)+(2-I) (x^2)+6x-(I/2)

soln=NSolve[R==0,x]

xsolns=x/.soln

Expand[-9I Product[x-xsolns[[k]],{k,1,7}]]

Systems of Equations

Now we look at some examples of solving systems of equations both linear and non-linear.  We begin by finding the intersection point of two lines in the xy-plane followed by finding the intersection points of two different pairs of circles.

Eqn1=5 x + 3 y ==-4
Eqn2=-7 x + 2 y ==6

Soln=Solve[{Eqn1, Eqn2}, {x,y}]

Eqn1/.Soln
Eqn2/.Soln

approxsolns=N[Soln,10]

xsolns=x/.approxsolns
ysolns=y/.approxsolns

IntersPt={xsolns[[1]], ysolns[[1]]}

ContourPlot[{5 x + 3 y ==-4,-7 x + 2 y ==6, (((x-IntersPt[[1]]))^2)+(((y-IntersPt[[2]]))^2)==(1/16)}, {x,-5,5}, {y,-5,5}, ContourStyle->{Directive[Red,Thick],Directive[Blue,Thick],Directive[Black,Thick]}]

Eqn1=(((x-7))^2) + (((y-2))^2) ==25
Eqn2=(((x-1))^2) + (((y-5))^2) ==16

Soln=Solve[{Eqn1, Eqn2}, {x,y}]

approxsolns=N[Soln,15]

solns={x,y}/.approxsolns

IntersPt1=solns[[1]]
IntersPt2=solns[[2]]

ContourPlot[{(((x-7))^2) + (((y-2))^2) ==25, (((x-1))^2) + (((y-5))^2) ==16}, {x,-5,15}, {y,-5,15}, ContourStyle->{Directive[Red,Thick],Directive[Blue,Thick]}]

ContourPlot[{(((x-7))^2) + (((y-2))^2) ==25, (((x-1))^2) + (((y-5))^2) ==16, (((x-IntersPt1[[1]]))^2)+(((y-IntersPt1[[2]]))^2)==(1/8), (((x-IntersPt2[[1]]))^2)+(((y-IntersPt2[[2]]))^2)==(1/8)}, {x,-5,15}, {y,-5,15}, ContourStyle->{Directive[Red,Thick],Directive[Blue,Thick], Directive[Black,Thick], Directive[Black,Thick]}]

Eqn1=(((x-7))^2) + (((y-2))^2) ==25
Eqn2=(((x+2))^2) + (((y-5))^2) ==16

Soln=Solve[{Eqn1, Eqn2}, {x,y}]

N[Soln,10]

Chapter 7: Basic Differential and Integral Calculus

This chapter will cover some of the basic differential and integral Calculus topics covered in the first (or perhaps second) course in Calculus.  The first section on differential calculus will look at differentiation (difference quotient and derivative rules) and its application to limits using L'Hopital's limit rule, tangent lines and Taylor polynomials.  The second section will look at integration (antidifferentiation rules and Riemann sums) which has applications to geometric and physical quantities.

Differentiation

We begin by graphing the difference quotient ((f(x+H)-f(x))/H) as a function of both x and H for a function y = f(x) and the value of H = (10^(-5)).   This difference quotient will approximate very accurately the derivative function (dy/dx)as H approaches 0, that is, (Underscript[Lim, H->0]) ((f(x+H)-f(x))/H)=(dy/dx).  The command D carries out differentiation with D[f[x], {x, K}] giving the Kth derivative of the function f(x).  The difference quotient will have 0.1 added to it so that its graph will not be identical to the derivative's graph, but shifted up by 0.1 and parallel to the derivative's graph.

f[x_]=Sin[x]

Clear[H]

diffquot[x_,H_]=((f[x+H]-f[x])/H)

D[f[x],x]

Plot[{Cos[x], diffquot[x,(10^(-5))]+0.1}, {x,-2π,2π}, PlotStyle->{Directive[Red,Thick], Directive[Blue,Thick]}]

Animate[Plot[{diffquot[x,-H],Cos[x]},{x,-2π,2π},PlotStyle->{Directive[Red,Thick],Directive[Blue,Thick]}],{H,-3,-(10^(-5))},AnimationRunning->False]

Limit[diffquot[x,H], H->0]

ProductRule=D[a[x] b[x],x]
QuotientRule=D[((a[x])/(b[x])),x]
GeneralPowerRule=D[((a[x])^N),x]
ChainRule=D[Composition[u,v][x],x]
ExponentialRule=D[(a^(g[x])),x]

Let's now graph the tangent line (first order Taylor polynomial) along with some Taylor polynomials to a function y = f(x) at the center x = c.  The command Normal[Series[f, {x, c, N}]] will generate the Nth degree Taylor polynomial to the function f(x) at x = c.  We use the := assignment below to delay computation of the Taylor polynomials while still defining them as a function of their variable x and their degree K.

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

c=1

dydx[x_]=D[f[x],x]

TangLine=y==f[c]+dydx[c] (x-c)

TaylorPoly1=Normal[Series[f[x],{x,c,1}]]

TaylorPoly5=Normal[Series[f[x],{x,c,5}]]

TPoly5=Sum[((D[f[x],{x,K}]/.x->c)/(K!))(((x-c))^K),{K,0,5}]

TPoly5/.x->2

TayPoly[x_,K_]:=Normal[Series[f[z],{z,c,K}]]/.z->x

TayPoly[2.0,5]

Animate[Plot[{TayPoly[x,K],f[x]},{x,-3,5},PlotRange->{-2,2},PlotStyle->{Directive[Red,Thick],Directive[Blue,Thick]}],{K,1,25,1},AnimationRunning->False]

Now we do the last topic under differentiation, L'Hopital's limit rule for the indeterminate form (0/0).

Limit[(((x^2)((Sin[A x])^2))/(Sin[B (x^4)])),x->0]

Limit[((D[(x^2)((Sin[A x])^2),x])/(D[Sin[B (x^4)],x])),x->0]

Limit[((D[(x^2)((Sin[A x])^2),{x,2}])/(D[Sin[B (x^4)],{x,2}])),x->0]

Integration

Integration (or antidifferentiation) is essentially the process of reversing differentiation, that is, an antiderivative of f(x) is a new function F(x) where the derivative of F(x) is f(x) back.  We will see that Mathematica has integration rules built into it as well as sophisticated ways of approximating definite integrals which can give area, volume and many other physical and geometric quantities.

Let's begin by seeing that Mathematica can compute general antiderivatives as well as definite integrals.  Note that in computing an antiderivative, Mathematica does not give an arbitrary constant C added onto it.  The method of partial fractions requires a rewriting of the quotient of two polynomials you want to integrate into new quotients where the denominators are factors of the original denominator, the Apart command in Mathematica will do this rewriting for you.

Integrate[(x^K),x]

Integrate[((x^2)/(1+(x^4))),x]

Integrate[(x/(1+(x^4))),x]

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

Integrate[((x^3)/((5x-1)(x+2))),x]

Apart[((x^3)/((5x-1)(x+2)))]

Apart[((x^7)/((5x-1)(x+2)((((x^2)+1))^2)))]

Apart[((x^3)/((A x+B)(C x+D)))]

Now its time to do a few definite integrals.  If we use the command NIntegrate with a definite integral, then Mathematica approximates its value while otherwise it tries to apply the Fundamental Theorem of Calculus to get an exact value.  You can increase the precision and the number of digits in your answer from NIntegrate by using the option WorkingPrecision.  Let's also do a midpoint (or middle) Riemann sum approximation of a definite integral.

Integrate[Sin[x] (E^(-x)),{x,1,∞}]

NIntegrate[Abs[((Sin[x])/(x^4))],{x,1,10}, WorkingPrecision-> 25]

Integrate[((x^7)/((5x-1)(x+2)((((x^2)+1))^2))),{x,1,2}]

N[%,10]

Let's set up a midpoint Riemann sum for a function y = f(x) on an interval [a, b] for K equal subintervals of [a, b].  This midpoint Riemann sum will be created as a function of f, a, b and K in this order using the assignment command := instead of =.  In order for WorkingPrecision to work in MidPtRiemannSum below, you need to give a, b and f exactly without decimals.

MidPtRiemannSum[f_,a_,b_,K_]:= NSum[(((b-a))/K)f/.x->a+(J-(1/2))(((b-a))/K), {J,1,K}, WorkingPrecision-> 25]

MidPtRiemannSum[((Sin[(E^(x/5))])^4),0,10,1000]

Table[MidPtRiemannSum[((Sin[(E^(x/5))])^4),0,10,L],{L,5,505,50}]

NIntegrate[((Sin[(E^(x/5))])^4), {x,0,10}, WorkingPrecision-> 25]

Plot[((Sin[(E^(x/5))])^4),{x,0,10}, PlotStyle-> {Red, Thick}]

Chapter 8: Solving Single as well as Systems of Ordinary Differential Equations

Now we turn our attention to solving and plotting the solutions to single and systems of ordinary differential equations with or without initial conditions.

Solving a Single Ordinary Differential Equation

We are now interested in solving for, evaluating and plotting the solution function y(x) to a single ordinary differential equation.  You should already (or be about to) study how to solve simple ordinary differential equations by hand when they are exact, first order linear or homogeneous - Mathematica has all of these methods and many more built into it in the command DSolve.  The command NDSolve does numerical solving of a differential equation while DSolve tries to find an exact closed form solution.  In using DSolve or NDSolve, you can indicate a derivative by using the prime notation or using the command D, so y''[x] and D[y[x], {x,2}] are both the second derivative of y[x].

Let's begin by solving and plotting the solution to the initial condition problem y'(x) + 3 y(x) = 5 sin(3x), y(0) = - 9 followed by a 4th order problem.

DSolve[{y'[x]+a y[x]==b Sin[3x],y[0]==c},y[x],x]

soln=DSolve[{y'[x]+3y[x]==5 Sin[3x],y[0]==-9},y[x],x]

Solny=Simplify[y[x]/.soln[[1]][[1]]]

Solny/.x-> 3.0

Plot[Solny, {x,0,6}, PlotStyle-> {Red, Thick}, PlotRange-> {-9,2}]

soln=DSolve[{D[y[x],{x,4}]+3y[x]==5 Sin[3x],y[0]==-9, y'[0]==2, y''[0]==-5, y'''[0]==10},y[x],x]

Solny=Simplify[y[x]/.soln[[1]][[1]]]

Solny/.x-> 1.5

Plot[Solny, {x,0,3}, PlotStyle-> {Red, Thick}, PlotRange-> {-10,100}]

Next, we look at solving an ODE using NDSolve to get a numerical solution to the problem y'(x) + 3 cos(3x) y(x) = 5 sin(3x), y(0) = - 9.  We will then use DSolve to find the exact solution to the same problem.  The exact form solution y[x] to this problem is somewhat time consuming for Mathematica to evaluate and so we create a list of 31 equally spaced points on its graph and plot them instead of a full plot since Mathematica would by default plot at least 50 such points in a full plot.

soln=NDSolve[{y'[x]+2 Cos[3x]y[x]==5 Sin[3x],y[0]==-9},y[x],{x,0,3}]

Solny=y[x]/.soln[[1]][[1]]

Solny/.x-> 2.1

Plot[Solny, {x,0,3}, PlotStyle-> {Red, Thick}, PlotRange-> {-10,1}]

soln2=DSolve[{y'[x]+2 Cos[3x]y[x]==5 Sin[3x],y[0]==-9},y[x],x]

Soln2y=y[x]/.soln2[[1]][[1]]

N[Soln2y/.x-> 1.5]

SolnPoints=Table[{(k/10),N[Soln2y/.x-> (k/10.0)]}, {k,0,30}]

ListPlot[SolnPoints, PlotStyle-> {Red, Thick}, PlotRange-> {-10,1}, Joined-> True]

Solving a System of Ordinary Differential Equations

The first example we look at is the Lorenz system of three coupled ordinary differential equations.  It is an example of sensitivity to initial conditions which says that a very small change in the initial conditions of a system can result in very different (although perhaps similar looking) behavior between the two solutions for each set of initial conditions.  The plot we will get of the solution to the Lorenz system is a parametric curve in space called the Lorenz strange attractor since if the initial conditions are near the attractor then the solution curve is pulled towards it and strange refers to its sensitivity to initial conditions.  The second plot below is the graph of two solutions to the Lorenz system for nearby initial conditions - the two solutions start out close together but over time they end up in different locations although both approach the Lorenz strange attractor.

LorenzSystem={x'[t]==10(-x[t]+y[t]), y'[t]==28x[t]-y[t]-x[t]z[t],z'[t]==((-8)/3)z[t]+x[t]y[t],x[0]==3, y[0]==0, z[0]==1}

SolnLorz=NDSolve[LorenzSystem,{x,y,z},{t,25}, WorkingPrecision-> 10, MaxStepSize->0.005]

ParametricPlot3D[{x[t],y[t],z[t]}/.SolnLorz, {t,0,25}, PlotStyle-> {Red, Thick}]

LorenzSystem2={x'[t]==10(-x[t]+y[t]), y'[t]==28x[t]-y[t]-x[t]z[t],z'[t]==((-8)/3)z[t]+x[t]y[t],x[0]==4, y[0]==-2, z[0]==0}

SolnLorz2=NDSolve[LorenzSystem2,{x,y,z},{t,25}, WorkingPrecision-> 10, MaxStepSize->0.005]

ParametricPlot3D[{{x[t],y[t],z[t]}/.SolnLorz,{x[t],y[t],z[t]}/.SolnLorz2}, {t,0,25}, PlotStyle-> {Directive[Red, Thick],Directive[Blue, Thick]}]

Our second and last example is of the Van der Pol second order ordinary differential equation solved in two different ways - first as a single ordinary differential equation and then secondly as a system of two first order ordinary differential equations.  The Van der Pol second order ordinary differential equation is y''[x] + y[x] - 0.1 (1 - ((y[x])^2)) y'[x] = 0 with initial conditions y[0] = 1, y'[0] = 0 which can also be written as the system of two first order ordinary ODE's given by w'[x] = z[x], z'[x] = 0.1 (1 - ((w[x])^2)) z[x] - w[x] with initial conditions w[0] = 1, z[0] = 0 where we have created the two new functions w[x] = y[x] and z[x] = y'[x].

The two plots below are often called phase portraits as they parametrically plot the solution function y[x] against one of its derivatives.

VanderPolsingle:={y''[x]+y[x]-(1/10)(1-((y[x])^2)) y'[x]==0,y[0]==1,y'[0]==0}

SolnVanderPol1=NDSolve[VanderPolsingle,{y,y'},{x,25}, WorkingPrecision-> 10, MaxStepSize->0.005]

ParametricPlot[{y[x],y'[x]}/.SolnVanderPol1, {x,0,25}, PlotStyle-> {Red, Thick}]

VanderPolsystem={w'[x]==z[x],z'[x]==(1/10)(1-((w[x])^2))z[x]-w[x],w[0]==1,z[0]==0}

SolnVanderPol2=NDSolve[VanderPolsystem,{w,z},{x,25}, WorkingPrecision-> 10, MaxStepSize->0.005]

ParametricPlot[{w[x],z[x]}/.SolnVanderPol2, {x,0,25}, PlotStyle-> {Blue, Thick}]

Chapter 9: Some Basic Mathematical Programming

In this final chapter, we will discuss some basic mathematical programming concepts such as the For loop command, the use of commands such as While and Do, and the general concept of a procedure in Mathematica which is a series of commands with semicolons between each step of the procedure.  A procedure's steps will be executed in order from left to right if they appear all on the same line or from first to last if they are on successive lines in the same cell.  In general, by design, the output of a procedure is the result of its final step.

Loops, Simple Procedures and Other Commands

Let's begin by doing a few examples involving For loops such as printing values, doing sums and products, and using the commands Do and If.  The first For loop will print (z^k)+z starting with z = (x^2)as k goes from 1 to 4 where at each step it takes the previous value of z and plugs it into (z^k)+z for the next value of k.  The semicolons in the For loop separate the main parts of the procedure inside the For command, which are first the initialization of the loop variable k, the instructions to be carried out successively as k increases in implementing the loop, and finally the output you desire at each step (value of k) of the loop.

The second For loop will compute and display the sums of the reciprocals of the factorial function.  The results get closer and closer to the number e which is approximately 2.71828.

The two procedures after the two For loops will do an addition of the first 10 positive integers and then a multiplication of the same integers using the Do loop command.

For[k=1;z=(x^2),k≤4,k++,z=(z^k)+z;Print[z]]

For[k=1,k≤10,k++,Print[NSum[(1/(n!)),{n,0,k}, WorkingPrecision-> 10]]]

sum=0;Do[sum=sum+k, {k,1,10}];sum

Sum[k, {k,1,10}]

prod=1;Do[prod=k prod, {k,1,10}];prod

10!

prod=1;
Do[prod=prod k, {k,1,10}];
prod

Creating Functions Using Procedures or Nested Commands

Sometimes the most convenient way to create a somewhat complicated function is to use a procedure for its rule.  We will first look at two simple examples which create the factorial function and the add function using the procedures of the previous section to do it.  Remember to place parentheses around all of the steps in your procedure.  The output of the function will be the result of the last step of your procedure which is the function's rule.  The variable n in both the factorial and add functions below takes on only a positive integer value.

factorial[n_]:= (prod=1;Do[prod=k prod, {k,1,n}];prod)

factorial[10]

add[n_]:=(sum=0;Do[sum=sum+k, {k,1,n}];sum)

add[10]

Now we turn to a more sophisticated example involving piecewise functions and the use of the nested If command.  Mathematica has a Piecewise command to build piecewise functions, but it is not always sufficient or convenient for all purposes.  Let's create a function called quadrant which takes a point in the plane and tells us which quadrant it is in or if it is on an axis and so is not in any quadrant.  Note that && is the Mathematica logical And for joining two conditions so that both must be true for them together to be true.  The logical Or is | | with no space between the two vertical bars.  The Return command used below will return an output as well as terminate the If, While, For or Do loop you might be in.  There is also a command Break which terminates loops.

Be very careful in Mathematica to NOT use the letters D and N as variables in your functions or anywhere else as they are Mathematica command names.  The best thing is to stick to using small letters as your variable names and output for any function or procedure.

quadrant[point_]:=If[point[[1]]>0&&point[[2]]>0,Return[quadrant1],If[point[[1]]<0&&point[[2]]>0,Return[quadrant2],If[point[[1]]<0&&point[[2]]<0,Return[quadrant3],If[point[[1]]>0&&point[[2]]<0,Return[quadrant4],axis]]]]

quadrant[{-3,9}]
quadrant[{0,9}]