Rooting Around in MATLAB - Part 1 | Loren on the Art of MATLAB

June 11th, 2009

Rooting Around in MATLAB - Part 1

I've been interested in teaching for a long time, including ways to use MATLAB. One concept that students might need to understand early in their college careers is that of finding roots or zeros of functions. To understand at least some of the algorithms, you might want to teach the students about fixed points for functions. It's the basis for some methods of solving equations or finding roots, algorithms such as Newton's method, finding square roots, and more.

Example Function
Find the Roots or Zeros
The Series of Posts

Example Function

Let's start with a simple cubic polynomial f.

Here's one common way to represent this polynomial in MATLAB, using coefficients of descending powers of the independent variable.

p = [1 0 1 -1];

I can then use polyval to evaluate the polynomial. And then I can plot it.

x = -2:0.1:2;
y = polyval(p,x);
plot(x,y)
title f
grid

I could also represent the polynomial as an anonymous function and plot it with fplot.

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

Find the Roots or Zeros

I have at least 2 choices in MATLAB for finding a zero or root of this polynomial. The first is to use roots to get all the possible zeros.

rsolution = roots([1 0 1 -1])

rsolution =
     -0.34116 +     1.1615i
     -0.34116 -     1.1615i
      0.68233

You can see that this polynomial has one real root between 0 and 1, and 2 complex roots.

You can also use fzero, one of MATLAB's optimization functions, to find the value. Here we'll select 0.5 as the initial guess.

fzsolution = fzero(f,0.5)

fzsolution =
      0.68233

In the next post, I'll describe a way to solve this same problem using an algorithm based on fixed point iteration.

The Series of Posts

In addition to this post, there will be two more. Until they are published, the following two links will not be available.

Get the MATLAB code

Published with MATLAB® 7.8

By Loren Shure

You can follow any responses to this entry through the RSS 2.0 feed. You can skip to the end and leave a response. Pinging is currently not allowed.

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

Yi Cao replied on June 12th, 2009 at 14:46 UTC :

Another way to get the root:

x0=0.5;
x1=0;
while abs(x1-x0)>1e-6,
   x1=x0;
   x0=(1-x1)^(1/3);
end
disp(x0)

Loren replied on June 12th, 2009 at 14:55 UTC :

Yao-

True, many ways, but I am thinking of how to derive and explain an algorithm to students. How would someone know know to use (1-x)^(1/3)? Wait for my next post to my ideas on how to explore and teach some of these concepts.

–Loren

Yi Cao replied on June 12th, 2009 at 17:06 UTC :

Loren,

I do not know what you are going to write for Part 2 and Part 3. For me, the logic to use (1-x)^(1/3) is clear since the original equation x^3+x-1=0 is equivalent to x=(1-x)^(1/3). Therefore, the iteration is naturally derived as:

x_k+1 = (1-x_k)^(1/3)

Certainly, there are many ways to work out an iteration like this, e.g. from x = 1 - x^3, we can have iteration

x_k+1 = 1 - x_k^3

However, this iteration is unstable for 3*x_k^2 > 1, i.e. x_k > 0.5774. Since the solution is x=0.68233 > 0.5774, the second iteration is not suitable for this problem.

The advantage of my solution is that

1. It works for fractional order, e.g. x^3.25 + x - 1 =0, which cannot be solved by roots function.
2. It can be vectorised for different constants other than 1:

x^3.25 + x - a = 0, for a = 1:10.

a=1:10;
x0=0.5;
x1=0;
while norm(x1-x0)>1e-6,
   x1=x0;
   x0=(a-x1).^(1/3.25);
end
disp(x0)
err=x0.^3.25+x0-a

If we use fzero:

fz=zeros(1,10);
for a=1:10
  f=@(x)x.^3.25+x-a;
  fz(a)=fzero(f,0.5);
end
disp(fz)
err=fz.^3.25+fz-(1:10)

Only the first a=1 is solved. Other a’s end up with NaN!

Loren on the Art of MATLAB

June 11th, 2009