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.
