Newton's Method Question

"Use your code to locate the root of x^3 + x = 2*x^2 + 3, starting with x0 = 3 and print a table of at least 5 iterates for epsilon = 1e − 7."

I'm having trouble getting my code right. I know that there's better methods out there, but this is along the lines of how it's supposed to be written for an assignment. I'm new to using Matlab, so it may not even make sense + it's my first time encountering the "while" function. Looking for some feedback on how to make this code more understandable.

%Newton's Method

function root = newtonr(f, fp, x, epsilon, nmax)

n=1;

x=3;

f=x^3-2*x^2+x-3;

fp=3*x^2-4*x+1;

nmax=10;

epsilon=1.0e-15;

while (abs(error)>epsilon && n <= nmax)

      denom=fp;
      if (denom==0)
          disp('denom is zero');
          return;
      else
          x=x-f/denom;
      end
      n=n+1;
      error=x-x;
      x=x;
  end
  root=x;

Answer 1

[James Tursa]

Some hints (I leave it to you to insert them properly into your code):

I assume the function handle and function derivative handle are passed into your newtonr routine, so the f and fp lines in your code should evaluate these functions, not overwrite them. E.g., I would assume something like this would be in the calling routine:

f = @(x)x^3-2*x^2+x-3;
fp = @(x)3*x^2-4*x+1;

And then f and fp would be input arguments to the newtonr routine. So inside the newtonr routine these lines as written overwrite the f and fp inputs:

f=x^3-2*x^2+x-3;
fp=3*x^2-4*x+1;

To make your newtonr routine find roots of a function that is passed in via f and fp, don't overwrite them. E.g., do something like this instead to evaluate the function and its derivative at the value of x:

fx = f(x);
fpx = fp(x);  % <-- but not really needed at the front of newtonr

Also, don't use error for a variable name since it is a built-in function name:

while (abs(error)>epsilon && n <= nmax)

The above line should be testing for change in x instead.

Your denominator evaluation line:

      denom=fp;

Should look something like this instead (re-evaluate the function derivative at the current x):

      denom = fp(x);

And this line:

      x=x-f/denom;

Don't overwrite x immediately. First save the result in a new variable, e.g. xnew, so that you can calculate the change in x later. E.g.,

      xnew = x - fx/denom;

Then the change in x (needed for your convergence test) can be:

      xchange = xnew - x;

Then you can update x:

      x = xnew;
      fx = f(x);

[Zachary Grant]

function root = newtonr(f, fp, x, epsilon, nmax)
n=1;
x=3;
nmax=10;
epsilon=1.0e-15;
while (abs(xchange)>epsilon && n <= nmax)
        denom=fp(x);
        if (denom==0)
            disp('denom is zero');
            return;
        else
            xnew=x-f/denom;
        end
        n=n+1;
        xchange=xnew-x;
        x=xnew;
        f=f(x);
        fp=fp(x);
  end
  end
  root=xnew;

Would it look more like the above? I saved "f" and "fp" as separate .m files in the same directory. But then it says that "f" is an undefined function/variable

[Steven Lord]

Closer, but there are still some issues.

1) Your assignments to x, nmax, and epsilon on lines 3-5 of newtonr overwrite the values of x, nmax, and epsilon passed into your function by the user. This is a little bit of an advanced maneuver but you can assign default values to these variables if the user didn't specify them by checking the number of inputs as returned by the NARGIN function. For instance, the code below will only run the body of the IF statement if the user called the function with 4 or fewer input arguments.

 if nargin < 5 % user did not specify nmax
   nmax = 10;
 end

If you don't want to try that, just remove those three lines from newtonr.

2) In your ELSE block you're trying to divide the function handle f by the denominator. That won't work; you want to divide the value you receive when you evaluate the function handle f by the denominator instead. As a simple example:

 f = @(x) x.^2;
 g = f./2 % will not work
 h = f(4)/2 % will work and will return 8 (= 4^2/2)

3) The lines "f = f(x)" and "fp = fp(x)" overwrite the function handles with the values returned by evaluating them. The next time you try to evaluate them they're no longer function handles, and so those lines will likely error. You never need to change f and fp; you need to evaluate them.

4) You need to call this function with at the very least 2 or 3 inputs. [If you want to use 3 as a default starting point for the method you can require only two inputs; if you want to require the user to specify the starting point you need to require at least three.] Show us how you're calling your function; my guess is that you're calling it with 0 inputs (either by typing just its name at the Command Window prompt or by pressing the green triangle button in the Editor.) That's not going to work; newtonr won't know the function whose root you're trying to compute.

[Zachary Grant]

function x = newtonr(f, fp, x, e, nmax)
n=1;
if nargin < 5
    nmax=10;
end
while (abs(xchange)>e && n <= nmax)
      denom=fp(x);
      if (denom==0)
          disp('denom is zero');
          return;
      else
          xnew=x-f(3)/denom;
      end
      n=n+1;
      xchange=xnew-x;
      x=xnew;
end
end

"x=newtonr(f,fp,x,e,nmax)" is how I'm calling the function. When you say "evaluate," what exactly do you mean? Sorry for my lack of knowledge, I've been trying to teach myself this stuff by just searching around the internet. Once I've changed around some stuff, it now says "Undefined function or variable "xchange"." How do I define it?

[Steven Lord]

How do you define f and fp? Do you define them as anonymous functions, function handles to named functions, or vectors of data? For simple functions, an anonymous function (like James created in the first code block in his Answer) is probably good enough to define the function to be evaluated. But if you have a more complicated function, writing it as a function file and passing a function handle to it into your newtonr function to be evaluated when newtonr needs a value from your function is fine as well.

 fh = @why; % Define a function handle to the WHY function
 fh(1) % This is the same as why(1), but without hard-coding WHY into the code

 fh = @sin; % Change the function to be evaluated
 fh(1) % Now returns sin(1)

As for the reason for the "Undefined function or variable xchange" error -- where is the first place you define that variable in your code? Where is the first place you use that variable? Which one comes first, definition or usage?

[Zachary Grant]

I guess I could just leave it as an anonymous function, like James suggested. I don't define "xchange" until the while function, and I think the usage comes before the definition.