ECM6 Computational Methods :

Slide 1 of 12

Lecture 3
Introduction to Algorithms
Brian G. Higgins
Department of Chemical Engineering & Materials Science
University of California, Davis
April 2014, Hanoi, Vietnam

ECM6Lecture3Vietnam_2014.nb

Topics for Lecture 3

In this lecture we give a brief over view of the following topics

1. Concept of an algorithm
2. Example involving summations
3. Algorithm for approximating an integral (trapezoidal rule) using Table
4. Algorithm for approximating integral using For loop

ECM6Lecture3Vietnam_2014.nb

Introduction to Algorithms
An algorithm is a step-by-step procedure that produces a solution in a finite number of steps. The basic
components of an algorithm are
(i) List of the input parameters that must be supplied
(ii) List of operations that must be performed and in what order
(iii) Output of the algorithm.
When the output is a numerical result we call the procedure a numerical method. Numerical methods by
their nature do not produce exact solutions. There is always some error associated with the method.
Thus one goal of a numerical method is to assess how close the "solution" is to the exact solution.
Matters are complicated because the exact solution is often not known. If it were there would be no
reason to use a numerical method.
Here are some simple example of an algorithm

Example
Evaluate the mean and standard deviation of a given set of data (measurements). Suppose the data is
in the form of a sequence
data = 8x1 , x2 , x3 , , xn <
The quantities that we want to evaluate are given by the following formulas
mean :

Xx\ =

1
n

xi
i=1

Standard Deviation : s =

1
n Hn - 1L

x2i
i=1

- xi
i=1

ECM6Lecture3Vietnam_2014.nb

Algorithm
A word description of the steps in an algorithm
Given Inputs: xi
Step 1:
Calculate n (number of data points)
Step 2:
Initialize variables for algorithm:
xsum=0, x2sum=0
Step 3:
for i from 1 to n
add xi to xsum
add xi2 to x2sum
end
Step 4:
calculate xbar= xsum/n
Step 5:
calculate s=Sqrt(n(x2sum)- HxsumL2 )/(n(n-1)
Output:
xbar and s

ECM6Lecture3Vietnam_2014.nb

Mathematica Code for Algorithm

This algorithm can be implemented in many types of programming languages. Here we will show how it
is done in Mathematica. Let us suppose our data is given by.
x={13.247,12.911,7.821,9.283,12.874,12.880,2.472,10.782,6.605,13.506,12.686,10.317,6.387,
5.886,9.815}

Then the Mathematica code is

Input: Define the list of data

In[1]:=

x = 813.247, 12.911, 7.821, 9.283, 12.874, 12.880,

2.472, 10.782, 6.605, 13.506, 12.686, 10.317, 6.387, 5.886, 9.815<;

Step 1: Determine n numebr of data points

In[2]:=
Out[2]=

n = Length@xD
15

Step 2: Intialize variables

In[3]:=

xsum = 0;
x2sum = 0;

Step 3: Compute xsum and xs2sum using a For loop

Information about For loop
In[5]:=

? For
For@start, test, incr, bodyD executes start, then repeatedly evaluates body and incr until test fails to give True.

In[6]:=

ForAi = 1, i < n + 1, i ++, xsum = xsum + x@@iDD; x2sum = x2sum + x@@iDD2 E

Step 4: Calculate mean : xbar

In[7]:=
Out[7]=

xbar = xsum n
9.83147

Step 5: Calculate standard deviation: s

In[8]:=

Out[8]=

s=

1
n Hn - 1L

0.232434

n x2sum - xsum2

ECM6Lecture3Vietnam_2014.nb

Output: xbar, s
In[9]:=

Out[9]=

Out[10]=

xbar
s
9.83147
0.232434

Final Comments
This is not the only way to do the calculation in Mathematica. Here are some alternative methods that
make use of Mathematica's built-in functions

ECM6Lecture3Vietnam_2014.nb

Development of Algorithm without For Loop

In this example we use Mathematicas Sum function rather than use a For loop:

Example 1:
In[11]:=

x = 813.247, 12.911, 7.821, 9.283, 12.874, 12.880, 2.472,

10.782, 6.605, 13.506, 12.686, 10.317, 6.387, 5.886, 9.815<;
n = Length@xD;
xsum = Sum@x@@iDD, 8i, 1, n<D;
x2sum = SumAx@@iDD2 , 8i, 1, n<E;
xbar = xsum n
1
s=
n x2sum - xsum2
n Hn - 1L

Out[15]=

9.83147

Out[16]=

0.232434

Example 2:
Here is another method that makes use of functional programming concepts. In this example we use
the Apply function.
In[17]:=

? Apply
Apply@ f , exprD or f expr replaces the head of expr by f .
Apply@ f , expr, 81<D or f expr replaces heads at level 1 of expr by f .
Apply@ f , expr, levelspecD replaces heads in parts of expr specified by levelspec.

In[18]:=

x = 813.247, 12.911, 7.821, 9.283, 12.874, 12.880, 2.472,

10.782, 6.605, 13.506, 12.686, 10.317, 6.387, 5.886, 9.815<;
n = Length@xD;
xbar = Apply@Plus, xD n
1
s=
n ApplyAPlus, x2 E - HApply@Plus, xDL2
n Hn - 1L

Out[20]=

9.83147

Out[21]=

0.232434

Example 3:
Here is another method that makes use of functional programming concepts. In this example we use
the Fold function.

ECM6Lecture3Vietnam_2014.nb

In[22]:=

? Fold
Fold@ f , x, listD gives the last element of FoldList@ f , x, listD.

In[23]:=

x = 813.247, 12.911, 7.821, 9.283, 12.874, 12.880, 2.472,

10.782, 6.605, 13.506, 12.686, 10.317, 6.387, 5.886, 9.815<;
n = Length@xD;
xbar = Fold@Plus, 0, xD n
1
s=
n FoldAPlus, 0, x2 E - HFold@Plus, 0, xDL2
n Hn - 1L

Out[25]=

9.83147

Out[26]=

0.232434

ECM6Lecture3Vietnam_2014.nb

Algorithm for Evaluating an Integral

In this example we want to approximate a definite integral
b

Int = f HxL x
a

Our algorithm will be based on

Defining a partition for the x- axis

a = x0 < x1 < x2 < < xn-1 < xn = b
where xi = a = i h for each i and h = Hb - aL n.

Approximating the integral with the trapezoidal rule

b

Tn = f HxL x
a

h
2

n-1

Bf HaL + 2 f Hxi L + f HbLF

i=1

Algorithm
Given Inputs: Limits of integration a,b, integrand f(x)
and number of sub intervals
Step 1:
compute h; set sum=0
Step 2:
Step 3:
Step 4:
Output:

for i from 1 to n-1

add fHa + i hL to sum
end
multiply sum by 2
add f(a) and f(b) to sum
(h/2) sum

10

ECM6Lecture3Vietnam_2014.nb

Mathematica Code
Let us evaluate the following definite integral
5

Cos HxL x
0

We will partition the interval into 100 equally spaced points.

Here is a plot of the integrand
In[27]:=

Remove@x, a, b, h, n, i, mysum, f, T, eD

In[28]:=

Plot@Cos@xD, 8x, 0, 5<, PlotStyle Blue, Frame True, FrameLabel 8"x", "fHxL"<D
1.0

fHxL

0.5

0.0

Out[28]=

-0.5

-1.0
0

Mathematica Code for Evaluating the Integral

In[29]:=

Out[32]=

a = 0; b = 5; f@x_D := Cos@xD; n = 100;

h = Hb - aL n;
mysum = Sum@f@a + i hD, 8i, 1, n - 1<D;
h
T@100D =
Hf@aD + 2 mysum + f@bDL N
2
- 0.958724

Let us evaluate the integral using Mathematica's NIntegrate function

In[33]:=
Out[33]=

NIntegrate@Cos@xD, 8x, 0, 5<D

- 0.958924

Assessing Accuracy of Algorithm

We can asses the accuracy of the trapezoidal approximation by evaluating the error
b

e@nD = T@nD - f HxL x

a

ECM6Lecture3Vietnam_2014.nb

Let us see how the solution varies with the value of n. We will define a function for the error
In[34]:=

Remove@TD;
a = 0; b = 5; f@x_D := Cos@xD;
Hb - aL
T@n_D :=
Hf@aD + 2 Sum@f@a + i Hb - aL nD, 8i, 1, n - 1<D + f@bDL N
2n
b

In[37]:=

e@n_D := NB T@nD - Cos@xD xF

a

Let us examine what the error is when n=100

In[38]:=
Out[38]=

e@100D
0.000199784

Displaying the Data with TableForm

We can display the data as a function of n, using the function TableForm
In[39]:=

? TableForm
TableForm@listD prints with the elements of list arranged in an array of rectangular cells.

In[40]:=

TableForm@Table@8n, T@nD, Abs@e@nDD, Abs@e@nDD Abs@e@n 2DD<, 8n, 2, 32, 2<D,

TableHeadings -> 8None, 8"n", "Tn ", "en ", "en en2 "<<,
TableAlignments -> CenterD

Out[40]//TableForm=

Tn
- 0.398281
- 0.830687
- 0.902778
- 0.927504
- 0.938863
- 0.945011
- 0.94871
- 0.951108
- 0.95275
- 0.953925
- 0.954793
- 0.955453
- 0.955967
- 0.956375
- 0.956704
- 0.956973

n
2
4
6
8
10
12
14
16
18
20
22
24
26
28
30
32

en
0.560643
0.128237
0.0561464
0.0314201
0.0200613
0.0139136
0.0102144
0.00781648
0.00617387
0.00499961
0.00413116
0.00347084
0.00295709
0.00254952
0.00222076
0.00195173

en en2
0.134509
0.228733
0.240992
0.245015
0.246834
0.24781
0.248395
0.248773
0.249032
0.249216
0.249353
0.249456
0.249537
0.249601
0.249652
0.249694

Displaying the Data with Grid

In[41]:=

? Grid
Grid@88expr11 , expr12 , <, 8expr21 , expr22 , <, <D is
an object that formats with the exprij arranged in a two-dimensional grid.

11

12

ECM6Lecture3Vietnam_2014.nb

In[42]:=

In[43]:=

Out[43]=

mytable = Table@8n, T@nD, Abs@e@nDD, Abs@e@nDD Abs@e@n 2DD<, 8n, 2, 32, 2<D;
mytable2 = Prepend@mytable, 8"n", "Tn ", "en ", "en en2 "<D;
Grid@mytable2, Frame AllD
n
2
4
6
8
10
12
14
16
18
20
22
24
26
28
30
32

Tn
- 0.398281
- 0.830687
- 0.902778
- 0.927504
- 0.938863
- 0.945011
- 0.94871
- 0.951108
- 0.95275
- 0.953925
- 0.954793
- 0.955453
- 0.955967
- 0.956375
- 0.956704
- 0.956973

en
en en2
0.560643
0.134509
0.128237
0.228733
0.0561464
0.240992
0.0314201
0.245015
0.0200613
0.246834
0.0139136
0.24781
0.0102144
0.248395
0.00781648 0.248773
0.00617387 0.249032
0.00499961 0.249216
0.00413116 0.249353
0.00347084 0.249456
0.00295709 0.249537
0.00254952 0.249601
0.00222076 0.249652
0.00195173 0.249694

The function Grid takes on many options for embellishing the output of the table.
In[44]:=
Out[44]=

Options@GridD
8Alignment 8Center, Baseline<, AllowedDimensions Automatic,
AllowScriptLevelChange True, AutoDelete False, Background None,
BaselinePosition Automatic, BaseStyle 8<, DefaultBaseStyle Grid,
DefaultElement , DeleteWithContents True, Dividers None,
Editable Automatic, Frame None, FrameStyle Automatic, ItemSize Automatic,
ItemStyle None, Selectable Automatic, Spacings Automatic<

Here is an example of what you can do:

ECM6Lecture3Vietnam_2014.nb

In[45]:=

Grid@mytable2, Frame Darker@Gray, .6D, ItemStyle 14, Spacings 8Automatic, .8<,

Dividers 88Darker@Gray, .6D, 8Lighter@Gray, .5D<, Darker@Gray, .6D<,
8Darker@Gray, .6D, Darker@Gray, .6D, 8False<, Darker@Gray, .6D<<,
Alignment 88Left, Right, 8Left<<<, Background
8None, 8Lighter@Yellow, .9D, 8White, Lighter@Blend@8Blue, Green<D, .8D<<<D

en en2

-0.398281 0.560643

0.134509

-0.830687 0.128237

0.228733

-0.902778 0.0561464

0.240992

-0.927504 0.0314201

0.245015

10 -0.938863 0.0200613

0.246834

12 -0.945011 0.0139136

0.24781

14
Out[45]=

Tn en

-0.94871 0.0102144

0.248395

16 -0.951108 0.00781648 0.248773

18

-0.95275 0.00617387 0.249032

20 -0.953925 0.00499961 0.249216

22 -0.954793 0.00413116 0.249353
24 -0.955453 0.00347084 0.249456
26 -0.955967 0.00295709 0.249537
28 -0.956375 0.00254952 0.249601
30 -0.956704 0.00222076 0.249652
32 -0.956973 0.00195173 0.249694

13

14

ECM6Lecture3Vietnam_2014.nb

Algorithm for Iteration

In this example we will discuss the concept of a repeated substitution algorithm for the function
g HxL =

1
a

Ha - 1L x +

b
x

The essential idea is as follows. We are given a function g(x), and we want to iterate based on the
following algorithm
xn+1 = g Hxn L
starting with the value x0 .
We say the function g(x) has a fixed point if we can find a value x = xp such that
xp = g Ixp M
where xp is called the fixed point of g(x).
In[46]:=

Remove@g, x, a, i, x0, xp, g, resultListD

Remove::remal : Symbol Removed@gD already removed.

ECM6Lecture3Vietnam_2014.nb

15

Exploring the Function with Mathematica

Consider the following function
In[47]:=

g@a_, b_, x_D :=

1
a

Ha - 1L x +

b
x

This function is used to compute the square root of a real number a. The parameter a in the recursive
algorithm influences the order of convergence. When xn+1 = xn we say that the iteration has converged
to a stable fixed point.
Let us determine the fixed points of this function for b=78.4, a=5

Plot of function g(x)

In[48]:=

Plot@g@2, 78.4, xD, 8x, - 15, 15<, Frame True, PlotStyle 8Thick, Blue<,
FrameLabel 8Style@"x", 16D, Style@"gHx<", 16D<, AspectRatio 1D
40

gHx<

20

Out[48]=

-20

-40
-15

-10

-5

10

15

16

ECM6Lecture3Vietnam_2014.nb

Plot of function g(x)=x

In[49]:=

Plot@8g@2, 78.4, xD, x<, 8x, - 15, 15<,

Frame True, PlotStyle 88Thick, Blue<, 8Thick, Red<<,
FrameLabel 8Style@"x", 16D, Style@"gHx<", 16D<, Exclusions 80<, AspectRatio 1D
40

gHx<

20

Out[49]=

-20

-40
-15

-10

-5

10

15

x
Let us determine the fixed points of this function for b=78.4, a=5 with
In[50]:=

Solve@g@a, 78.4, xD x, xD
Solve::ratnz : Solve was unable to solve the system with inexact coefficients. The
answer was obtained by solving a corresponding exact system and numericizing the result.

Out[50]=

88x - 8.85438<, 8x 8.85438<<

Note in this example the fixed points do not depend on the parameter a. Further
Thus the fixed points are the values

78.4 = 8.85438.

b . How do we program Mathematica to produce the iterates?

ECM6Lecture3Vietnam_2014.nb

17

Fixed Point Algorithm

How do we program Mathematica to produce the iterates? We can use a For loop. The syntax for the
For loop is
For[start, test, incr, body]

The For loop executes start, then repeatedly evaluates body and incr until test fails to give True.
Here is an example where the arguments of the For function are defined as
start block:
test block:
incr block:
body block:

i = 1; x0 = 0.2
i<20
i++
xp=g[1.5,x0 ];Print[xp];x0=xp

Note that the final value returned by For is Null. For this reason we use the Print function to determine
the values of x during the iteration.

Mathematica Code
In[51]:=

For@i = 1; x0 = 0.2, i < 15, i ++, xp = [email protected], 78.4, x0 D; Print@xpD; x0 = xpD

261.4
87.3333
29.7096
11.6624
8.3691
9.03489
8.79661
8.87389
8.8479
8.85654
8.85366
8.85462
8.8543
8.8544

Rather than print the results of the iteration at each step we can store the intermediate results in a list
In[52]:=

For@resultList = 8<; i = 1; x0 = 0.2, i < 25, i ++,

xp = [email protected], 78.4, x0 D; AppendTo@resultList, x0D; x0 = xpD

18

ECM6Lecture3Vietnam_2014.nb

In[53]:=
Out[53]=

resultList
80.2, 261.4, 87.3333, 29.7096, 11.6624, 8.3691, 9.03489, 8.79661,
8.87389, 8.8479, 8.85654, 8.85366, 8.85462, 8.8543, 8.8544, 8.85437,
8.85438, 8.85438, 8.85438, 8.85438, 8.85438, 8.85438, 8.85438, 8.85438<

It is clear the iterates in this example have stabilized (converged) to a value of 8.85438 =

78.4 .

If we give a negative value for the initial guess, we converge to the other root!
In[54]:=

Out[55]=

For@resultList = 8<; i = 1; x0 = - 0.2, i < 25,

i ++, xp = [email protected], 78.4, x0 D; AppendTo@resultList, x0D; x0 = xpD;
resultList
8- 0.2, - 261.4, - 87.3333, - 29.7096, - 11.6624, - 8.3691, - 9.03489, - 8.79661, - 8.87389,
- 8.8479, - 8.85654, - 8.85366, - 8.85462, - 8.8543, - 8.8544, - 8.85437, - 8.85438,
- 8.85438, - 8.85438, - 8.85438, - 8.85438, - 8.85438, - 8.85438, - 8.85438<

ECM6Lecture3Vietnam_2014.nb

Final Comments
Questions we would like to know are:
Does the algorithm always converge?
(What is the order of convergence, i.e what is the change in the error at each iteration?
These topics will be more fully discussed in future Lectures.

19