Loren on the Art of MATLAB

June 23rd, 2009

Rooting Around in MATLAB - Part 3

I recently wrote a pair of posts (1 and 2) about finding roots of equations, and how you might use MATLAB to explore this topic with emphasis on the method of fixed point iteration.

Contents

Set Up

clear all
close all

Example Reminder

Let me restate the example function. Let's start with a simple cubic polynomial f(x)

which we define using an anonymous function.

f = @(x) x.^3 + x - 1;
fplot(f,[-2 2])
title f
grid on

First Fixed Point Iteration Attempt

Previously we rewrote f(x)=0 so that x was alone on one side of the equation.

g1 = @(x) 1 - x.^3;

This function, g1(x) did not help find the root between 0 and 1 - every step took us further away from the solutions we found with roots and fzero.

fzsolution = fzero(f,0.5)
fzsolution =
      0.68233

Second Way to Rewrite Equation

There's another way we can write the equation for the fixed point. Remember, we want to rearrange f(x)=0 into something like g2(x)=x. We have already tried g1 = 1 - x^3. We can also rewrite the equation as

g2 = @(x) 1./(x.^2+1);
fplot(g2,[0 1]);
hold on
straightLine = @(x) x;
fplot(straightLine, [0 1], 'g')
legend('g2','x','Location','SouthEast')
grid on
axis equal, axis([0 1 0 1])

Let's Iterate

Following the same idea from the last post, we now iterate, finding successive guesses of x, computing the value of g2(x), setting this value to be the next estimate of x, etc.

Let's Try to "Zero" in on the Solution

First, we'll select a starting value x = 0.5, and y, from the straight line also at 0.5.

x(1) = 0.5;
y(1) = x(1);
x(2) = x(1);
y(2) = g2(x(2));

And let's plot the trajectory of the solutions, starting with the first point.

plot(x,y,'r',x(1),y(1),'r*')

First Iteration

As described in the previous post, I now set the current value of g2 to the new x value, sliding onto the straight line, and then calculate the next value of g2 from this new value for x. And plot it.

x(3) = y(2);
y(3) = y(2);
x(4) = x(3);
y(4) = g2(x(4));
plot(x,y,'r')

Second Iteration

Here's the second iteration.

n = 5;
x(n) = y(n-1);
y(n) = y(n-1);
x(n+1) = x(n);
y(n+1) = g2(x(n+1));
plot(x,y,'r')

Third Iteration

And the third.

n = 7;
x(n) = y(n-1);
y(n) = y(n-1);
x(n+1) = x(n);
y(n+1) = g2(x(n+1));
plot(x,y,'r')

A Few More Rounds

Let's do a few more iterations and preallocate the arrays instead of growing them. And yes, I know in this post I kept overplotting lines successively. I didn't want us to get distracted by managing the plots.

x(9:22) = NaN;
y(9:22) = NaN;
for n = 9:2:21
    x(n) = y(n-1);
    y(n) = y(n-1);
    x(n+1) = x(n);
    y(n+1) = g2(x(n+1));
end
plot(x,y,'r')
hold off

Fixed Point Found!

We appear to be converging on a fixed point, since after iterating, the final values for x and g2(x) (which is y) are getting quite close to each other.

final = [ x(end) y(end) fzsolution]
final =
      0.68039      0.68356      0.68233

Wrapping Up the Series of Posts

This post completes this series of posts on fixed point iteration. There is, of course, more that you could do in a class. For example, you might explore what characteristics of the functions g1 and g2 made the fixed point iteration strategy behave differently. Perhaps look at first derivatives?

Here's the recap of the series

Is this sort of tutorial relevant to your work, especially if you are teaching? How about the incorporation of the visualization during the exploration? Let me know your thoughts here.


Get the MATLAB code

Published with MATLAB® 7.8

4 Responses to “Rooting Around in MATLAB - Part 3”

  1. matt fig replied on :

    I often find with these types of things that a simple animation is worth a thousand words. For instance:

     
    g2 = @(x) 1./(x.^2+1);
    fplot(g2,[.4 .9]);
    hold on
    straightLine = @(x) x;
    fplot(straightLine, [.4 .9], 'g')
    legend('g2','x','Location','SouthEast')
    grid on
    axis equal, axis([.4 .9 .4 .9])
    
    xo = .5; % Initial Guess.
    yo = .5;
    plot(xo,yo,'*r','markersize',6)
    T = title(['Root: ',num2str(yo,'%.6f')],'fonts',14,'fontw','b');
    pause(.4)
    
    for n = 1:22
        x = yo;
        y = g2(xo);
        plot(x,y,'*r','markersize',5)
        set(T,'string',['Root: ',num2str(y,'%.6f')])
        xo = x;
        yo = y;
        pause(.4)
    end
    
    hold off
     
  2. Loren replied on :

    Matt-

    Animations can definitely help. I was trying to illustrate the steps for teaching purposes.

    –loren

  3. Ben replied on :

    Loren, as you said there are many ways to write the fixed point equation. Another interesting way is to run the iterations backwards if possible. Here we can solve g1inv=@(x) (1-x).^(1/3) which converges for x0 in (0,1).

    This technique is also useful for finding unstable fixed points and limit cycles of dynamical systems.

  4. Loren replied on :

    Ben-

    Yes, you have listed the 3rd way to factor the equation. Yi Cao mentioned it as well in a comment on the first part of this series.

    –Loren

Leave a Reply

Wrap code fragments inside <pre> tags, like this:

<pre class="code">
a = magic(3);
sum(a)
</pre>

If you have a "<" character in your code, either follow it with a space or replace it with "&lt;" (including the semicolon).


Loren Shure works on design of the MATLAB language at The MathWorks. She writes here about once a week on MATLAB programming and related topics.

  • Loren: Dirga- You might consider looking at the function interp2. –loren
  • Dirga: I do not know whether the term I used is correct or not. I have set of measurement data. It is the density of...
  • Loren: Alex- You could overload sort in the same way that you create any class method. That would get you the shorter...
  • Alan Hester: can you put infomation in the interlaced images that the subconscious will pick up on and remember that...
  • Andrew Kraev: Hello, here is simple code to calculate for an optional vector a[i], i=1,..,n the sum by i of products,...
  • Alex: It would be nice to be able to sort an array of objects of a custom class, provided the class implements lt()...
  • Stu: Excellent post. No, I wasn’t aware of that property of for...end structures. Thank you!!
  • Loren: Eleanor- You either need an algorithm to determine which 0s should be replaced by NaN values or you will need...
  • Eleanor: Loren- I have a different sort of problem with NaNs and 0s. My data has zeros where data is missing, but...
  • Aslak Grinsted: Here’s another fun piece of code inspired by cellular automata. This one generates what looks...

These postings are the author's and don't necessarily represent the opinions of The MathWorks.