MATLAB: How to use Newton-Ramphson method to find an equation root

iterationmatrix codenewton's method

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;

Best Answer

λ 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;
Related Question