Home > K25104 > NonlinearSystems > regula_falsi.m

regula_falsi

PURPOSE ^

Implements regula falsi algorithm to find root of f(x)=0

SYNOPSIS ^

function [ x,k ] = regula_falsi( f,a,b,tol,max )

DESCRIPTION ^

 Implements regula falsi algorithm to find root of f(x)=0
 Input arguments:
   f, function handle or vector of polynomial coeffiencts
   a,b, end points of initial search interval
   tol, tolerance
   max, maximum number of iterations
 Output arguments:
   x, solution
   k, number of iterations to convergence

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ x,k ] = regula_falsi( f,a,b,tol,max )
0002 % Implements regula falsi algorithm to find root of f(x)=0
0003 % Input arguments:
0004 %   f, function handle or vector of polynomial coeffiencts
0005 %   a,b, end points of initial search interval
0006 %   tol, tolerance
0007 %   max, maximum number of iterations
0008 % Output arguments:
0009 %   x, solution
0010 %   k, number of iterations to convergence
0011 
0012 % first check user inputs
0013 if a>=b;    
0014     error('specifed interval invalid');
0015 elseif tol<=0;
0016     error('tol must be positive');
0017 elseif max<=0;
0018     error('max must be positive');
0019 elseif isa(f,'function_handle');
0020     % do nothing
0021 elseif isa(f,'double'); 
0022     [n,m]=size(f);  % find the size of f
0023     if n~=1 && m~=1; % then f is not a vector
0024         error(['input must be a function handle or vector holding ',...
0025             'polynomial coefficients']);
0026     else 
0027         v=poly2sym(f);  % this converts the vector of coefficients to an
0028                         % expression for the polynomial
0029         f=matlabFunction(v);    % this converts the expression to a
0030                                 % function handle
0031     end
0032 else
0033     error(['f must be a function handle or vector holding ',...
0034             'polynomial coefficients']);
0035 end
0036 
0037 % start of algorithm
0038 fa=feval(f,a);  % evaluate f(a)
0039 fb=feval(f,b);  % evaluate f(b)
0040 k=0;            % initialize counter
0041 
0042 if fa==0;           % then a is a root
0043     x=a; return;
0044 elseif fb==0;       % then b is a root
0045     x=b; return;    
0046 elseif fa*fb>0;     % then f(a),f(b) have the same sign
0047     error('function values at endpoints have the same sign');
0048 end
0049 
0050 while abs(fa)>tol && abs(fb)>tol && k<=max;   
0051     
0052     m=(fb*a-fa*b)/(fb-fa);  % calculate estimate of root
0053     fm=feval(f,m);  % evaluate f(m)
0054     
0055     if fm==0;       % then m is a root
0056         x=m; return;
0057     end
0058     
0059     if fm*fa<0;     % then f(a) and f(m) have opposite signs
0060         b=m;        % update b
0061         fb=fm;      % update f(b)
0062     else            % f(m) and f(b) have opposite signs
0063         a=m;        % update a
0064         fa=fm;      % update f(a)
0065     end
0066     
0067     k=k+1;      % increment counter
0068 end
0069 
0070 if abs(fa)<tol;
0071     x=a;    % set x
0072 else
0073     x=b;    % set x
0074 end
0075 
0076 if k==max+1;
0077     disp('maximum number of iterations reached');
0078 end
0079 
0080 end
0081

Generated on Mon 18-Jan-2016 10:25:49 by m2html © 2005