Main Page EEng School DCU


Section 3: Iterative Techniques and Applications (Matlab Examples).

    ·    Iterative Techniques
    ·    Input-Output Approach
    ·    Newton’s Method (Newton-Raphson method)
    ·    Secant Method
    ·    Iterative Calculation of Inverses

In this section iterative techniques are used to find solution of equations. Because formulas that give exact numerical values of the solution to equations will exist only in very simple situations, in most cases we have to use an approximation method. Iteration methods start from an initial guess  (which may be poor) and compute step-by-step (in general, better and better) approximations of an unknown solution of equation.
 

approot.m

We are going to start with a simple iterative algorithm proposed in 1999 by John H. Mathews and Kurtis D. Fink (NUMERICAL METHODS Using MATLAB). The matlab code has been improved here in order to be student-friendly. The code is given below. The algorithm evaluates a vector of approximate locations for roots of a test function  f_test which is saved in a different .m file. The student is free to change this function and approximate its roots!

% Input - f is object function saved as an M-file named f.m
%       - X is the vector of abscissas
%       - epsilon is the tolerance
% Output - R is the vector of approximate locations for roots

while 1

 clear, clf,clc
   fprintf( '\n============================================\n' );
   fprintf( ' Approximative Root Method  \n' );
   i = 1;
   fprintf( '** Give the X abscissa values:\n');

   while 1
      km = input( '' );
      X(i) = km;
      i=i+1;
      kont = input( 'Type 1 to continue, or 0 to stop adding X values: ' );
  if kont == 0, break; end
   end

   while 1
      fprintf( '\n Input the tolerance value, epsilon = ')
      epsilon = input('');
      if epsilon < 1 & epsilon > 0;  break; end
      fprintf( ' \n Epsilon should be greater than 0 and less than 1!!! ')
      fprintf( ' \n Repeat your input for tolerance value! ')
   end

 Y=f_test(X);
 axes;
 plot(X,Y,'r-o');
 centaxes

 yrange = max(Y)-min(Y);
 epsilon2 = yrange*epsilon;
 n=length(X);
 m=0;
 X(n+1)=X(n);
 Y(n+1)=Y(n);

 for k=2:n,
  if Y(k-1)*Y(k) <= 0,
   m=m+1;
   R(m)=(X(k-1)+X(k))/2;
  end
  s=(Y(k)-Y(k-1))*(Y(k+1)-Y(k));
  if (abs(Y(k)) < epsilon2) & (s <= 0),
   m=m+1;
   R(m)=X(k);
  end
 end

 fprintf(' The solutions are: \n');
 if length(R) ~= 0
    for k=1:length(R),
       fprintf(' R(%u) = %6.2f \n', k, R(k));
    end
 else
    fprintf(' No root values! ')
 end

 fprintf('\n============================================')
 kont = input( 'Type 1 to continue, or 0 to stop: ' );
 if kont == 0, break; end
end

We take the following  f_test example: f(x) = x3 + x2 - 1.25x - 0.75 which was used throughout in the previous chapter's examples and in the course notes. The figure drawn by the algorithm is given below:

 The approximative abscissas values where the f(x) function is approaching the X axe are the folowing (R is the vector where they are stored):

R(1) =  -1.75
R(2) =  -1.25
R(3) =  -1.00
R(4) =  -0.75
R(5) =  -0.25
R(6) =   0.75
R(7) =   1.25

It does not mean that function f(x) has 7 roots. The algorithm finds out in which coordinates the f(x) function given in the following X abscissa points:
X0 = -3, X1 = -2.5, X2 = -2, X3 = -1.5, X4 = -1, X5 = -0.5, X6 = 0, X7 = 0.5, X8 = 1, X9 = 1.5, X10 = 2, X11 = 2.5, X12 = 3, X13 = 3.5
with a tolerance epsilon = 0.01
 

secant_alg.m

Next the Secant iterative algorithm is experimented. The algorithm given for it is improved from the function coded in the same "NUMERICAL METHODS Using MATLAB" by John H. Mathews and Kurtis D. Fink. The improved matlab code is given below:

%function [p1,err,k,y]=secant(f,p0,p1,delta,epsilon,max1)
%Input - f is the object function input as a string 'f'
%      - p0 and p1 are the initial approximations to a zero of f
%     - delta is the tolerance for p1
%     - epsilon is the tolerance for the function values y
%     - max1 is the maximum number of iterations
%Output - p1 is the secant method approximation to the zero
%      - err is the error estimate for p1
%      - k is the number of iterations
%      - y is the function value f(p1)

while 1

   clear, clf,clc
   fprintf( '\n============================================\n' );
   fprintf( ' Secant Method \n' );

   while 1
      fprintf( ' \n The absciss extremes a and b (a < b)')
      fprintf( ' \n a = ');  a = input('');
      fprintf( ' b = ');  b = input('');
      if a < b;  break; end
      fprintf( ' \n a should be less than b! ')
      fprintf( ' \n Repeat your input for a and b! ')
   end

   fprintf( ' \n Now the function is plotted \n')
   x=a:0.1:b;
   f = feval('f306',x);
   plot(x, f, 'b-');
   centaxes
   hold on;

   fprintf( '** Give the Po, first initial abscissa value approximation to a zero of f:\n');
   p0 = input('Initial approximation Po = ');

   fprintf( '** Give the P1, second initial abscissa value approximation to a zero of f:\n');
   p1 = input('Initial approximation P1 = ');

   while 1
      fprintf( '\n Input the tolerance value for P1, delta = ')
      delta = input('');
      if delta < 1 & delta > 0;  break; end
      fprintf( ' \n Delta should be greater than 0 and less than 1!!! ')
      fprintf( ' \n Repeat your input for tolerance value! ')
   end

   while 1
      fprintf( '\n Input the tolerance value for the function Y values, epsilon = ')
      epsilon = input('');
      if epsilon < 1 & epsilon > 0;  break; end
      fprintf( ' \n Epsilon should be greater than 0 and less than 1!!! ')
      fprintf( ' \n Repeat your input for tolerance value! ')
   end

   while 1
      fprintf( '\n Input the maximum number of iterations, N_max = ')
      max1 = input('');
      if max1 > 1;  break; end
      fprintf( ' \n The number of iterations should be greater than 1!!! ')
      fprintf( ' \n Repeat your input for the number of iterations value! ')
   end

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

   fprintf('\n==============RESULTS================')

   fprintf('\n P1 is the Secant approximation to zero = %12.6f', p1);
   fprintf('\n Err is the error estimate for P1 = %12.6f', err);
   fprintf('\n K is the number of iterations = %u', k);
   fprintf('\n Y is the function value f(P1) = %12.6f', y);

   fprintf('\n============================================')
   kont1 = input( 'Type 1 to continue, or 0 to stop: ' );
   if kont1 == 0, break; end
end

Explanations on some of the input variables given to the Secant algorithm are written in the comments at the begining of the matlab code. Next, the algorithm is run for the following function example: f(x) = x-cos(x) given in f306.m. The function is given for the a = -10, b = 10 interval and is depicted below:


The other input parameters are: P0 = -1.67, P1 = -2.21, delta = 0.002, epsilon = 0.012, N_max = 100 and the results are taken from the main (command) Matlab window:

P1 is the Secant approximation to zero  = 0.737000
Err is the error estimate for P1 = 0.064575
K is the number of iterations = 5
Y is the function value f(P1) = -0.003488

The initial Secant function from whichthe above code was developed is given in secant.m.
 

newt_n.m

The following iterative method presented is Newton. Its matlab code is devoled on the code written by Nakamura: newt_n_func.m. The whole code is given below. The algorithm is run for the same function:  f(x) = x-cos(x) given in f306.m. Another Newton function code is written in fnewton.m.

% Newt_n(f_name, x0) finds a root of a function by
% Newton iteration
% f_name: the function name that defines the equation
%         to solve
% x0: an initial guess.
% improved after Copyright S. Nakamura, 1995
% function x = Newt_n(f_name, x0)
while 1
   clear, clf,clc

   fprintf( '\n============================================\n' );
   fprintf( ' Newton Iteration Method (without graphics) \n' );
   fprintf( '** Give the Xo abscissa value approximation:\n');
   x0 = input('Initial guess Xo = ');

   x = x0; xb=x-999;
   n=0;  del_x = 0.01;

 while abs(x-xb)>0.000001
     n=n+1;   xb=x;
     if n>300 break; end
     y=feval('f306', x);
     y_driv=(feval('f306', x+del_x) - y)/del_x;
     x = xb - y/y_driv;
     fprintf(' n=%3.0f, x=%12.5e, y = %12.5e, ', n,x,y)
     fprintf(' yd = %12.5e \n', y_driv)
 end

   fprintf('\n     Final answer = %12.6e\n', x);
   fprintf('\n============================================')
   kont = input( 'Type 1 to continue, or 0 to stop: ' );
   if kont == 0, break; end
end

The initial X0 abscissa value approximation and the results are given below exactly the same way they are listed by the algorithm in the Matlab command window:

============================================
 Newton Iteration Method (without graphics)
** Give the Xo abscissa value approximation:
Initial guess Xo = -1.23
 n=  1, x=2.51938e+001, y = -1.56424e+000,  yd = 5.91981e-002
 n=  2, x=2.49622e+000, y = 2.41957e+001,  yd = 1.06600e+000
 n=  3, x=4.33557e-001, y = 3.29509e+000,  yd = 1.59750e+000
 n=  4, x=7.66219e-001, y = -4.73920e-001,  yd = 1.42463e+000
 n=  5, x=7.39300e-001, y = 4.56807e-002,  yd = 1.69701e+000
 n=  6, x=7.39086e-001, y = 3.59955e-004,  yd = 1.67745e+000
 n=  7, x=7.39085e-001, y = 8.07423e-007,  yd = 1.67730e+000

     Final answer = 7.390851e-001

============================================Type 1 to continue, or 0 to stop:
 

newt_g.m

The above Newton method is re-analyzed in this code example but this time with graphical features as well. The code is given below:

% Newt_g(f_name, x0, xmin, xmax, n_points)
% finds a root of a function f_name by
% newton iteration (with graphics)
% improved after Copyright S. Nakamura, 1995
while 1

   clear, clf,clc
   fprintf( '\n============================================\n' );
   fprintf( ' Newton Iteration Method (with graphics) \n' );
   fprintf( '** Give the Xmin and Xmax abscissa values:\n');

   while 1
      xmin = input( 'Xmin = ' );
      xmax = input( 'Xmax = ' );
      if xmin < xmax;  break; end
      fprintf( ' \n Xmin has to be less than Xmax!!! ')
      fprintf( ' \n Repeat your input for Xmin and Xmax values! ')
   end

   fprintf( '\n Input the number of points\n ')
   n_points = input('N_points = ');

   fprintf( '\n Input an initial guess\n ')
   x0 = input('Xo = ');

   clf, hold off
   % xmin, xmax:  min and max for plotting
   del_x=0.001;
   wid_x = xmax - xmin;  dx = (xmax- xmin)/n_points;
   xp=xmin:dx:xmax;   yp=feval('f306', xp);
   plot(xp,yp); xlabel('x');ylabel('f(x)');
   title('Newton Iteration'),hold on
   ymin=min(yp); ymax=max(yp);wid_y = ymax-ymin;
   yp=0.*xp;  plot(xp,yp)
   x = x0;   xb=x+999; n=0;
   while abs(x-xb)>0.000001
      if n>300 break; end
      y=feval('f306', x);
      plot([x,x],[y,0],'--');
      plot(x,0,'o')
      fprintf(' n=%3.0f,  x=%12.5e,  y=%12.5e\n', n,x,y);
      xsc=(x-xmin)/wid_x;
      if n<4,
        text(x, -wid_y/20, [ num2str(n)], 'FontSize', [16]),
     end
     y_driv=(feval('f306', x+del_x) - y)/del_x;
     xb=x;
     x = xb - y/y_driv; n=n+1;
     plot([xb,x],[y,0],':')
  end
  plot([x x],[0.02*wid_y 0.2*wid_y])
  text( x, 0.25*wid_y, 'Final solution','FontSize', [16])
  plot([x (x-wid_x*0.004)],[0.02*wid_y  0.05*wid_y])
  plot([x (x+wid_x*0.004)],[0.02*wid_y  0.05*wid_y])
  axis([xmin,xmax,ymin*1.2,ymax*1.2])

  fprintf('\n============================================')
  kont = input( 'Type 1 to continue, or 0 to stop: ' );
  if kont == 0, break; end
end

As an example the above algorithm is run having the input and results depicted below exactly like in the Matlab command window:

============================================
 Newton Iteration Method (with graphics)
** Give the Xmin and Xmax abscissa values:
Xmin = -10
Xmax = 10

 Input the number of points
 N_points = 30

 Input an initial guess
 Xo = -1.56
 n=  0,  x=-1.56000e+000,  y=-1.57080e+000
 n=  1,  x=2.46019e+004,  y=2.46029e+004
 n=  2,  x=-2.79821e+003,  y=-2.79763e+003
 n=  3,  x=1.23465e+004,  y=1.23455e+004
 n=  4,  x=7.16039e+002,  y=7.15068e+002
 n=  5,  x=-2.26678e+002,  y=-2.27563e+002
 n=  6,  x=1.98154e+002,  y=1.99127e+002
 n=  7,  x=-6.11256e+001,  y=-6.09906e+001
 n=  8,  x=-3.04890e+001,  y=-3.10893e+001
 n=  9,  x=-1.32181e+001,  y=-1.40131e+001
 n= 10,  x=2.23598e+001,  y=2.32926e+001
 n= 11,  x=-1.40845e+001,  y=-1.41371e+001
 n= 12,  x=9.98146e+003,  y=9.98227e+003
 n= 13,  x=-1.39418e+004,  y=-1.39427e+004
 n= 14,  x=-4.87296e+003,  y=-4.87202e+003
 n= 15,  x=-1.23669e+003,  y=-1.23715e+003
 n= 16,  x=-5.81733e+002,  y=-5.80875e+002
 n= 17,  x=-1.97658e+002,  y=-1.96692e+002
 n= 18,  x=6.80764e+001,  y=6.75690e+001
 n= 19,  x=-4.19595e+002,  y=-4.19785e+002
 n= 20,  x=-2.07767e+002,  y=-2.08679e+002
 n= 21,  x=1.45265e+002,  y=1.44534e+002
 n= 22,  x=5.93973e+001,  y=6.03547e+001
 n= 23,  x=1.25486e+001,  y=1.15488e+001
 n= 24,  x=7.97007e-001,  y=9.81568e-002
 n= 25,  x=7.39794e-001,  y=1.18581e-003
 n= 26,  x=7.39085e-001,  y=4.46683e-007

============================================Type 1 to continue, or 0 to stop:

The graphical representation of the algorithm run is given next:


newton_raphson_polynom.m

The next method proposed here is the one proposed by Newton-Raphson. The below Matlab code is an extension of the function proposed in newton_raphson.m (proposed in "NUMERICAL METHODS Using MATLAB" by John H. Mathews and Kurtis D. Fink) and is dedicated to the particular case of polynomial functions because their analytical first derivatives can be calculated easily.

%function [p0,err,k,y]=newton(f,df,p0,delta,epsilon,max1)
%Input - f is the object function input as a string 'f'
%      - df is the derivative of f input as a string 'df'
%      - 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)

while 1
   clear, clf,clc
   fprintf( '\n============================================\n' );
   fprintf( ' Newton-Raphson Method. Polynomial examples \n' );

   while 1
      km = input( '** The degree of the polynom ? N = ' );
      if km>1 break; end
      fprintf(' Input is invalid: Repeat.\n')
   end

   while 1
      fprintf( '\n Input the polynomial coefficients ')
      fprintf( '\n like [x^N x^N-1 ... x^1 x^0]');  coef = input('');
      if length(coef) == km+1;  break; end
      fprintf( ' \n Number of coefficients do not match with the polynomial degree ')
      fprintf( ' \n Repeat your input for coefficients! ')
   end

   while 1
      fprintf( ' \n The absciss extremes a and b (a < b)')
      fprintf( ' \n a = ');  a = input('');
      fprintf( ' b = ');  b = input('');
      if a < b;  break; end
      fprintf( ' \n a should be less than b! ')
      fprintf( ' \n Repeat your input for a and b! ')
   end

   fprintf( ' \n Now the polynomial function is plotted ')
   x=a:0.1:b;
   f = polyval(coef,x);
   plot(x, f, 'b-');
   centaxes
   hold on;

   for i=1:km+1
     pow_def(i) = km - i + 1;  % diagonal elements
   end;

   dev1_coef=pow_def.*coef;
   dev1_coef=dev1_coef(1:km)
   df=polyval(dev1_coef, x)
   fprintf( ' \n Now the 1st derivative of the polynomial function is plotted \n')
   plot(x, df,'c-');  % the analytical 1st derivative
   hold on;

   fprintf( '** Give the Po initial abscissa value approximation to a zero of f:\n');
   p0 = input('Initial approximation Po = ');

   while 1
      fprintf( '\n Input the tolerance value for Po, delta = ')
      delta = input('');
      if delta < 1 & delta > 0;  break; end
      fprintf( ' \n Delta should be greater than 0 and less than 1!!! ')
      fprintf( ' \n Repeat your input for tolerance value! ')
   end

   while 1
      fprintf( '\n Input the tolerance value for the function Y values, epsilon = ')
      epsilon = input('');
      if epsilon < 1 & epsilon > 0;  break; end
      fprintf( ' \n Epsilon should be greater than 0 and less than 1!!! ')
      fprintf( ' \n Repeat your input for tolerance value! ')
   end

   while 1
      fprintf( '\n Input the maximum number of iterations, N_max = ')
      max1 = input('');
      if max1 > 1;  break; end
      fprintf( ' \n The number of iterations should be greater than 1!!! ')
      fprintf( ' \n Repeat your input for the number of iterations value! ')
   end

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

   fprintf('\n==============RESULTS================')

   fprintf('\n Po is the Newton-Raphson approximation to zero = %12.6f', p0);
   fprintf('\n Err is the error estimate for Po = %12.6f', err);
   fprintf('\n K is the number of iterations = %u', k);
   fprintf('\n Y is the function value f(Po) = %12.6f', y);

   fprintf('\n============================================')
   kont1 = input( 'Type 1 to continue, or 0 to stop: ' );
   if kont1 == 0, break; end
end
 
 

As an example we used the same polynomial function used before: f(x) = x3 + x2 - 1.25x - 0.75 . The inputs and the outputs of the algorithm captured from the  Matlab command window is given below:

============================================
 Newton-Raphson Method. Polynomial examples
** The degree of the polynom ? N = 3

 Input the polynomial coefficients
 like [x^N x^N-1 ... x^1 x^0][1 1 -1.25 -0.75]

 The absciss extremes a and b (a < b)
 a = -10
 b = 10

 Now the 1st derivative of the polynomial function is plotted
** Give the Po initial abscissa value approximation to a zero of f:
Initial approximation Po = -2.34

 Input the tolerance value for Po, delta = 0.034

 Input the tolerance value for the function Y values, epsilon = 0.023

 Input the maximum number of iterations, N_max = 10

==============RESULTS================
 Po is the Newton-Raphson approximation to zero = -1.500155
 Err is the error estimate for Po = 0.010475
 K is the number of iterations = 4
 Y is the function value f(Po) = -0.000386
============================================Type 1 to continue, or 0 to stop:

The figure window depicts the function f(x) and its first derivative:


fzero_root.m

This application is simple graphical comparison of the root approximated by means of the cursor and the exact one. The code given below uses the plotapp Matlab function to give the user the power to approximate graphically the root. Plotapp plots a function and allows the user to approximate a particular root using the cursor. Fzero Matlab function is the one used to determine the exact root. Below is depicted the figure window and the command window for the example executed for the function f306.m

============================================
Approximate root is     0.75
Exact root is  0.73906
============================================Type 1 to continue, or 0 to stop:
 

mat_root.m

This simple application uses the fzero Matlab function to calculate the exact root of the f306.m function. Students are encouraged to executed these algorithms writting their own function in the f306.m file.
 

Other iterative method to approximate the root of a function is the Bairstow's method. It is beyond the scope of this course but is given though in bairstow.m.


Valentin Muresan, Dublin City University, muresanv@eeng.dcu.ie