Introduction to Numerical Methods in Mathematica

A Big Calculator:

Mathematica can be used to perform arithmetic operations for a variety of different types of object.  Ordinary numbers can be used together with arithmetic and transcendental operations (don’t forget to “shift-enter” after the line):

2+5

3 * 7

Sin[Pi/4]

When a number appears by itself, like “2” or “5”, Mathematica tries to keep the output in symbolic, integer form - so that, for example, when we add fractions, we get back a rational number:

1/2 + 2/3

The symbolic nature of Mathematica is useful in some situations, but when used for actual calculations, the symbolic arithmetic operations are very slow.  In addition, we will be working with numbers that have ~15 digits of accuracy, and do not want Mathematica to interfere with that accuracy.  To force Mathematica to evaluate expressions numerically, you can either wrap the entire statement in the “N[]” command:

N[1/2 + 2/3]

or use a period after a number in the expression:

1./2 + 2/3

This approach will also work for transcendental functions:

Sin[6]

versus

Sin[6.]

You can also multiply and add vectors:

5*{1,2,3} + {4,5,6}

Be careful - vectors also work with dot products and matrix-vector multiplication, but the syntax is slightly different.  For a dot product, use a period between two vectors:

{1,2,3}.{4,5,6}

If you apply an operation like Cosine to a vector, Mathematica applies cosine to each term in the vector:

Cos[{Pi/2, Pi/4, Pi/3}]

Variables

You can assign expression to “variable” names.  If we want to give “x” the value “5.0”, we use:

x = 5.0

now we can use “x” and its value will be replaced by 5.0:

x + 2

x/10

Print["The value of x = ", x]

You can assign any type of object to a variable name - for example, we’ll let “u” and “v” be two vectors (of the same length, of course):

u = {1,2,3};
v = {4,5,6};

Notice how we combined two assignments in one entry - the semicolon at the end of the line suppresses output.  We can use the variables to execute arithmetic operations as above:

u+ v

u.v

We can also make matrices - these are just nested vectors from Mathematica’s point of view:

A = N[{{1,2,3},{4,5,6},{7,4,5}}];

You can view the matrix in familiar form using the “MatrixForm” command:

MatrixForm[A]

A dot is used for matrix multiplication:

A.u

and you can also use linear algebra operations on the matrix:

Det[A]

Eigenvalues[A]

Variables retain their value once they have been set.  You can tell if a variable name has not been set by entering it - at first, there is no value for “y”:

y

no value currently, the variable name is returned - if we set y:

y = 10

then we recover its value when y is entered:

y

To clear a variable’s value, you can use the “Clear” command,

Clear[y]

now y is again undefined:

y

A shorthand for “Clear” is “y = .”.

Tables

A table is a vector of entries (the entries themselves could be vectors) - you can make them directly by assigning specific things to a variable:

hodgepodge = {1.,"hello", Sin[Pi/4],{a,b}}

you reference elements of the table using double brackets:

hodgepodge[[1]]

hodgepodge[[2]]

and you can act on the elements in the usual way:

hodgepodge[[1]] +5.

N[hodgepodge[[3]]]

When tables include tables, you have to use multiple indices to get the element of interest.  For example:

hodgepodge[[4]]

has two elements, and if we want the second one, we would use:

hodgepodge[[4]][[2]]

for which there is a shorthand:

hodgepodge[[4,2]]

When a table has more structure to it, you can use the “Table” command to generate it - suppose we wanted a table of the square of integers from 1 to 10:

Table[j^2,{j,1,10}]

or we might want to sample sine over some interval:

sinetable = Table[Sin[2 Pi x],{x,0,1,.1}]

You can plot the elements of a table using the “ListPlot” command:

ListPlot[sinetable]

but ListPlot uses the entry index (starting at 1) for the x-axis of the plot - if we want the correct x values along the axis, we can provide tuples:

sinetable = Table[{x,Sin[2 Pi x]},{x,0,1,.1}];
ListPlot[sinetable]

The points can be connected using “ListLinePlot”:

ListLinePlot[sinetable]

You can also plot pure functions (although Mathematica will discretize them and plot lines connecting points, as above) using “Plot”:

Plot[Sin[2 Pi x],{x,0,1}]

Control Structures

We will use comparison operations to determine outcomes, our primary tool will be the “If” statement.  An “If” statement performs one action if a certain comparison is “True” and another if that comparison is “False”.  The implementation in Mathematica is “If[test, do if test is true, else do if test is false]”.  As an example:

If[3 < 5, Print["Three is less than five."], Print["Three is not less than five."]]

Another important command for controlling the flow of a program is the “For” loop.  In these, an index counts from some provided starting value to a designated ending value, performing some set of operations at each value of index along the way.  The syntax in Mathematica is 
“For[index = start, while index <= end, index = index+1, operations to perform for each value of index]”  Here is an example of a For loop used to print out a sequence of numbers, the index is “j” and it goes from 1 to 5:

For[j = 1, j ≤ 5, j = j+1,
Print[j];
]

You can change the increment, this time, we skip the even numbers:

z = 8;

For[j = 1, j ≤ 5, j = j+2,
Print[j];
]

For loops are useful when you have to perform a set of commands an integer number of times, and where the commands depend on the value of the indexing integer.  Here is an example in which we sum up the numbers from one to ten, depositing the sum in the varialbe “sumup”:

sumup = 0;
For[j = 1, j ≤ 10, j = j+1,
sumup = sumup + j;
];
Print[sumup]

Functions

Mathematica has a stack-like system of definition - once you have entered an expression, it remains “live” until cleared or until you quit Mathematica.  Functions can be defined and used, allowing you to extend the built-in functionality of the program.

We will use the “Module” command to box up series of operations that need to be performed again and again.  Module takes two pieces of input - first the “local” variables that will be used in the function, and then the set of steps to execute (in order).  You name the function, and specify its input parameters via assignment - here is an example of a function that takes in an argument and returns the square of it:

squareit[argone_] := Module[{x},
x = argone*argone;
Return[x];
]

The name of the function is “squareit”, and it takes the argument “argone” (you can name this anything you like - the underscore after the name identifies this as an argument of the function).  Then the “Module” command is used with a local variable “x” and two operations - the first assigns x the value argone*argone, squaring the input, and the second operation returns x to the outside world.  Here is an example of the function in action:

y = squareit[6.]

y

The notion of a “local” variable is important - suppose someone wanted to use the “squareit” function we just defined, but already has a variable “x” alive in the notebook - should squareit be allowed to change the user’s value for x?  No, the function squareit should stand alone and not interfere with the user’s other data - so we make a local variable “x” for use inside squareit (it is defined only inside the function, and preempts any “global” value “x” may have).

Here is a function that takes two vectors and adds them together, returning the result:

addvecs[veca_, vecb_] := Module[{retval, k},
retval = Table[0.,{j,1,Length[veca]}];
For[k = 1, k ≤ Length[veca], k = k+1,
retval[[k]] = veca[[k]] + vecb[[k]];
];
Return[retval];
]

This function takes in two arguments, “veca” and “vecb”, then defines a local variable “retval” to be a table of zeroes of length “Length[veca]” (the length of the vector “veca”) - there is then a “For” loop, that runs k (a local variable, to be safe) from 1 to the length of the vector “veca” in steps of one, and for each value of k, puts the sum of the k^th elements of veca and vecb in the k^th slot of retval.  After the “For” loop has finished, retval is “Returned”.  Here is an example:

x = {1,2,3};
y = {4,5,6};
addvecs[x,y]

There is a problem with our “addvecs” function, though - what if “veca” and “vecb” have different lengths?  There is no sensible addition function in that case, and we’d like to check before blindly adding the vectors together.  We can do that with an “If” statement:

addvecs[veca_, vecb_] := Module[{retval, k},
If[Length[veca] ≠ Length[vecb], 
Print["The two vectors must be the same length"];
Return["error"];
];
retval = Table[0.,{j,1,Length[veca]}];
For[k = 1, k ≤ Length[veca], k = k+1,
retval[[k]] = veca[[k]] + vecb[[k]];
];
Return[retval];
]

Now, if the two vectors are not of the same length (using the “not equal” to expression, obtained via “!=”), an error message is printed, and the function returns the string “error” (not much chance that could be confused with a vector sum!).

x = {1,2,3,4};
y = {4,5,6};
addvecs[x,y]

Exercises:

Here are some things to try out - use the above information as a resource.

1. Mathematica has a built-in help area - you can type “?” and then the name of a command to get more information.  Use “?For”, “?If” and “?While” to learn about the structure of the For, If and While commands.

2. Find the series form of the sine function (you can look it up, or develop it directly from the Taylor expansion).  Write a function that takes in an angle “theta” and a number “Nterms”, and returns the value of “sine(theta)” from the series expansion truncated at Nterms.  Set theta = 1., and run your function for Nterms = 1, 2, 3, 4, 5.  Make a plot of Nterms versus the value your function returns at that truncation level (for theta = 1.).  How close is your truncated value at Nterms =5 to the “true” value determined by Mathematica using Sin[1.] (subtract the two)?

3. Look up the “FindFit” function - this function can be used to fit data to a model equation - you provide the data (a set of {x, f(x)} pairs, in a table, and the model equation with some parameters in it, then “FindFit” will find the “best” values of those parameters for the provided data.  Make some example data, like

data= Table[{x, 5 Exp[-2 x^2]},{x,0,1,.2}];

and use FindFit with model function “a Exp[b x^2]” with this data to find the values of the parameters {a,b}.

4.  Write a function that takes two vectors and returns the dot product (there is a built-in version of this function, but make your own).  Make sure to check for valid vector lengths.  Check it with vectors

a = {1.,2.,3.,4.,5.};
b= {2.,4.,6.,8.,10.};

and

a = {1.,2.,3.,4.,5.};
b= {2.,4.,6.,8.,10.,11.};