Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Aug 2016
Lecture topic: Programming in Mathematica Introduction

Mathematica's programming syntax

We will not be doing a tremendous amount of programming in this class, but you will want to know how to set up simple programming structures in Mathematica to help you with some of your assignments. More complicated programming structures we will look at in class.

Looping

Do, While and For do not print any output; they simply do what they are asked to do. After the calculation you will need to output whatever you are interested in. Table is very similar to Do, but it forms a list of the output. In many situations Table is all you will need. Other commands that may prove useful are Sum and Product. Notice the syntax of the semicolon after each line in the loop.

Do Loops

Clear[a,b,c]

a=0;
b=0;
c=0;
Do[
a=a+n;
b=a+b;
c=b^2,
{n,1,4}]

{a,b,c}

While Loops

Clear[a,b,c,n]

a=b=c=0; n=1;
While[n≤4,
a=a+n;
b=a+b;
c=b^2;
n=n+1]

{a,b,c}

For Loops

Clear[a,b,c,n]

For[a=b=c=0; n=1,
n≤4,
n=n+1,
a=a+n;
b=a+b;
c=b^2]

{a,b,c}

Table Loops

Clear[a,b,c]

a=b=c=0;
Table[
a=a+n;
b=a+b;
c=b^2,
{n,1,4}]

{a,b,c}

Branching

Branching in Mathematica can be accomplished via the commands If, Which, or Switch. The most common one you will use will probably be IF.

If

f[x_] = If[x≤0,x^2,x];
Plot[f[x],{x,-3,3},AxesLabel->{"x","y"},PlotStyle->{Blue}]

Flow Control

You can control the flow of programming in Mathematica using the commands Continue, Break, Return, and Goto. Here is a basic example which shows how these work.

f[n_] := 
Do[
Which[
n==1,Print["n=1 ","k= ",k]; Continue[],
n==2,Print["n=2 ","k= ",k];Break[],
n==3,Return["That's All Folks"],
True,Goto[jump];
Label[jump];Return[Print["n= ",n]]]
,{k,1,3}]

f[1]

f[2]

f[3]

f[4]

Selecting Elements from a List

To select elements from a list, use the command Part.

Clear[a,b,c]
data = Table[{a[n],b[n],a[n]*b[n]^2},{n,0,5}];
TableForm[data,TableHeadings->{None,{"1","2","3"}}]

data2 = Table[{Part[data,n+1,2],Part[data,n+1,3]},{n,0,5}];
TableForm[data2,TableHeadings->{None,{"2","3"}}]

A shorthand does exist for Part, and it is done by typing <esc> [[ <esc> and then the part of the item you want followed by <esc> ]] <esc>. Here is an example, try it yourself!

The full list of data:

data

The third element in data:

data[[3]]

The second element in the third element in data:

data[[3,2]]

Numerical Output

If you simply leave Mathematica to its defaults, it will output numbers to 6 digit precision, or significant digits precision. For example:

N[Pi]
N[Pi,10]
N[3/4]
N[3/4,8]

If we want to include more digits in the output, we should use the command SetPrecision when we output the number. We would still like Mathematica to work with high digits of precision internally, but we want to be able to see as many digits as we like when needed.

SetPrecision[Pi,9]
SetPrecision[3/4,1]
SetPrecision[3/4,5]

The machine precision of Mathematica can be determined:

$MachinePrecision

If we use the symbolic representation of numbers in Mathematica to do arithmetic (variable precision arithmetic) we get no round off error, but our calculations will be slow.

Pi-Pi
Timing[Sum[(1-24)^n ,{n,0,2003}]]

If we use the hardware to do our arithmetic (fixed precision arithmetic), we get results much faster, but we introduce roundoff error:

N[Pi,25]-N[Pi,25]
Timing[Sum[(1.0-24.0)^n ,{n,0,2003}]]

This may not seem like a big time difference for the accuracy you get out at the end, and it really isn't for this simple case. However, when you inadvertantly use arbitrary precision numbers inside a looping structure, you can end up with code which takes days to run instead of minutes. The lesson: know your expected error and work within it!

Putting It All Together: Euler's Method

For solutions to assignment problems, I am not expecting you to do high level programming in Mathematica. For example, here are some ways of implementing Euler's method to solve an initial value problem. I would expect you to be able to create the recursion relation method, or the While loop method, but not the module method (although you probably will start using Modules at some point during the semester). You should feel free to explore some of the more advanced syntax of Mathematica if it pleases you!

Here we want to solve the initial value problem y'(t) = f(y,t), y(α)=β for y(t) in a given range β≤t≤γ. Euler's method generates approximations to the solution y(t) using the following relations, once you have picked a value for the step size h:

(Subscript[t, n])=(Subscript[t, n-1])+h
(Subscript[y, n])=(Subscript[y, n-1])+h *f((Subscript[y, n-1]),(Subscript[t, n-1])), n=1,2,3,...
(Subscript[t, 0])=α,(Subscript[y, 0])=y(α)=β.

Let f(y,t)=y-(t^2)+1, y(0)=1/2 and 0≤t≤2. I've included some extra tables just to show you what Mathematica can do.

Euler's Method: While Loop

Clear[f,y,t,a,b,α,data,num]

(*initialization of Euler's method*)
f[t_,y_]=y-(t^2)+1;
a=0.;
b=2.;
α=0.5;
num=10;
h=((b-a)/num);
yi=α;
ti=a;

(*result stored in data*)
data={{ti,yi}};
n=1;

(*Euler's method*)
While[n≤num,{ti,yi}=N[{t+h,y+h f[t,y]}/.{t->ti,y->yi}];data=Append[data,{ti,yi}];n=n+1]

SetPrecision[data,9]

ListPlot[data,PlotRange->{{0,2},{0,5}},Joined->True,Frame->True,FrameLabel->{"t","y"},PlotStyle->{Blue,Thick,Dashed}]

Euler's Method: Recursion Relations

Clear[f,y,t,a,b,α,data]

(*initializations*)
f[t_,y_]=y-(t^2)+1;
a=0.;
b=2.;
α=0.5;
num=10;
h=((b-a)/num);

(*Euler's emthod*)
t[0]=a;
y[0]=α;
y[n_]:=y[n]=y[n-1]+h f[t[n-1],y[n-1]]
t[n_]:=t[n]=a+n h

(*results stored in data*)
data=Table[{n,t[n],N[y[n],20]},{n,0,num}];
SetPrecision[data,9]
data2=Table[{data[[n+1,2]],data[[n+1,3]]},{n,0,num}];
SetPrecision[data2,9]

phase=Table[{data[[n+1,2]],data[[n+1,3]]},{n,0,num}];

ListPlot[phase,PlotRange->{{0,2},{0,5}},Joined->True,Frame->True,FrameLabel->{"t","y"},PlotStyle->{Red,Thick,Dashed}]

Euler's Method: Modules

Modules allow you to create executable blocks of code, with variables that are internal to the module.

You should add some comments which explain what the variables are, and which variables are local. This is just good programming practice, but I often forget to go back and do it! This is something I need to be more conscientous at doing. The reason we do this is because often code you write will be used in new situations, by you or your collaborators.

Clear[f,y,t,a,b,α,data]

euler[f_,t_,y_,t0_,y0_,h_,n_]:=
Module[{ti=t0,yi=y0,p=8},
(*Euler's method module
Inputs 
f rate function in IVP
y dependent variable
t independent variable
t0 initial t value
y0 initial y value
h step size
n number of steps
p is precision (local variable)
ti is local t variable
yi is local y variable

Output
data set approximation to solution to IVP
*)

Prepend[Table[{ti,yi}=SetPrecision[{t+h,y+h*f} /.{t->ti,y->yi},p],{k,1,n}],{t0,y0}]
]

Here is the output from executing the module:

approx=euler[ y-t^2+1,t,y,0.0,0.5,0.2,10]

Here the output plotted.

ListPlot[approx,PlotRange->{{0,2},{0,5}},Joined->True,Frame->True,FrameLabel->{"t","y"},PlotStyle->{Blue,Thick,Dashed},BaseStyle->20]