Rooting Around in MATLAB – Part 3

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

Contents

Set Up

clearallcloseall

例子提醒

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

which we define using an anonymous function.

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

First Fixed Point Iteration Attempt

Previously we rewrotef(x)=0so thatxwas 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 withrootsandfzero.

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 rearrangef(x)=0into something likeg2(x)=x. We have already triedg1 = 1 - x^3. We can also rewrite the equation as

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

Let's Iterate

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

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

First, we'll select a starting valuex = 0.5, andy, 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 ofg2to the newxvalue, sliding onto the straight line, and then calculate the next value ofg2from this new value forx. 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;forn = 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));endplot(x,y,'r') holdoff

Fixed Point Found!

We appear to be converging on a fixed point, since after iterating, the final values forxandg2(x)(which isy) 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 functionsg1andg2made 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 thoughtshere.

Get the MATLAB code

Published with MATLAB® 7.8