## Laboratory 4

This laboratory will investigate two topics:

1. Iteration, the fixed-point theorem and convergence (MATH 211); and
2. Interpolating polynomials and polynomials in Matlab (MATH 215).

# 4.1 Fixed-point Iteration

Given any real-valued function of a real variable, , if we start with an initial point and then iterate

for , then one of two possible results may occur:

• Either the sequence will converge to a point , or
• It will not converge.

If it converges, it converges to a solution to a solution of the equation

.

In order to quickly demonstrate this, consider the cosine function: If we plot and , we note that they cross somewhere around . If we start with and iterate, we get

Notice that it does seem to be approaching the point which is the point at which .

## 4.1.1 Implementation in Matlab

Save the function

function [x] = fp( f, x1 )
x = x1;

for k = 1:100
x_prev = x;

x = f(x);

if abs( x - x_prev ) < 1e-10
return;
end
end

throw( MException( 'MATLAB:non_convergence', ...
'The fixed-point iteration did not converge' ) );
end


This function introduces three new new structures:

#### 4.1.1.1 If Statements

The if statement in Matlab is similar to the C/C++/Java/C# if statement: a Boolean-valued predicate is tested and if the predicate returns true, the body of the if statement is executed. In the above example, the predicate was a test if .

A more general structure is the if-else-end which is similar to the if (...) {...} else {...} format:

	if predicate
% body of code to be executed if the predicate is true
else
% body of code to be executed if the predicate is false
end


It is also possible to have more than one predicate:

	if predicate #1
% body of code to be executed if the predicate #1 is true
elseif predicate #2
% body of code to be executed if predicate #1 is false but predicate #2 is true
elseif predicate #3
% body of code to be executed if predicates #1 and #2 are false but predicate #3 is true
else
% body of code to be executed if all three predicates are false
end


You can also type help if at the Matlab command prompt and you can also view the Matlab help page.

#### 4.1.1.2 The Line Continuation Operator ...

Matlab assumes that one line contains one command. If a command is to be spread over multiple lines, the programmer should do one of two things:

1. Determine if the command is too complex, in which case, it can be broken down into a sequence of commands, or
2. Before the end of any line other than the last line in the command, append three dots ... to indicate to Matlab that the command continues on the next line.

#### 4.1.1.3 The throw Function and Exceptions

Throwing an exception in Matlab is similar to throwing exceptions in C++/Java/C#. The format of the exception object is two strings MException( identifier, message ).

• The first argument, the identifier, is of the form MATLAB:string which identifies the exception. In this course, the identifier will always be provided.
• The second argument is a string which is to be printed to the screen if the exception is not caught.

For example, the iteration for the function does not converge with the initial point :

fp( @tan, 1 )
??? Error using ==> fp at 12
The fixed-point iteration did not converge


For the purposes of this course, the message will always be a fixed string; however, it is possible to use a format similar to printf for those of you who have already learned native C (see Matlab's help page for the exception class).

## 4.1.2 Examples of Fixed-Point Iterations

If we iterate the cosine function starting with , it does converge to the solution of shown in Figure 1.

Figure 1. The plot of the cosine function and .

The Matlab code is here:

format long
fp( @cos, 1 )
ans =
0.739085133245110

cos( ans )
ans =
0.739085133194986


The actual value of the root, to 20 decimal digits, is .

Two examples of non-convergence include the tangent function (as shown above) but also the sine function:

fp( @sin, 1 )
??? Error using ==> fp at 12
The fixed-point iteration did not converge


If you were to remove the semicolon from the line

		x = f(x)


each point would be printed to the screen. You would see the sequence of points

fp( @sin, 1 )
0.841470984807897
0.745624141665558
0.678430477360740
0.627571832049159
0.587180996573431
0.554016390755630
0.526107075502842
0.502170676268555
0.481329355262346
0.462957898537812
0.446596593386982
0.431898433265822
0.418595660998210
0.406477764984069
0.395376546988835
0.385155713096526
0.375703446847715
0.366927001062791
0.358748688107810
0.351102858963179
0.343933594328648
0.337192916936002
0.330839391070358
0.324837013640270
0.319154327474343
0.313763705915325
0.308640770822766
0.303763915468900
0.299113910636710
0.294673577256255
0.290427512659243
0.286361860348455
0.282464115317791
0.278722958597813
0.275128115968017
0.271670236763068
0.268340789473641
0.265131971453282
0.262036630528281
0.259048196695839
0.256160622408316
0.253368330194057
0.250666166570835
0.248049361375991
0.245513491775248
0.243054450326041
0.240668416565461
0.238351831671369
0.236101375810769
0.233913947844445
0.231786647103119
0.229716756989405
0.227701730192909
0.225739175333946
0.223826844875308
0.221962624161997
0.220144521466446
0.218370658931835
0.216639264319163
0.214948663474997
0.213297273446584
0.211683596179494
0.210106212740356
0.208563778013676
0.207055015827396
0.205578714466746
0.204133722540338
0.202718945166210
0.201333340448950
0.199975916221949
0.198645727031511
0.197341871341849
0.196063488942078
0.194809758538159
0.193579895514375
0.192373149850419
0.191188804181423
0.190026171989499
0.188884595916346
0.187763446187463
0.186662119139339
0.185580035841746
0.184516640807959
0.183471400786333
0.182443803627253
0.181433357219929
0.180439588494031
0.179462042481508
0.178500281434370
0.177553883994511
0.176622444411982
0.175705571808417
0.174802889482536
0.173914034254936
0.173038655849547
0.172176416309364
0.171326989444225
0.170490060308578
0.169665324707324
0.168852488727981

??? Error using ==> fp at 12
The fixed-point iteration did not converge


It is clear that is a solution to the equation and it does look like the sequence above is converging to 0; however, the rate of convergence is very slow.

## 4.1.3 Possible Outcomes Fixed-Point Iterations

Thus, for fixed point iterations, we have three possibilities:

• The iterations converge to a value close to a solution of ,
• The iterations do not converge, and
• The iterations appear to converge, but the convergence is so slow as to be useless for practical purposes.

We will use iteration for most numerical methods in this course, and each iteration is—in a sense—a fixed-point iteration. In each case, we will have a definite structure to the routines. We will start with one or more initial points (depending on the algorithm) and continue iterating until:

• The difference between iterative steps is sufficiently small: we will call this or eps_step, in which case, the last result of the iteration will be taken to be the approximation, or
• We will stop iterating after some fixed number of iterations .

# 4.2 Finding Interpolating Polynomials

In the previous Laboratory, we worked on Euler's method for solving first-order initial value problems. Before we can proceed to finding better algorithms and solving more difficult problems, we must first look at interpolation.

## 4.2.1 The Theory

### 4.2.1.1 Interpolating Two Points

Given two points and , it is possible to find a straight line which passes through the two points so long as . This is shown in Figure 2 where the polynomial passes through the points and .

Figure 2. A line interpolating two points.

Recall from Linear Algebra that to find the straight line which passes through two points, we must satisfy two equations:

which may be written as

.

A quick inspection shows that the determinant is which is zero if and only if .

### 4.2.1.2 Interpolating Three Points

Given three points , and , it seems reasonable, also, that a quadratic function could be passed through the three points. For example, the quadratic function passes through the three points , and , as is shown in Figure 3.

Figure 3. A quadratic function interpolating three points.

Again, from Linear Algebra, to find such a quadratic function, we must satisfy the three equations

which may be written as

.

A quick inspection shows that the determinant is which, again, is zero if and only if , , or .

### 4.2.1.3 Interpolating Points

Given points , is it possible to find a polynomial of degree which passes through those points?

The answer, as you may have guessed, is yes, provided the x-values are all different. The process is similar, a system of linear equations is created and solved. Without proof, we will state that such a system always has a unique solution and therefore there always exists such an interpolating polynomial.

## 4.2.2 Interpolating in Matlab

In order to interpolating the two points, and , we must solve the linear system

.

In Matlab, this is done as follows:

[3.2 1; 5.0 1] \ [4.7 1.8]'

ans =

-1.6111
9.8556


and therefore the interpolating line is .

In order to interpolating the three points, , and , we must solve the linear system

.

[3.2^2 3.2 1; 5.0^2 5.0 1; 9.4^2 9.4 1] \ [4.7 1.8 6.5]'
ans =
0.4321
-5.1547
16.7699



and therefore the interpolating line is .

By this point, you should be aware that this would be very frustrating do on a regular basis, and therefore we will introduce two new features in Matlab: element-wise exponentiation of vectors and matrices and more powerful matrix constructors.

### 4.2.2.1 Element-wise Exponentiation

Normally, if you consider the Matlab statement M^2, you would quickly realize that M must be a square matrix for this to make sense, in which case, M^2 is the same as M*M. Also, it does not simply square the entries of the matrix:

[1 2; 3 4]^2
ans =
7    10
15    22


If you want to simply raise each entry of a matrix to a power, you can do this using the element-wise exponentiation operator .^ as is demonstrated here:

[1 2; 3 4].^2
ans =
1     4
9    16


This works with matrices of an arbitrary size including vectors:

xs = [3.2 5.0 9.4]'
xs =
3.2000
5.0000
9.4000

xs.^2
ans =
10.2400
25.0000
88.3600


### 4.2.2.2 Joining Matrices

Previously, we saw how we can define the entries of a matrix by just listing them using semicolons where appropriate, e.g., [1 2; 3 4]. This notation, however, also allows you to join matrices of the appropriate size:

M = [1 2 3; 4 5 6; 7 8 9];
u = [13 14 15]';
v = [10 11 12];
[M u]
ans =
1     2     3    13
4     5     6    14
7     8     9    15

[M; v]
ans =
1     2     3
4     5     6
7     8     9
10    11    12

[M u; v 16]
ans =
1     2     3    13
4     5     6    14
7     8     9    15
10    11    12    16


### 4.2.2.3 Interpolation with Matlab

Thus, it should be apparent that the following will work:

xs = [1 3 5]';
ys = [2 1 4]';
M = [xs.^2 xs.^1 xs.^0]
M =
1     1     1
9     3     1
25     5     1

M \ ys
ans =
0.5000
-2.5000
4.0000


which gives us the polynomial shown in Figure 3.

With our second example, we have

xs = [3.2 5.0 9.4]';
ys = [4.7 1.8 6.5]';
M = [xs.^2 xs.^1 xs.^0]
M =
10.2400    3.2000    1.0000
25.0000    5.0000    1.0000
88.3600    9.4000    1.0000

M \ ys
ans =
0.4321
-5.1547
16.7699


### 4.2.2.4 Plotting the Points and the Interpolating Quadratic

Suppose we have the solution

xs = [1 3 5]';
ys = [2 1 4]';
M = [xs.^2 xs.^1 xs.^0];
c = M \ ys;    % The polynomial is c(1)*x^2 + c(2)*x + c(3)
c =
0.5000
-2.5000
4.0000


and we wish to plot it.

xr = linspace( 0, 6, 100 );
yr = c(1)*xr.^2 + c(2)*xr + c(3);
plot( xs, ys, 'o' );
hold on
plot( xr, yr, 'r' );


The result is the plot shown in Figure 4.

Figure 4. The Matlab plot.

### 4.2.3 Actual Polynomial Tools in Matlab

The previous examples have shown how we can create the necessary matrices using Matrix constructors. Matlab also provides a similar tool; the Matlab vander function.

xs = [1 3 5]';
M1 = vander( xs )
M1 =
1     1     1
9     3     1
25     5     1

M = [xs.^2 xs.^1 xs.^0]
M2 =
1     1     1
9     3     1
25     5     1


In addition, the polyfit function performs polynomial interpolation:

xs = [1 3 5]';
ys = [2 1 4]';
M = [xs.^2 xs.^1 xs.^0];
c1 = M \ ys;    % The polynomial is c(1)*x^2 + c(2)*x + c(3)
c1 =
0.5000
-2.5000
4.0000
c2 = polyfit( xs, ys, 2 )
c2 =
0.5000   -2.5000    4.0000


At this point, you should realize that Matlab uses vectors for defining polynomials. The format is straight-forward: a vector c of length defines the polynomial

.

Thus, as an example, the vector v = [1 5 3 9] (or v = [1 5 3 9]') defines the polynomial

.

The roots of a polynomial are found with (what else?) the roots function:

v = [1 5 3 9]';
roots( v )
ans =
-4.7667
-0.1166 + 1.3691i
-0.1166 - 1.3691i


### 4.2.4 Evaluating Polynomial in Matlab

The evaluation of a polynomial could proceed as follows:

Calculate each term independently and add them together.

This is exactly what happens when you previously calculated

xr = linspace( 0, 6, 100 );
yr = c(1)*xr.^2 + c(2)*xr + c(3);


Unfortunately, this operation is ; that is, if you have a polynomial of degree , it takes more than linear time to solve it. This is because it becomes more and more expensive to calculate for higher and higher values of . A naive algorithm for calculating is and therefore, an algorithm as poorly written as

function [y] = evaluate_1( p, x )
n = length( p );

y = 0;

for k=(n - 1):-1:0
% Calculate x^k
xk = 1;

for ell=1:k
xk = xk .* x;
end

% Add the term to the solution
y = y + p(n - k) * xk;
end
end


would run in time.

It is possible to calculate x^k in time using the algorithm

function [xk] = my_power( x, k )
if k == 0
xk = 1;
return;
elseif rem( k, 2 ) == 0
xk = my_power( x, k/2 );
xk = xk .* xk;
else
% rem( k, 2 ) == 1
xk = my_power( x, (k - 1)/2 );
xk = xk .* xk .* x;
end
end


We could now write the polynomial evaluation routine as

function [y] = evaluate_2( p, x )
n = length( p );

y = 0;

for k=(n - 1):-1:0
y = y + p(n - k)*my_power( x, k );
end
end


but this would still, therefore, take time.

We can finally reduce this

function [y] = evaluate_3( p, x )
n = length( p );

y = p(1);

for k=2:n;
y = y .* x + p(k);
end
end


In order to demonstrate the effect on the run times, consider evaluating an array of size 1000 at a polynomial with random coefficients of degree 1000, 2000, 4000, ..., 128 0000. The values are shown in Table 1.

DegreeRun Time (s)
evaluate_1evaluate_2evaluate_3
1 000 0.752073 0.2079830.005031
2 000 5.518213 0.3673520.005797
4 000 58.493659 0.7260440.007369
8 000425.452002 1.4801130.014426
16 000n/a 2.8718500.029272
32 000n/a 5.7164110.057208
64 000n/a 11.3753100.114740
128 000n/a 22.6735460.228557

To demonstrate the asymptotic properties, consider Figures 5 through 7 which show these points plotted. While an analysis of the algorithm suggests the run-time of evaluate_1 is quadratic, the image in Figure 5 is actually cubic. Thus, there must be additional complexity in the Matlab functions. The data in Figure 6 is slightly super-linear while the data in Figure 7 appears to be linear.

Students interested in understanding how we can fit a cubic polynomial to the data in Figure 5 are invited to read fitting functions, an optional sub-topic. We will discuss best-fitting least-squares curves in Laboratory 9.

Figure 5. The run time for evaluate_1 for polynomials of degree 1000, 2000, 4000 and 8000.

Figure 6. The run time for evaluate_2 for polynomials of degree 1000, 2000, 4000 through 128 000.

Figure 7. The run time for evaluate_3 for polynomials of degree 1000, 2000, 4000 through 128 000.

It follows that the use of the correct algorithm can have significant impact on the run time.

An easier way of evaluating a polynomial is the Matlab polyval function (a play on the phonetic overlap of polynomial evaluate).

xr = linspace( 0, 6, 100 );
yr = polyval( c, xr );           % evaluate the polynomial defined by the coefficients of c


For example, if we wanted to evaluate the polynomial at each point in the vector xr, we would use

xr = linspace( 0, 6, 100 );
yr = polyval( [4.9e-2 -7.3 8.1 6.5], xr );