Skip to the content of the web site.

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, tex:$$f(x)$$, if we start with an initial point tex:$$x_1$$ and then iterate

tex:$$x_k = f(x_{k - 1})$$

for tex:$$k = 2, 3, 4, \ldots$$, then one of two possible results may occur:

  • Either the sequence will converge to a point tex:$$x^*$$, or
  • It will not converge.

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

tex:$$x = f(x)$$.

In order to quickly demonstrate this, consider the cosine function: If we plot tex:$$x$$ and tex:$$\cos(x)$$, we note that they cross somewhere around tex:$$x = 0.75$$. If we start with tex:$$x_1 = 0.75$$ and iterate, we get

tex:$$x_1 = 0.75$$
tex:$$x_2 = \cos(x_1) \approx 0.731688868873821$$
tex:$$x_3 = \cos(x_2) \approx 0.744047084788764$$
tex:$$x_4 = \cos(x_3) \approx 0.735733618187236$$
tex:$$x_5 = \cos(x_4) \approx 0.741338598887922$$
tex:$$x_6 = \cos(x_5) \approx 0.737565296339266$$

Notice that it does seem to be approaching the point tex:$$x^* \approx 0.739085133215161$$ which is the point at which tex:$$x = cos(x)$$.

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 tex:$$|x - x_{\rm prev}| < 10^{-10}$$.

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 tex:$$tan(x)$$ function does not converge with the initial point tex:$$x_1 = 1$$:

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 tex:$$x_1 = 1$$, it does converge to the solution of tex:$$x = cos(x)$$ shown in Figure 1.


Figure 1. The plot of the cosine function and tex:$$y(x) = x$$.

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 tex:$$x^* = 0.73908513321516064166$$.

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 tex:$$x^* = 0$$ is a solution to the equation tex:$$\sin(x) = x$$ 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 tex:$$x = f(x)$$,
  • 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 tex:$$\epsilon_{step}$$ 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 tex:$$N$$.

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 tex:$$(x_1, y_1)$$ and tex:$$(x_2, y_2)$$, it is possible to find a straight line which passes through the two points so long as tex:$$x_1 \ne x_2$$. This is shown in Figure 2 where the polynomial tex:$$p(x) = 0.5x + 1.5$$ passes through the points tex:$$(1,2)$$ and tex:$$(5,4)$$.


Figure 2. A line interpolating two points.

Recall from Linear Algebra that to find the straight line tex:$$p(x) = a_1x + a_0$$ which passes through two points, we must satisfy two equations:

tex:$$a_1 x_1 + a_0 = y_1$$
tex:$$a_1 x_2 + a_0 = y_2$$

which may be written as

tex:$$\left( \matrix{ x_1 & 1 \cr x_2 & 1 } \right )\left( \matrix{ a_1 \cr a_0 } \right ) = \left( \matrix{ y_1 \cr y_2 } \right )$$.

A quick inspection shows that the determinant is tex:$$x_1 - x_2$$ which is zero if and only if tex:$$x_1 = x_2$$.

4.2.1.2 Interpolating Three Points

Given three points tex:$$(x_1, y_1)$$, tex:$$(x_2, y_2)$$ and tex:$$(x_3, y_3)$$, it seems reasonable, also, that a quadratic function could be passed through the three points. For example, the quadratic function tex:$$p(x) = 0.5x^2 - 2.5x + 4$$ passes through the three points tex:$$(1, 2)$$, tex:$$(3, 1)$$ and tex:$$(5, 4)$$, 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

tex:$$a_2 x_1^2 + a_1 x_1 + a_0 = y_1$$
tex:$$a_2 x_2^2 + a_1 x_2 + a_0 = y_2$$
tex:$$a_2 x_3^2 + a_1 x_3 + a_0 = y_3$$

which may be written as

tex:$$\left( \matrix{ x_1^2 & x_1 & 1 \cr x_2^2 & x_2 & 1 \cr x_3^2 & x_3 & 1 } \right )\left( \matrix{ a_2 \cr a_1 \cr a_0 } \right ) = \left( \matrix{ y_1 \cr y_2 \cr y_3 } \right )$$.

A quick inspection shows that the determinant is tex:$$(x_1 - x_2)(x_1 - x_3)(x_2 - x_3)$$ which, again, is zero if and only if tex:$$x_1 = x_2$$, tex:$$x_1 = x_3$$, or tex:$$x_2 = x_3$$.

4.2.1.3 Interpolating tex:$$n$$ Points

Given tex:$$n$$ points tex:$$(x_1, y_1), (x_1, y_1), \ldots, (x_n, y_n)$$, is it possible to find a polynomial of degree tex:$$n - 1$$ 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, tex:$$(3.2, 4.7)$$ and tex:$$(5.0, 1.8)$$, we must solve the linear system

tex:$$\left( \matrix{ 3.2 & 1 \cr 5.0 & 1 } \right )\left( \matrix{ a_1 \cr a_0 } \right ) = \left( \matrix{ 4.7 \cr 1.8 } \right )$$.

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 tex:$$p(x) = -1.6111x + 9.8556$$.

In order to interpolating the three points, tex:$$(3.2, 4.7)$$, tex:$$(5.0, 1.8)$$ and tex:$$(9.4, 6.5)$$, we must solve the linear system

tex:$$\left( \matrix{ 3.2^2 & 3.2 & 1 \cr 5.0^2 & 5.0 & 1 \cr 9.4^2 & 9.4 & 1 } \right )\left( \matrix{ a_2 \cr a_1 \cr a_0 } \right ) = \left( \matrix{ 4.7 \cr 1.8 \cr 6.5 } \right )$$.

[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 tex:$$p(x) = 0.4321x^2 - 5.1547 + 16.7699$$.

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 tex:$$n$$ defines the polynomial

tex:$$c_1 x^{n - 1} + c_2 x^{n - 2} + \mdots + c_{n - 1} x + c_n$$.

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

tex:$$x^3 + 5 x^2 + 3 x + 9$$.

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 tex:$$\omega(n)$$; that is, if you have a polynomial of degree tex:$$n$$, it takes more than linear time to solve it. This is because it becomes more and more expensive to calculate tex:$$x^k$$ for higher and higher values of tex:$$k$$. A naive algorithm for calculating tex:$$x^k$$ is tex:$$\Theta(k)$$ 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 tex:$$\Theta(n^2)$$ time.

It is possible to calculate x^k in tex:$$\Theta(ln(k))$$ 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 tex:$$\Theta(n \ln(n))$$ 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 tex:$$0.049x^3 - 7.3x^2 + 8.1x + 6.5$$ 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 );