How to use Newton-Ramphson method to find an equation root?

Hi, I want to make all of them one M-file by making the second and the third M-file as anonymous functions in the main M-file( f and df are defined as anonymous functions).

THE CURRENT MAIN CODE

        function [p0,err,k,y]=newtonlab(f,df,p0,delta,epsilon,max1)

%Input    - f is the object function 
%            - df is the derivative of f 
%            - p0 is the initial approximation to a zero of f
%           - delta is the tolerance for p0
%           - epsilon is the tolerance for the function values y
%           - max1 is the maximum number of iterations
%Output - p0 is the Newton-Raphson approximation to the zero
%           - err is the error estimate for p0
%           - k is the number of iterations
%           - y is the function value f(p0)

%If f and df are defined as M-file functions use the @ notation
% call [p0,err,k,y]=newton(@f,@df,p0,delta,epsilon,max1).
%If f and df are defined as anonymous functions use the
% call  [p0,err,k,y]=newton(f,df,p0,delta,epsilon,max1).

%  NUMERICAL METHODS: Matlab Programs

for k=1:max1  
  p1=p0-f(p0)/df(p0);  
  err=abs(p1-p0);
  relerr=2*err/(abs(p1)+delta);
  p0=p1;
  y=f(p0);
  if (err<delta)|(relerr<delta)|(abs(y)<epsilon),break,end
end
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

*SECOND M_FILE *_ITALIC TEXT_

      function y=f(x)

%%input to function

    format short e;
            Params = load('saved_data.mat');
            theta = pi/2;
            zeta = cos(theta);
            I = eye(Params.n,Params.n);
            Q = zeta*I-Params.p*Params.p';

%T is a matrix(5,5)

            Mroot = Params.M.^(1/2);  %optimization
            T = Mroot*Q*Mroot;
  %find the eigen values

            E =real( eig(T));

%find the negative eigen values -- %find the smallest negative eigen value

            gamma = min(E);  

%now solve for lambda
          M_inv = inv(Params.M);  %optimization
          zm = Params.zm;

           y= zm'*M_inv*inv(M_inv+x.*Q)*M_inv*zm;

THE THIRD M_FILE:

      function y1=df(x)

%%%%input to the equation
   Params = load('saved_data.mat');
          theta = pi/2;
          zeta = cos(theta);
          I = eye(Params.n,Params.n);
          Q = zeta*I-Params.p*Params.p';

%T is a matrix(5,5)

            Mroot = Params.M.^(1/2);  %optimization
            T = Mroot*Q*Mroot;

%find the eigen values

            E =real( eig(T));
    %find the negative eigen values -- %find the smallest negative eigen value

            gamma = min(E);  
  %now solve for lambda

             M_inv = inv(Params.M);  %optimization
            zm = Params.zm;

            y1=-zm'*M_inv*inv(M_inv+x.*Q)*Q*M_inv*zm;

[Walter Roberson]

Is this question asking something different than your previous question, http://www.mathworks.com/matlabcentral/answers/22711-newton-raphson-method ? If not then the two should be merged together and one of them deleted.

Answer 1

[Walter Roberson]

λ is not a valid variable name in MATLAB. This code would never have passed m-lint. Please do basic error checking on your code before you post it.

It is fairly unlikely that a program would use both fzero() and Newton Raphson.

[EDIT: Dec 4 2011 17:42 CDT - put in explicit merged code]

function newtonlab_driver
  [p0, err, k, y] = newtonlab(@f, @df, 0, 100, 1e-8, 1000);
  fprintf(1, 'Best answer was %g at lambda = %g\n', p0, y);

        function [p0,err,k,y]=newtonlab(f,df,p0,delta,epsilon,max1)

%Input    - f is the object function 
%            - df is the derivative of f 
%            - p0 is the initial approximation to a zero of f
%           - delta is the tolerance for p0
%           - epsilon is the tolerance for the function values y
%           - max1 is the maximum number of iterations
%Output - p0 is the Newton-Raphson approximation to the zero
%           - err is the error estimate for p0
%           - k is the number of iterations
%           - y is the function value f(p0)

%If f and df are defined as M-file functions use the @ notation
% call [p0,err,k,y]=newton(@f,@df,p0,delta,epsilon,max1).
%If f and df are defined as anonymous functions use the
% call  [p0,err,k,y]=newton(f,df,p0,delta,epsilon,max1).

%  NUMERICAL METHODS: Matlab Programs

for k=1:max1  
  p1=p0-f(p0)/df(p0);  
  err=abs(p1-p0);
  relerr=2*err/(abs(p1)+delta);
  p0=p1;
  y=f(p0);
  if (err<delta)|(relerr<delta)|(abs(y)<epsilon),break,end
end

      function y=f(x)

%%input to function

    format short e;
            Params = load('saved_data.mat');
            theta = pi/2;
            zeta = cos(theta);
            I = eye(Params.n,Params.n);
            Q = zeta*I-Params.p*Params.p';

%T is a matrix(5,5)

            Mroot = Params.M.^(1/2);  %optimization
            T = Mroot*Q*Mroot;
  %find the eigen values

            E =real( eig(T));

%find the negative eigen values -- %find the smallest negative eigen value

            gamma = min(E);  

%now solve for lambda
          M_inv = inv(Params.M);  %optimization
          zm = Params.zm;

           y= zm'*M_inv*inv(M_inv+x.*Q)*M_inv*zm;

      function y1=df(x)

%%%%input to the equation
   Params = load('saved_data.mat');
          theta = pi/2;
          zeta = cos(theta);
          I = eye(Params.n,Params.n);
          Q = zeta*I-Params.p*Params.p';

%T is a matrix(5,5)

            Mroot = Params.M.^(1/2);  %optimization
            T = Mroot*Q*Mroot;

%find the eigen values

            E =real( eig(T));
    %find the negative eigen values -- %find the smallest negative eigen value

            gamma = min(E);  
  %now solve for lambda

             M_inv = inv(Params.M);  %optimization
            zm = Params.zm;

            y1=-zm'*M_inv*inv(M_inv+x.*Q)*Q*M_inv*zm;

[zayed]

I have checked the code by m-lint,but it gives error that xi is undefined ,because i don't how to guess initial value

[Walter Roberson]

You might as well guess (-1/gamma)/2

[zayed]

It's always give :

??? Subscript indices must either be real positive integers or logicals.

[Walter Roberson]

fx and dfx need to be functions, not variables, and you should not be overwriting the fx and dfx that are passed in to the routine (probably in your code you should not be passing anything in to your driver routine.)

Your functions have lost lambda which should still be what you are trying to solve for.

Your dfx is the same as your fx, which is wrong. Your dfx needs to calculate the derivative at the current point, and by using the same calculation for both you are implicitly defining an exponential function as only the exponential function can have a derivative that is equal to its function value.

I really do not know how you plan to calculate the derivative for your function. It might be do-able given specific matrices. If you were going to go through that trouble, though, you might as well calculate the derivative of f(x)^2 and solve that derivative for 0, as that calculation is a way of solving for the roots of f(x).

[zayed]

I have edited the code in the thread,but when i run it gives complex outputs for x in a continuous way.Is this code in his way to find the roots of f.Also I have put f & df in the code as you see they are different only in negative sign and Q in df.

[Walter Roberson]

In your previous threads, the version with the Q and the minus sign was what you were trying to find the zero of. Are you saying now that that expression is for the derivative, and thus before you were solving for the derivative being 0, but now you are instead solving for "the original function" being 0 ?

Your calculation of the initial value of f and df is equivalent to calculating with lambda = 0, which you had determined in your previous thread leads to a complex result. As your code has lost its lambda, it is not surprising that your f/df stays complex and thus that you are always dealing with complex.

Your code should have the lambda still in it, and the x of the Newton formula should be rewritten with a change of variable to use lambda.

Really, though, it would have been easier if you had stayed with the earlier code version that calculated all the way down to the Lz function handle, and passed that handle in to a generic Newton-Raphson routine. You would have had to construct a second function to calculate the derivative, but at least then you would not be getting so confused in your code about what x is and what xi is.

Your code really should not be adding 1 to x !

[zayed]

I have edited the code and I did my best to follow your information,But Newtons-Raphsons method is easy to implement in Mathematica but in Matlab it seems a bit difficult. I don't get if I can pass a function to a function and how to use the derivative as a function.

[Walter Roberson]

You cannot pass a function handle in MATLAB and then take the derivative of the represented function. You would have to pass a symbolic function and then take the derivative of the symbolic function, symbolically.

Record your f and df lines as

f = @(lambda) (zm'*M_inv*inv(M_inv+lambda.*Q)*M_inv*zm);

dfdx = @(lambda) (zm'*M_inv*inv(M_inv+lambda.*Q)*M_inv*zm);

but please check double-check the formulae.

My calculations do not agree with the derivative function that you have given.

[zayed]

Thank you for your response.I have edited the code because there was an error in writing the equation,but it gives error :

??? Input argument "x0" is undefined.

Error in ==> newton at 26

x = x0;

[Walter Roberson]

You wrote your function recursively. Your function "newton" tries to call upon "newton" to do the work. That's a recursive call which would never end if it started going at all. But it doesn't start going because you are trying to launch the entire process by calling newton without providing any arguments.

You had a previous version of your code that had the loop with f/dfdx and so on. You need use code similar to what you had in that posting.

[zayed]

I have edited the code with three separated ones,and I have executed them from the command window and it's work,but I need to merge them in one M-file,so how can call the second and the third from the main one ,make f and df defined as anonymous functions),please

[Walter Roberson]

What is your current code? Showing the breaks between the files.

[zayed]

My code is in the thread ,I have edited the code and I have separated the M-files the have to be merged as anonymous,please clarify,thanks in advance.

[Walter Roberson]

Create a new file, such as newtonlab_driver.m

In that file, put a function newtonlab_driver with no parameters.

Have that function call newtonlab with all appropriate parameters, and have it do whatever displaying of results that you want.

In that same file after the end of the driver, put the code from the other three files. The order will not matter for execution purposes.

Note: your f and df functions do not use gamma and so do not need to compute gamma -- which means you also do not need to compute the eigenvalues in those routines, and that means in turn you do not need to compute T or Mroot.

It is, of course, very inefficient to load() from the file and do the matrix inversions each time, but that's just an efficiency matter and as long as you get the routine working I do not care whether it takes 3 seconds or 3 days to compute the result.

[zayed]

Actually I tried to do it as you said, but when I execute it it does nothing ,so I need to you to see it ,can put it in the thread that I have edited.

[Walter Roberson]

Merged code posted in my Answer. I used some WAGs for the parameters such as tolerance.

I did not make any of the efficiency improvements: time enough for those after you get everything working.

[zayed]

It gives the following error :

"""?? Error: File: C:\MATLAB7\work\project\NEWTON\newtonlab_driver.m Line: 65 Column: 54

The function "newtonlab_driver" was closed

with an 'end', but at least one other function definition was not.

To avoid confusion when using nested functions,

it is illegal to use both conventions in the same file.""""

[Walter Roberson]

Fixed.

[zayed]

What this word mean,please?

[Walter Roberson]

I removed the 'end' statement that was causing the problem. The code is repaired.

http://dictionary.reference.com/browse/fixed

fixed [...], noun verb (used with object)

1. to repair; mend.

[zayed]

The following error :

""??? Undefined command/function 'fprint'.""

[bym]

the command is fprintf

[zayed]

It's working now ,thank you very much ,so how can I get rid of rhe duplication in the code to make it faster,thanks.

[Walter Roberson]

As I wrote above:

Note: your f and df functions do not use gamma and so do not need to compute gamma -- which means you also do not need to compute the eigenvalues in those routines, and that means in turn you do not need to compute T or Mroot.

It is, of course, very inefficient to load() from the file and do the matrix inversions each time, but that's just an efficiency matter and as long as you get the routine working I do not care whether it takes 3 seconds or 3 days to compute the result.

====

For more hints on how to speed up your code look back at some of my previous responses such as http://www.mathworks.com/matlabcentral/answers/22585-code

[zayed]

But how can can I make the initial point (p0) to changes with gamma,so if I want p0 to be (-/gamma)/2,how can I do this in main function,thanks

[Walter Roberson]

Then in your main function you would load the file and calculate gamma, just like you do in your current code.

[zayed]

So I should put 'load file and do the next step involving calculating gamma and also (I,Q,M_inv) to get rid of duplication, on the line after the first line (function newtonlab_driver) in the edited code

[zayed]

How can I call the result of the FUNCTION (newtonlab_driver) -in the thread below-into another code.Also I need to change the intial value guess p0 to changed as gamma changed .

http://www.mathworks.com/matlabcentral/answers/22787-how-to-use-newton-ramphson-method-to-find-an-equation-root

Thanks