Skip to the content of the web site.
## Laboratory 4

# 4.1 Fixed-point Iteration

## 4.1.1 Implementation in Matlab

#### 4.1.1.1 If Statements

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

#### 4.1.1.3 The `throw` Function and Exceptions

## 4.1.2 Examples of Fixed-Point Iterations

## 4.1.3 Possible Outcomes Fixed-Point Iterations

# 4.2 Finding Interpolating Polynomials

## 4.2.1 The Theory

### 4.2.1.1 Interpolating Two Points

### 4.2.1.2 Interpolating Three Points

### 4.2.1.3 Interpolating Points

## 4.2.2 Interpolating in Matlab

### 4.2.2.1 Element-wise Exponentiation

### 4.2.2.2 Joining Matrices

### 4.2.2.3 Interpolation with Matlab

### 4.2.2.4 Plotting the Points and the Interpolating Quadratic

### 4.2.3 Actual Polynomial Tools in Matlab

### 4.2.4 Evaluating Polynomial in Matlab

This laboratory will investigate two topics:

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

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 .

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:

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:

ifpredicate% 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:

ifpredicate #1% body of code to be executed if the predicate #1 is true elseifpredicate #2% body of code to be executed if predicate #1 is false but predicate #2 is true elseifpredicate #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.

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:

- Determine if the command is too complex, in which case, it can be broken down into a sequence of commands, or
- 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.

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:`which identifies the exception. In this course, the identifier will always be provided.*string* - 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).

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.

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 .

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.

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 .

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 .

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.

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.

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

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

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

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.

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

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.

Degree | Run Time (s) | ||
---|---|---|---|

evaluate_1 | evaluate_2 | evaluate_3 | |

1 000 | 0.752073 | 0.207983 | 0.005031 |

2 000 | 5.518213 | 0.367352 | 0.005797 |

4 000 | 58.493659 | 0.726044 | 0.007369 |

8 000 | 425.452002 | 1.480113 | 0.014426 |

16 000 | n/a | 2.871850 | 0.029272 |

32 000 | n/a | 5.716411 | 0.057208 |

64 000 | n/a | 11.375310 | 0.114740 |

128 000 | n/a | 22.673546 | 0.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 *poly*nomial *eval*uate).

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 );