Numerical Mathematics with

MATLAB

Reza Abazari
2008
Contents
1. Rootfinding for Nonlinear Equations ...…......………….………………….1-38

2. Numerical Linear Algebra ……………………………….……..…………..40-76

3. Polynomial Interpolation ………………………………….….……………78-102

4. Numerical Integration ………………………..…………….…………..…104-137

5. Nonlinear Systems and Numerical Optimization ……..….………….139-169

6. Numerical solution of ODE ………………….…………….…………….171-199

7. Numerical solution of PDE …………………..……………..……………201-235

2
 1. Rootfinding for Nonlinear Equations
function R = approot (f,X,epsilon)

% Input - f is object function

% - X is the vector of abscissas
% - epsilon is the tolerance
% Output - R is the vector of approximate locations for roots

% If f is an M-file function call R = approot (@f,X,epsilon).

% If f is an anonymous function call R = approot (f,X,epsilon).

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

Y=f(X);
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

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

% bisec_g(f_name, a,c, xmin, xmax, n_points)

% Bisection method with graphics
% a, c : end points of initial interval
% xmin and xmax : limits on the graphic plot.
% n_points: number of points used in plotting.
%
function bisec_g(f_name, a,c, xmin, xmax, n_points)
% it_limit : limit of iteration number
% Y_a, Y_c : y values of the current end points
% fun_f(x) : functional value at x

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

clg, hold off;

hold on
clear Y_a, clear Y_c
wid_x = xmax - xmin; dx = (xmax- xmin)/n_points;
xp=xmin:dx:xmax; yp=feval(f_name, xp);
plot(xp,yp); xlabel('x');ylabel('f(x)');
title('Bisection Method')
ymin=min(yp); ymax=max(yp);wid_y = ymax-ymin;

3
 yp=0.*xp; plot(xp,yp)
fprintf( 'Bisection Scheme\n\n' );
tolerance = 0.000001; it_limit = 30;
fprintf( ' It. a b c f(a) ');
fprintf( ' f(c) abs(f(c)-f(a))/2\n' );
it = 0;
Y_a = feval(f_name, a ); Y_c = feval(f_name, c );
plot([a,a],[Y_a,0]); text(a,-0.1*wid_y,'x=a')
plot([c,c],[Y_c,0]); text(c,-0.1*wid_y,'x=c')
if ( Y_a*Y_c > 0 ) fprintf( ' f(a)f(c) > 0 \n' );
else
while 1
it = it + 1;
b = (a + c)/2; Y_b = feval(f_name, b );
plot([b,b],[Y_b,0],':'); plot(b,0,'o')
if it<4, text(b, wid_y/20, [num2str(it)]), end
fprintf('%3.0f %10.6f, %10.6f', it, a, b );
fprintf('%10.6f, %10.6f, %10.6f', c, Y_a, Y_c );
fprintf( ' %12.3e\n', abs((Y_c - Y_a)/2));
if ( abs(c-a)<=tolerance )
fprintf( ' Tolerance is satisfied. \n' );break
end
if ( it>it_limit )
fprintf( 'Iteration limit exceeded.\n' ); break
end
if( Y_a*Y_b <= 0 ) c = b; Y_c = Y_b;
else a = b; Y_a = Y_b;
end
end
fprintf('Final result: Root = %12.6f \n', b );
end
x=b;
plot([x x],[0.05*wid_y 0.2*wid_y])
text( x, 0.25*wid_y, 'Final solution')
plot([x (x-wid_x*0.004)],[0.05*wid_y 0.09*wid_y])
plot([x (x+wid_x*0.004)],[0.05*wid_y 0.09*wid_y])

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

% bisec_n(f_name, a,c): bisection method without graphics

% f_name: definition of the equation to solve
% a , c : end points of initial interval
%
function bisec_n(f_name, a,c)
f_name
% tolerance : tolerance
% it_limit : limit of iteration number
% Y_a, Y_c : y values of the current end points

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

fprintf( 'Bisection Scheme\n\n' );

tolerance = 0.000001; it_limit = 30;
fprintf( ' It. a b c f(a) ');
fprintf( ' f(b) f(c)\n' );
it = 0;
Y_a = feval(f_name, a ); Y_c = feval(f_name, c );
if ( Y_a*Y_c > 0 )
fprintf( '\n\n Stopped because f(a)f(c) > 0 \n' );
fprintf( '\n Change a or b and run again.\n' );

4
 else
while 1
it = it + 1;
b = (a + c)/2; Y_b = feval(f_name, b );
fprintf('%3.0f %10.6f, %10.6f', it, a, b );
fprintf('%10.6f, %10.6f, %10.6f, %10.6f\n', c, Y_a, Y_b, Y_c );
if ( abs(c-a)<=tolerance )
fprintf( ' Tolerance is satisfied. \n' );break
end
if ( it>it_limit )
fprintf( 'Iteration limit exceeded.\n' ); break
end
if( Y_a*Y_b <= 0 ) c = b; Y_c = Y_b;
else a = b; Y_a = Y_b;
end
end
fprintf('Final result: Root = %12.6f \n', b );
end.

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [c,err,yc]=bisect(f,a,b,delta)
%Input - f is the function
% - a and b are the left and right endpoints
% - delta is the tolerance
%Output - c is the zero
% - yc= f(c)
% - err is the error estimate for c
%If f is defined as an M-file function use the @ notation
% call [c,err,yc]=bisect(@f,a,b,delta).
%If f is defined as an anonymous function use the
% call [c,err,yc]=bisect(f,a,b,delta).

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

ya=f(a);
yb=f(b);
if ya*yb > 0,return,end
max1=1+round((log(b-a)-log(delta))/log(2));
for k=1:max1
c=(a+b)/2;
yc=f(c);
if yc==0
a=c;
b=c;
elseif yb*yc>0
b=c;
yb=yc;
else
a=c;
ya=yc;
end
if b-a < delta, break,end
end

c=(a+b)/2;
err=abs(b-a);
yc=f(c);

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

5
 function [k,p,err,P] = fixpt(g,p0,tol,max1)

%Input - g is the iteration function

% - p0 is the initial guess for the fixed-point
% - tol is the tolerance
% - max1 is the maximum number of iterations
%Output- k is the number of iterations
% - p is the approximation to the fixed-point
% - err is the error in the approximation
% - P' contains the sequence {pn}

%If g is defined as an M-file function use the @ notation

% call [k,p,err,P]=fixedpoint(@g,p0,tol,max1).
%If g is defined as an anonymous function use the
% call [k,p,err,P]=fixedpoint(g,p0,tol,max1).

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

P(1)= p0;

for k=2:max1
P(k)=g(P(k-1));
err=abs(P(k)-P(k-1));
relerr=err/(abs(P(k))+eps);
p=P(k);
if (err<tol) | (relerr<tol),break;end
end

P=P';

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [S,E,G]=golden(f,a,b,delta,epsilon)

%Input - f is the object function

% - a and b are the endpoints of the interval
% - delta is the tolerance for the abscissas
% - epsilon is the tolerance for the ordinates
%Output - S=(p,yp) contains the abscissa p and
% the ordinate yp of the minimum
% - E=(dp,dy) contains the error bounds for p and yp
% - G is an n x 4 matrix: the kth row contains [ak ck dk bk];
% the values of a, c, d, and b at the kth iteration

%If f is defined as an M-file function use the @ notation

% call [S,E,G]=golden(@f,a,b,delta,epsilon).
%If f is defined as an anonymous function use the
% call [S,E,G]=golden(f,a,b,delta,epsilon).

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

r1=(sqrt(5)-1)/2;
r2=r1^2;
h=b-a;
ya=f(a);
yb=f(b);
c=a+r2*h;
d=a+r1*h;

6
 yc=f(c);
yd=f(d);
k=1;
A(k)=a; B(k)=b; C(k)=c; D(k)=d;

while (abs(yb-ya)>epsilon)|(h>delta)
k=k+1;
if (yc<yd)
b=d;
yb=yd;
d=c;
yd=yc;
h=b-a;
c=a+r2*h;
yc=f(c);
else
a=c;
ya=yc;
c=d;
yc=yd;
h=b-a;
d=a+r1*h;
yd=f(d);
end
A(k)=a; B(k)=b; C(k)=c; D(k)=d;
end

dp=abs(b-a);
dy=abs(yb-ya);
p=a;
yp=ya;

if (yb<ya)
p=b;
yp=yb;
end

G=[A' C' D' B'];

S=[p yp];
E=[dp dy];

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [P0,y0,err,P]=grads(F,G,P0,max1,delta,epsilon,show)

%Input - F is the object function input as an M-file function

% - G = -(1/norm(grad F))*grad F; the search direction
% input as an M-file function
% - P0 is the initial starting point
% - max1 is the maximum number of iterations
% - delta is the tolerance for hmin in the single parameter
% minimization in the search direction
% - epsilon is the tolerance for the error in y0
% - show; if show==1 the iterations are displayed
%Output - P0 is the point for the minimum.
% - y0 is the function value F(P0)
% - err is the error bound for y0
% - P is a vector containing the iterations

%Use @F and @G to call the M-file functions F and G.

7
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

if nargin==5, show = 0; end

[mm n] = size(P0);
maxj = 10; big = 1e8; h = 1;
P=zeros(maxj,n+1);
len = norm(P0);
y0 = F(P0);
if (len>1e4), h = len/1e4; end
err = 1;cnt = 0;cond = 0;
P(cnt+1,:)=[P0 y0];

while (cnt<max1 & cond~=5 & (h>delta | err>epsilon))

%Compute search direction

S = G(P0);

%Start single parameter quadratic minimization

P1 = P0 + h*S;
P2 = P0 + 2*h*S;
y1 = F(P1);
y2 = F(P2);
cond = 0; j = 0;
while (j<maxj & cond==0)
len = norm(P0);
if (y0<y1)
P2 = P1;
y2 = y1;
h = h/2;
P1 = P0 + h*S;
y1 = F(P1);
else
if (y2<y1)
P1 = P2;
y1 = y2;
h = 2*h;
P2 = P0 + 2*h*S;
y2 = F(P2);
else
cond = -1;
end
end
j = j+1;
if (h<delta), cond=1; end
if (abs(h)>big | len>big), cond=5; end
end
if (cond==5)
Pmin = P1;
ymin = y1;
else
d = 4*y1 - 2*y0 - 2*y2;
if (d<0)
hmin = h*(4*y1-3*y0-y2)/d;
else
cond = 4;
hmin = h/3;
end

%Construct the next point

Pmin = P0 + hmin*S;

8
 ymin =F(Pmin);

%Determine magnitude of next h

h0 = abs(hmin);
h1 = abs(hmin-h);
h2 = abs(hmin-2*h);
if (h0<h), h = h0; end
if (h1<h), h = h1; end
if (h2<h), h = h2; end
if (h==0), h = hmin; end
if (h<delta), cond=1; end

%Termination test for minimization

e0 = abs(y0-ymin);
e1 = abs(y1-ymin);
e2 = abs(y2-ymin);
if (e0~=0 & e0<err), err = e0; end
if (e1~=0 & e1<err), err = e1; end
if (e2~=0 & e2<err), err = e2; end
if (e0==0 & e1==0 & e2==0), err = 0; end
if (err<epsilon), cond=2; end
if (cond==2 & h<delta), cond=3; end
end
cnt = cnt+1;
P(cnt+1,:)=[Pmin ymin];
P0 = Pmin;
y0 = ymin;
end

if (show==1)
disp(P);
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

% Abstract : This simulation illustrates the convergence of the

% bisection method of finding roots of an equation f(x)=0.
%
% INPUTS: Enter the following
% Function in f(x)=0
f = inline('x^3-0.165*x^2+3.993*10^(-4)');
% Lower initial guess
xl = 0.0;
% Upper intial guess
xu = 0.11;
% Lower bound of range of 'x' to be seen
lxrange = -0.02;
% Upper bound of range of 'x' to be seen
uxrange = 0.12;
% Maximum number of iterations
nmax = 30;
% Guess of root for Matlab root function
xrguess=0.05;
%
% SOLUTION

% The following finds the upper and lower 'y' limits for the plot based on the
% given

9
% 'x' range in the input section.
maxi = f(lxrange);
mini = f(lxrange);
for i=lxrange:(uxrange-lxrange)/10:uxrange
if f(i) > maxi
maxi = f(i);
end
if f(i) < mini
mini = f(i);
end
end
tot=maxi-mini;
mini=mini-0.1*tot;
maxi=maxi+0.1*tot;

% This calculates window size to be used in figures

set(0,'Units','pixels')
scnsize = get(0,'ScreenSize');
wid = round(scnsize(3));
hei = round(0.95*scnsize(4));
wind = [1, 1, wid, hei];

% This graphs the function and two lines representing the two guesses
figure('Position',wind)
clf
subplot(2,1,2),fplot(f,[lxrange,uxrange])
hold on
plot([xl,xl],[maxi,mini],'y','linewidth',2)
plot([xu,xu],[maxi,mini],'g','linewidth',2)
plot([lxrange,uxrange],[0,0],'k','linewidth',1)
title('Entered function on given interval with upper and lower guesses')

% This portion adds the text and math

subplot(2,1,1), text(0,0.8,['Check first if the lower and upper guesses bracket
the root of the equation'])
axis off
text(0.2,0.6,['f(xl) = ',num2str(f(xl))])
text(0.2,0.4,['f(xu) = ',num2str(f(xu))])
text(0.2,0.2,['f(xu)*f(xl) = ',num2str(f(xu)*f(xl))])
hold off

% True solution
% This is how Matlab calculates the root of a non-linear equation.
xrtrue=fzero(f,xrguess);

% Value of root by iterations

% Here the bisection method algorithm is applied to generate the values of the
roots, true error,
% absolute relative true error, approximate error, absolute relative
approximate error, and the
% number of significant digits at least correct in the estimated root as a
function of number of
% iterations.
for i=1:nmax
xr(i)=(xl+xu)/2;
if f(xu)*f(xr(i))<0
xl=xr(i);
else
xu=xr(i);
end
end

10
 n=1:nmax;

% Absolute true error

for i=1:nmax
Et(i)=abs(xrtrue-xr(i));
end

% Absolute relative true error

for i=1:nmax
et(i)=abs(Et(i)/xrtrue)*100;
end

% Absolute approximate error

for i=1:nmax
if i==1
Ea(i)=0;
else
Ea(i)=abs(xr(i)-xr(i-1));
end
end

% Absolute relative approximate error

for i=1:nmax
ea(i)=abs(Ea(i)/xrtrue)*100;
end

% Significant digits at least correct

for i=1:nmax
if (ea(i)>5 | i==1)
sigdigits(i)=0;
else
sigdigits(i)=floor((2-log10(ea(i)/0.5)));
end
end

% The graphs
figure('Position',wind)
plot(n,xr,'r','linewidth',2)
title('Estimated root as a function of number of iterations')
figure('Position',wind)
subplot(2,1,1), plot(n,Et,'b','linewidth',2)
title('Absolute true error as a function of number of iterations')
subplot(2,1,2), plot(n,et,'b','linewidth',2)
title('Absolute relative true error as a function of number of iterations')
figure('Position',wind)
subplot(2,1,1), plot(n,Ea,'g','linewidth',2)
title('Absolute relative error as a function of number of iterations')
subplot(2,1,2), plot(n,ea,'g','linewidth',2)
title('Absolute relative approximate error as a function of number of
iterations')
figure('Position',wind)
bar(sigdigits,'r')

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

% Abstract : This simulation shows how the bisection method works

% for finding roots of an equation f(x)=0

11
clear all
% INPUTS: Enter the following
% Function in f(x)=0
f = inline('x^3-0.165*x^2+3.993*10^(-4)');
% Upper initial guess
xu = 0.11;
% Lower initial guess
xl = 0.0;
% Lower bound of range of 'x' to be seen
lxrange = -0.02;
% Upper bound of range of 'x' to be seen
uxrange = 0.12;
%
% SOLUTION

% The following finds the upper and lower 'y' limits for the plot
% based on the given
% 'x' range in the input section.
maxi = f(lxrange);
mini = f(lxrange);
for i=lxrange:(uxrange-lxrange)/10:uxrange
if f(i) > maxi
maxi = f(i);
end
if f(i) < mini
mini = f(i);
end
end
tot=maxi-mini;
mini=mini-0.1*tot;
maxi=maxi+0.1*tot;

% This calculates window size to be used in figures

set(0,'Units','pixels')
scnsize = get(0,'ScreenSize');
wid = round(scnsize(3));
hei = round(0.95*scnsize(4));
wind = [1, 1, wid, hei];

% This graphs the function and two lines representing the two guesses
figure('Position',wind)
clf
fplot(f,[lxrange,uxrange])
hold on
plot([xl,xl],[maxi,mini],'y','linewidth',2)
plot([xu,xu],[maxi,mini],'g','linewidth',2)
plot([lxrange,uxrange],[0,0],'k','linewidth',1)
title('Entered function on given interval with initial guesses')

hold off

%
% Iteration 1
figure('Position',wind)
xr=(xu+xl)/2;
% This graphs the function and two lines representing the two guesses
clf
subplot(3,1,2),fplot(f,[lxrange,uxrange])
hold on
plot([xl,xl],[maxi,mini],'y','linewidth',2)
plot([xu,xu],[maxi,mini],'g','linewidth',2)

12
plot([xr,xr],[maxi,mini],'r','linewidth',2)
plot([lxrange,uxrange],[0,0],'k','linewidth',1)
title('Entered function on given interval with upper and lower guesses')

% This portion adds the text and math to the top part of the figure
% window
subplot(3,1,1), text(0,1,['Iteration 1'])
text(0.2,.8,['xr = (xu + xl)/2 = ',num2str(xr)])
text(0,.6,['Finding the value of the function at the lower and upper… guesses
and the estimated root'])
text(0.2,.4,['f(xl) = ',num2str(f(xl))])
text(0.2,.2,['f(xu) = ',num2str(f(xu))])
text(0.2,0,['f(xr) = ',num2str(f(xr))])
axis off
hold off

% Check the interval between which the root lies

if f(xu)*f(xr)<0
xl=xr;
else
xu=xr;
end

% This portion adds the text and math to the bottom part of the
% figure window
subplot(3,1,3), text(0,1,['Check the interval between which the root… lies. Does
it lie in ( xl , xu ) or ( xr , xu )?'])
text(0,.8,[''])
text(0.2,0.6,['xu = ',num2str(xu)])
text(0.2,0.4,['xl = ',num2str(xl)])
axis off
xp=xr;

%
% Iteration 2
figure('Position',wind)
xr=(xu+xl)/2;
% This graphs the function and two lines representing the two guesses
clf
subplot(3,1,2),fplot(f,[lxrange,uxrange])
hold on
plot([xl,xl],[maxi,mini],'y','linewidth',2)
plot([xu,xu],[maxi,mini],'g','linewidth',2)
plot([xr,xr],[maxi,mini],'r','linewidth',2)
plot([lxrange,uxrange],[0,0],'k','linewidth',1)
title('Entered function on given interval with upper and lower… guesses')

% This portion adds the text and math to the top part of the figure
% window
subplot(3,1,1), text(0,1,['Iteration 2'])
text(0.2,.8,['xr = (xu + xl) / 2 = ',num2str(xr)])
text(0,.6,['Finding the value of the function at the lower and upper… guesses
and the estimated root'])
text(0.2,.4,['f(xl) = ',num2str(f(xl))])
text(0.2,.2,['f(xu) = ',num2str(f(xu))])
text(0.2,0,['f(xr) = ',num2str(f(xr))])
axis off
hold off

% Check the interval between which the root lies

if f(xu)*f(xr)<0

13
 xl=xr;
else
xu=xr;
end

% Calculate relative approximate error

ea=abs((xr-xp)/xr)*100;

% This portion adds the text and math to the bottom part of the
% figure window
subplot(3,1,3), text(0,1,['Absolute relative approximate error'])
text(0,.8,['ea = abs((xr - xp) / xr)*100 = ',num2str(ea),'%'])
text(0,.4,['Check the interval between which the root lies. Does it… lie in ( xl
, xu ) or ( xr , xu )?'])
text(0.2,0.2,['xu = ',num2str(xu)])
text(0.2,0,['xl = ',num2str(xl)])
axis off
xp=xr;

%
% Iteration 3
figure('Position',wind)
xr=(xu+xl)/2;
% This graphs the function and two lines representing the two guesses
clf
subplot(3,1,2),fplot(f,[lxrange,uxrange])
hold on
plot([xl,xl],[maxi,mini],'y','linewidth',2)
plot([xu,xu],[maxi,mini],'g','linewidth',2)
plot([xr,xr],[maxi,mini],'r','linewidth',2)
plot([lxrange,uxrange],[0,0],'k','linewidth',1)
title('Entered function on given interval with upper and lower … guesses')

% This portion adds the text and math to the top part of the figure
% window
subplot(3,1,1), text(0,1,['Iteration 3'])
text(0.2,.8,['xr = (xu + xl) / 2 = ',num2str(xr)])
text(0,.6,['Finding the value of the function at the lower and upper… guesses
and the estimated root'])
text(0.2,.4,['f(xl) = ',num2str(f(xl))])
text(0.2,.2,['f(xu) = ',num2str(f(xu))])
text(0.2,0,['f(xr) = ',num2str(f(xr))])
axis off
hold off

% Check the interval between which the root lies

if f(xu)*f(xr)<0
xl=xr;
else
xu=xr;
end

% Calculate relative approximate error

ea=abs((xr-xp)/xr)*100;

% This portion adds the text and math to the bottom part of the
% figure window
subplot(3,1,3), text(0,1,['Absolute relative approximate error'])
text(0,.8,['ea = abs((xr - xp) / xr)*100 = ',num2str(ea),'%'])
text(0,.4,['Check the interval between which the root lies. Does it… lie in ( xl
, xu ) or ( xr , xu )?'])

14
text(0.2,0.2,['xu = ',num2str(xu)])
text(0.2,0,['xl = ',num2str(xl)])
axis off
xp=xr;

%
% Iteration 4
figure('Position',wind)
xr=(xu+xl)/2;
% This graphs the function and two lines representing the two guesses
clf
subplot(3,1,2),fplot(f,[lxrange,uxrange])
hold on
plot([xl,xl],[maxi,mini],'y','linewidth',2)
plot([xu,xu],[maxi,mini],'g','linewidth',2)
plot([xr,xr],[maxi,mini],'r','linewidth',2)
plot([lxrange,uxrange],[0,0],'k','linewidth',1)
title('Entered function on given interval with upper and lower guesses')

% This portion adds the text and math to the top part of the figure window
subplot(3,1,1), text(0,1,['Iteration 4'])
text(0.2,.8,['xr = (xu + xl) / 2 = ',num2str(xr)])
text(0,.6,['Finding the value of the function at the lower and upper guesses and
the estimated root'])
text(0.2,.4,['f(xl) = ',num2str(f(xl))])
text(0.2,.2,['f(xu) = ',num2str(f(xu))])
text(0.2,0,['f(xr) = ',num2str(f(xr))])
axis off
hold off

% Check the interval between which the root lies

if f(xu)*f(xr)<0
xl=xr;
else
xu=xr;
end

% Calculate relative approximate error

ea=abs((xr-xp)/xr)*100;

% This portion adds the text and math to the bottom part of the figure window
subplot(3,1,3), text(0,1,['Absolute relative approximate error'])
text(0,.8,['ea = abs((xr - xp) / xr)*100 = ',num2str(ea),'%'])
text(0,.4,['Check the interval between which the root lies. Does it lie in ( xl
, xu ) or ( xr , xu )?'])
text(0.2,0.2,['xu = ',num2str(xu)])
text(0.2,0,['xl = ',num2str(xl)])
axis off
xp=xr;

%
% Iteration 5
figure('Position',wind)
xr=(xu+xl)/2;
% This graphs the function and two lines representing the two guesses
clf
subplot(3,1,2),fplot(f,[lxrange,uxrange])
hold on
plot([xl,xl],[maxi,mini],'y','linewidth',2)
plot([xu,xu],[maxi,mini],'g','linewidth',2)
plot([xr,xr],[maxi,mini],'r','linewidth',2)

15
 plot([lxrange,uxrange],[0,0],'k','linewidth',1)
title('Entered function on given interval with upper and lower guesses')

% This portion adds the text and math to the top part of the figure window
subplot(3,1,1), text(0,1,['Iteration 5'])
text(0.2,.8,['xr = (xu + xl) / 2 = ',num2str(xr)])
text(0,.6,['Finding the value of the function at the lower and upper guesses and
the estimated root'])
text(0.2,.4,['f(xl) = ',num2str(f(xl))])
text(0.2,.2,['f(xu) = ',num2str(f(xu))])
text(0.2,0,['f(xr) = ',num2str(f(xr))])
axis off
hold off

% Check the interval between which the root lies

if f(xu)*f(xr)<0
xl=xr;
else
xu=xr;
end

% Calculate relative approximate error

ea=abs((xr-xp)/xr)*100;

% This portion adds the text and math to the bottom part of the figure window
subplot(3,1,3), text(0,1,['Absolute relative approximate error'])
text(0,.8,['ea = abs((xr - xp) / xr)*100 = ',num2str(ea),'%'])
text(0,.4,['Check the interval between which the root lies. Does it lie in ( xl
, xu ) or ( xr , xu )?'])
text(0.2,0.2,['xu = ',num2str(xu)])
text(0.2,0,['xl = ',num2str(xl)])
axis off
xp=xr;

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [p,y,err]=muller(f,p0,p1,p2,delta,epsilon,max1)

%Input - f is the object function

% - p0, p1, and p2 are the initial approximations
% - delta is the tolerance for p0, p1, and p2
% - epsilon the the tolerance for the function values y
% - max1 is the maximum number of iterations
%Output - p is the Muller approximation to the zero of f
% - y is the function value y = f(p)
% - err is the error in the approximation of p.

%If f is defined as an M-file function use the @ notation

% call [p,y,err]=muller(@f,p0,p1,p2,delta,epsilon,max1).
%If f is defined as an anonymous function use the
% call [p,y,err]=muller(f,p0,p1,p2,delta,epsilon,max1).

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

%Initalize the matrices P and Y

P=[p0 p1 p2];
Y=f(P);

%Calculate a and b in formula (15)

16
 for k=1:max1
h0=P(1)-P(3);h1=P(2)-P(3);e0=Y(1)-Y(3);e1=Y(2)-Y(3);c=Y(3);
denom=h1*h0^2-h0*h1^2;
a=(e0*h1-e1*h0)/denom;
b=(e1*h0^2-e0*h1^2)/denom;

%Suppress any complex roots

if b^2-4*a*c > 0
disc=sqrt(b^2-4*a*c);
else
disc=0;
end

%Find the smallest root of (17)

if b < 0
disc=-disc;
end

z=-2*c/(b+disc);
p=P(3)+z;
%Sort the entries of P to find the two closest to p
if abs(p-P(2))<abs(p-P(1))
Q=[P(2) P(1) P(3)];
P=Q;
Y=f(P);
end
if abs(p-P(3))<abs(p-P(2))
R=[P(1) P(3) P(2)];
P=R;
Y=f(P);
end

%Replace the entry of P that was farthest from p with p

P(3)=p;
Y(3) = f(P(3));
y=Y(3);

%Determine stopping criteria

err=abs(z);
relerr=err/(abs(p)+delta);
if (err<delta)|(relerr<delta)|(abs(y)<epsilon)
break
end
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [c,err,yc]=regula(f,a,b,delta,epsilon,max1)

%Input - f is the function

% - a and b are the left and right endpoints
% - delta is the tolerance for the zero
% - epsilon is the tolerance for the value of f at the zero
% - max1 is the maximum number of iterations
%Output - c is the zero
% - yc=f(c)
% - err is the error estimate for c

%If f is defined as an M-file function use the @ notation

% call [c,err,yc]=regula(@f,a,b,delta,epsilon,max1)
%If f is defined as an anonymous function use the

17
 % call [c,err,yc]=regula(f,a,b,delta,epsilon,max1)

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

ya=f(a);
yb=f(b);

if ya*yb>0
disp('Note: f(a)*f(b) >0'),
return,
end

for k=1:max1
dx=yb*(b-a)/(yb-ya);
c=b-dx;
ac=c-a;
yc=f(c);
if yc==0,break;
elseif yb*yc>0
b=c;
yb=yc;
else
a=c;
ya=yc;
end
dx=min(abs(dx),ac);
if abs(dx)<delta,break,end
if abs(yc)<epsilon, break,end
end

c;
err=abs(b-a)/2;
yc=f(c);

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [p1,err,k,y]=secant(f,p0,p1,delta,epsilon,max1)

%Input - f is the object function

% - 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)

%If f is defined as an M-file function use the @ notation

% call [c,err,yc]=bisect(@f,a,b,delta).
%If f is defined as an anonymous function use the
% call [c,err,yc]=bisect(f,a,b,delta).

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

for k=1:max1
p2=p1-f(p1)*(p1-p0)/(f(p1)-f(p0));
err=abs(p2-p1);
relerr=2*err/(abs(p2)+delta);

18
 p0=p1;
p1=p2;
y=f(p1);
if (err<delta)|(relerr<delta)|(abs(y)<epsilon),break,end
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [p0,err,k,y]=newton(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).

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

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 [V0,y0,dV,dy]=nelder(F,V,min1,max1,epsilon,show)

%Input - F is input as an M-file function

% - V is a 3xn matrix containing starting simplex
% - min1 & max1 are minimum and maximum number of iterations
% - epsilon is the tolerance
% - show == 1 displays iterations (P and Q)

%Output - V0 is the vertex forthe minimum

% - y0 is the function value F(V0)
% - dV is the size of the final simplex
% - dy is the error bound for the minimum
% - P is a matrix containing the vertex iterations
% - Q is an array containing the iterations for F(P)

%Call [V0,y0,dV,dy]=nelder(@F,V,min1,max1,epsilon,show)

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

19
if nargin==5,
show=0;
end

[mm n]=size(V);

% Order the vertices

for j=1:n+1
Z=V(j,1:n);
Y(j)=F(Z);
end

[mm lo]=min(Y);
[mm hi]=max(Y);
li=hi;
ho=lo;

for j=1:n+1
if(j~=lo&j~=hi&Y(j)<=Y(li))
li=j;
end
if (j~=hi&j~=lo&Y(j)>=Y(ho))
ho=j;
end
end

cnt=0;

% Start of Nelder-Mead algorithm

while (Y(hi)>Y(lo)+epsilon&cnt<max1)|cnt<min1
S=zeros(1:n);
for j=1:n+1
S=S+V(j,1:n);
end
M=(S-V(hi,1:n))/n;
R=2*M-V(hi,1:n);
yR=F(R);
if (yR<Y(ho))
if (Y(li)<yR)
V(hi,1:n)=R;
Y(hi)=yR;
else
E=2*R-M;
yE=F(E);
if (yE<Y(li))
V(hi,1:n)=E;
Y(hi)=yE;
else
V(hi,1:n)=R;
Y(hi)=yR;
end
end
else
if (yR<Y(hi))
V(hi,1:n)=R;
Y(hi)=yR;
end
C=(V(hi,1:n)+M)/2;
yC=F(C);
C2=(M+R)/2;

20
 yC2=F(C2);
if (yC2<yC)
C=C2;
yC=yC2;
end
if (yC<Y(hi))
V(hi,1:n)=C;
Y(hi)=yC;
else
for j=1:n+1
if (j~=lo)
V(j,1:n)=(V(j,1:n)+V(lo,1:n))/2;
Z=V(j,1:n);
Y(j)=F(Z);
end
end
end
end
[mm lo]=min(Y);
[mm hi]=max(Y);
li=hi;
ho=lo;
for j=1:n+1
if (j~=lo&j~=hi&Y(j)<=Y(li))
li=j;
end
if (j~=hi&j~=lo&Y(j)>=Y(ho))
ho=j;
end
end
cnt=cnt+1;
P(cnt,:)=V(lo,:);
Q(cnt)=Y(lo);
end

% End of Nelder-Mead algorithm

%Determine size of simplex

snorm=0;
for j=1:n+1
s=norm(V(j)-V(lo));
if(s>=snorm)
snorm=s;
end
end

Q=Q';
V0=V(lo,1:n);
y0=Y(lo);
dV=snorm;
dy=abs(Y(hi)-Y(lo));
if (show==1)
disp([P Q])
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

% this m-file calculates the real roots of the given polynomial using
% newton raphson [Link] m-file calls the functions in the two
% m-files named as syn_division.m and derivate.m.
% coeff_function is the 1xn matrix conatining the coeff of the

21
% polynomial.

% Copyright by Reza Abazari 2008-12-09

% Email : abazari-r@[Link]

function [final_roots,functionvalue] =
newton(coeff_function,initial_guess,error_tolerance,max_iterations)
iterations=0;
max_no_roots=size(coeff_function,2);
final_roots=0;
functionvalue=0;
for no_roots=1:max_no_roots-1
fun_root_new=initial_guess;
flag=1;
coeff_der_function=derivate(coeff_function);
order_fun=size(coeff_function,2);
order_der_fun=size(coeff_der_function,2);
while flag==1
fun_root_old=fun_root_new;
fx=0;
dfx=0;
nonzero=1;
while nonzero==1
powers=order_fun-1;
for index=1:order_fun
fx=fx+coeff_function(index)*fun_root_old^powers;
powers=powers-1;
end
powers=order_der_fun-1;
for index=1:order_der_fun
dfx=dfx+coeff_der_function(index)*fun_root_old^powers;
powers=powers-1;
end
if dfx==0
fun_root_old=fun_root_old+1;
else
nonzero=0;
end
end
iterations = iterations + 1;
fun_root_new = fun_root_old - fx/dfx;
if iterations >= max_iterations
flag=0;
elseif abs(fun_root_new-fun_root_old)<=error_tolerance
flag=0;
final_roots(no_roots)=fun_root_new;
functionvalue(no_roots)=fx;
end
end
coeff_function=syn_division(coeff_function,fun_root_new);
end

++++++++++++++++++++++++++++++++++++++++++

% This m-file calculates the derivative of the function, the

% limitation of
% this function is, it can calculate only the derivatives of
% power(x,n)....

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

22
 function coeff_derivative=derivate(coeff_function)
der_order=size((coeff_function),2)-1;
coeff_derivative=0;
for index=1:size((coeff_function),2)-1
coeff_derivative(index)=der_order*coeff_function(index);
der_order=der_order-1;
end

++++++++++++++++++++++++++++++++++++++++++

% This m-file takes care of synthetic division.

% By giving one polynomial and one root this function returns
% the polynomial formed with the other roots of the given polynomial
% excluding the given root.

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

function coeff_second=syn_division(coeff_function,fun_root_new)
order_fun=size((coeff_function),2);
coeff_second=0;
for index=1:size((coeff_function),2)-1
if index==1
coeff_second(index)=coeff_function(index);
else

coeff_second(index)=coeff_function(index)+fun_root_new*coeff_second(index-1);
end
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function A=remez(fun, fun_der,interval,order)

%This [Link] used findzero.m and err.m that’s given the end of this [Link].
% findzero.m

% Copyright by Reza Abazari 2008-12-09

% Email : abazari-r@[Link]

powers=ones(order+2,1)*([0:order]);% the powers of the polynomial

% repeated in rows (order +2)
% times
coeff_E =(-1).^[1:order+2];
coeff_E=coeff_E(:); % the coefficients of the E as a column array
t=1:order;
t=t(:); % the powers of the polynomial starting from 1 in a column.
% This is used when differntiation the polynomial

y=linspace(interval(1),interval(2),order+2); % the first choice of

% the (order+2) points

for i=1:10
y=y(:); % make the points array a column array
h=(y-interval(1))*ones(1,order+1); % repeat the points column
% minus the start of the
% interval (order +1) times

coeff_h=h.^powers; % raise the h matrix by the

% power matrix elementwise

23
 M=[coeff_h coeff_E]; % the matrix of the LHS of the
% linear system of equations

N= feval(fun,y); % the column vector of the RHS

% of the linear system of equations

A=M\N; % solution of the linear system of equations,

% first (order +1) element are the coefficients
% of the polynomial. Last element is the value
% of the error at these points

A1=A(1:end-1); % the coefficients only

A_der=A(2:end-1).*t; % the coeffcients of the derivative
% of the polynomial

z(1)=interval(1); % z(1) is the start point of the interval

z(order+3)=interval(2); % z(order+3) is the end point of the
% interval
% in between we fill in with the roots of the error function
for k=1: order+1
z(k+1)=findzero(@err,y(k),y(k+1),fun,A1,interval(1));
end

% between every two points in the array z, we seek the point that
% maximizes the magnitude of the error function. If there is an extreme
% point (local maximum or local minimum) between such two points of the
% z array then the derivative of the error function is zero at this
% extreme point. We thus look for the extreme point by looking for the
% root of the derivative of the error function between these two
% points. If the extreme point doesn't exist then we check the value of
% the error function at the two current points of z and pick the one
% that gives maximum magnitude.

for k=1:order+2
if sign(err(z(k),fun_der,A_der,interval(1) ))~=...
sign(err(z(k+1),fun_der,A_der,interval(1)))
% check for a change in sign
y1(k)=findzero(@err,z(k),z(k+1),fun_der,A_der,interval(1));
% the extreme point that we seek
v(k)=abs(err(y1(k),fun,A1,interval(1)));
% the value of the error function at the extreme point else
% if there is no change in sign therefore there is no extreme
% point and we compare the endpoints of the sub-interval
v1=abs(err(z(k),fun,A1,interval(1)));
% magnitude of the error function at the start of the sub-interval
v2=abs(err(z(k+1),fun,A1,interval(1)));
% magnitude of the error function at the end of the sub-interval
% pick the larger of the two
if v1>v2
y1(k)=z(k);
v(k)=v1;
else
y1(k)=z(k+1);
v(k)=v2;
end
end
end
[mx ind]=max(v);
% search for the point in the extreme points array that gives maximum
% magnitude for the error function
% if the difference between this point and the corressponding point
% in the old array is less than a certain threshold then quit the loop

24
 if abs(y(ind)-y1(ind)) <2^-30
break;
end
% compare it also with the following point if it is not
% the last point
if ind<length(y) & abs(y(ind+1)-y1(ind)) < 2^-30
break
end
% replace the old points with the new points
y=y1;
end

++++++++++++++++++++++++++++++++++++++++++

function y=findzero(fun,x0,x1,varargin)
% fun is the function that we need to compute its root.
% x0 and x1 are two arguments to the function such that the root that
% we seek lie between them and the function has different sign at
% these two points varargin are the other arguments of the function
% the value of the function at the first point
f0=feval(fun,x0,varargin{:});
%the value of the function at the second point
f1=feval(fun,x1,varargin{:});
% check that the sign of the function at x0 and x1 is different.
% Otherwise report an error
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

if sign(f0)==sign(f1)
error('the function at the two endpoints must be of opposite...
signs');
end

% find a closer point to the root using the method of chords. In the
% method of chords we simply connect the two points (x0, f(x0)) and
%(x1, f(x1))using a straight line and compute the intersection point
% with this line and the horizontal axis. This new point is closer to
% the desired root

x=x0 - f0 * ((x1-x0)/(f1-f0));

%evaluate the function at this new point

f=feval(fun,x,varargin{:});

% enter this root as long as the difference between the two points
% that sandwitch the desired root is larger than a certain threshold

while abs(f)>2^-52
% we keep one of the two old points that has a different sign than the
% new point and we overwrite the other old point with the new point

if sign(f)==sign(f0)
x0=x;
f0=f;
else
x1=x;
f1=f;
end
x=x0 - f0 * ((x1-x0)/(f1-f0));
f=feval(fun,x,varargin{:});
end

25
% at the end of the loop we reach the root with the desired precision and
% it is given by x
y=x;

++++++++++++++++++++++++++++++++++++++++++

function e= err(x,fun, A, first)

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

% the polynomial coefficients array , make it a column array

A=A(:);

% the argument array , make it a column array

x=x(:);

% order of the polynomial is equal to the number of coefficients

% minus one
order=length(A)-1;

% the powers out in a row and repeated for each argument to form a
% matrix for example if the order is 2 and we have 3 arguments in x
% then
% [0 1 2]
% powers= [0 1 2]
% [0 1 2]
powers=ones(length(x),1)*[0:order];

% each argument is repeated a number of times equal to the number of

% coefficients to form a row then each element of the resulting row
% is raised with the corresponding power in the powers matrix
temp=((x-first)*ones(1,order+1)).^powers;

% multiply the resulting matrix with the coefficients table in order

% to obtain a column array. Each element of the resulting array is
% equal to the polynomial evaluated at the distance between the
% corresponding argument and the start of the interval
temp=temp*A;

% the error vector is then given as the difference between the

% function evaluated at the argument array and the polynomial
% evaluated at the argument array
e=feval(fun,x)-temp;

Test:

some test for this [Link]:

%Testing Script%
% Example 1
fun=inline('exp(x)');
fun_der= inline('exp(x)');
interval=[0, 2^(-10)];
order =2;

A= remez(fun, fun_der, interval, order);

A1=A(1:end-1); % the 3 coefficients of the second order

polynomial
E=A(end); % the maximum approximation error

26
 % plotting the error of the whole interval
x=0: 2^-15:2^-10 ;
e=err(x,fun, A1, interval(1));
plot(x,e)
xlabel('x')
ylabel('e(x)=f(x)-p(x)')
title('Error function for when approximating exp(x)')

% Example 2
fun=inline('sin(x)');
fun_der= inline('cos(x)');
interval=[0, 2^(-10)];
order =2;

A= remez(fun, fun_der, interval, order);

A1=A(1:end-1); % the 3 coefficients of the second order

polynomial
E=A(end); % the maximum approximation error

% plotting the error of the whole interval

x=0: 2^-15:2^-10 ;
e=err(x,fun, A1, interval(1));
figure;
plot(x,e)
xlabel('x')
ylabel('e(x)=f(x)-p(x)')
title('Error function for when approximating sin(x)')

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function r = bisect(fun,xb,xtol,ftol,verbose)

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

% bisect Use bisection to find a root of the scalar equation f(x) = 0

%
% Synopsis: r = bisect(fun,xb)
% r = bisect(fun,xb,xtol)
% r = bisect(fun,xb,xtol,ftol)
% r = bisect(fun,xb,xtol,ftol,verbose)
%
% Input: fun = (string) name of function for which roots are sought
% xb = vector of bracket endpoints. xleft = xb(1), xright = xb(2)
% xtol = (optional) relative x tolerance. Default: xtol=5*eps
% ftol = (optional) relative f(x) tolerance. Default: ftol=5*eps
% verbose = (optional) print switch. Default: verbose=0, no printing
%
% Output: r = root (or singularity) of the function in xb(1) <= x <= xb(2)
if size(xb,1)>1, warning('Only first row of xb is used as bracket'); end
if nargin < 3, xtol = 5*eps; end
if nargin < 4, ftol = 5*eps; end
if nargin < 5, verbose = 0; end

xeps = max(xtol,5*eps); % Smallest tolerances are 5*eps

feps = max(ftol,5*eps);
a = xb(1,1); b = xb(1,2); % Use first row if xb is a matrix
xref = abs(b - a); % Use initial bracket in convergence test
fa = feval(fun,a); fb = feval(fun,b);
fref = max([abs(fa) abs(fb)]); % Use max f in convergence test

27
if sign(fa)==sign(fb) % Verify sign change in the interval
error(sprintf('Root not bracketed by [%f, %f]',a,b));
end

if verbose
fprintf('\nBisection iterations for %s.m\n',fun);
fprintf(' k xm fm\n');
end
k = 0; maxit = 50; % Current and max number of iterations
while k < maxit
k = k + 1;
dx = b - a;
xm = a + 0.5*dx; % Minimize roundoff in computing the midpoint
fm = feval(fun,xm);
if verbose, fprintf('%4d %12.4e %12.4e\n',k,xm,fm); end

if (abs(fm)/fref < feps) | (abs(dx)/xref < xeps) % True when root is found
r = xm; return;
end
if sign(fm)==sign(fa)
a = xm; fa = fm; % Root lies in interval [xm,b], replace a and fa
else
b = xm; fb = fm; % Root lies in interval [a,xm], replace b and fb
end
end
warning(sprintf('root not within tolerance after %d iterations\n',k));

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function Xb = brackPlot(fun,xmin,xmax,nx)

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

% brackPlot Find subintervals on x that contain sign changes of f(x)

%
% Synopsis: Xb = brackPlot(fun,xmin,xmax)
% Xb = brackPlot(fun,xmin,xmax,nx)
%
% Input: fun = (string) name of mfile function that evaluates f(x)
% xmin,xmax = endpoints of interval to subdivide into brackets.
% nx = (optional) number of samples along x axis used to test for
% brackets. The interval xmin <= x <= xmax is divided into
% nx-1 subintervals. Default: nx = 20.
%
% Output: Xb = two column matrix of bracket limits. Xb(k,1) is the left
% (lower x value) bracket and Xb(k,2) is the right bracket for
% the k^th potential root. If no brackets are found, Xb = [].

if nargin<4, nx=20; end

% --- Create data for a plot of f(x) on interval xmin <= x <= xmax
xp = linspace(xmin,xmax); yp = feval(fun,xp);
% --- Save data used to draw boxes that indicate brackets
ytop = max(yp); ybot = min(yp); % y coordinates of the box
ybox = 0.05*[ybot ytop ytop ybot ybot]; % around a bracket
c = [0.7 0.7 0.7]; % RGB color used to fill the box

% --- Begin search for brackets

x = linspace(xmin,xmax,nx); % Vector of potential bracket limits
f = feval(fun,x); % Vector of f(x) values at potential brackets

28
nb = 0; Xb = []; % Xb is null unless brackets are found
for k = 1:length(f)-1
if sign(f(k))~=sign(f(k+1)) % True if sign of f(x) changes in the interval
nb = nb + 1;
Xb(nb,:) = [x(k) x(k+1)]; % Save left and right ends of the bracket
hold on; fill([x(k) x(k) x(k+1) x(k+1) x(k)],ybox,c); % Add filled box
end
end
if isempty(Xb) % Free advice
warning('No brackets found. Check [xmin,xmax] or increase nx');
return; % return without drawing a plot
end
% --- Add plot here so that curve is on top of boxes used to indicate brackets
plot(xp,yp,[xmin xmax],[0 0]);
grid on; xlabel('x');
if isa(fun,'inline')
ylabel(sprintf('Roots of f(x) = %s',formula(fun))); % label is formul in f(x)
else
ylabel(sprintf('Roots of f(x) defined in %s',fun)); % label is name of m-file
end
hold off

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function y = bisection ( f, a, b, TOL )

%BISECTION bisection method for locating the root of a nonlinear function
% calling sequences:
% y = bisection ( 'f', a, b, TOL )
% bisection ( 'f', a, b, TOL )
% inputs:
% f string containing name of m-file defining function
% whose root is to be located
% a,b left and right endpoints, respectively, of interval
% known to contain root of f
% TOL absolute error convergence tolerance
% (iterations will be performed until the size of
% enclosing interval is smaller than 2*TOL)
% output:
% y approximate value of root
%
% NOTE:
% if BISECTION is called with no output arguments, the iteration
% number, the current enclosing interval and the current
% approximation to the root are displayed
%

% Copyright by Reza Abazari 2008-12-09

% Email : abazari-r@[Link]

sfa = sign(feval(f,a));
Nmax = floor ( log((b-a)/TOL) / log(2.0) ) + 1

for i = 1 : Nmax
p = ( a + b ) / 2.0;
sfp = sign(feval(f,p));

if ( nargout == 0 )
disp (sprintf('\t\t %3d \t (%.10f,%.10f) \t %.10f \n', i, a, b, p))
end
if ( (b-a)<2*TOL | fp == 0 )
if ( nargout == 1 )

29
 y = p;
end
return
elseif ( sfa * sfp < 0 )
b = p;
else
a = p;
sfa = sfp;
end
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function y = false_pos ( f, a, b, TOL, Nmax )

%FALSE_POS uses method of false position (regula falsi) to locate the root
% of a nonlinear function
%
% calling sequences:
% y = false_pos ( 'f', a, b, TOL, Nmax )
% false_pos ( 'f', a, b, TOL, Nmax )
%
% inputs:
% f string containing name of m-file defining function
% whose root is to be located
% a,b left and right endpoints, respectively, of interval
% known to contain root of f
% TOL absolute error convergence tolerance
% Nmax maximum number of iterations to be performed
%
% output:
% y approximate value of root
%
% NOTE:
% if FALSE_POS is called with no output arguments, the iteration
% number, the current enclosing interval and the current
% approximation to the root are displayed
%
% if the maximum number of iterations is exceeded, a message
% to this effect will be displayed and the most recent
% approximation will be returned in the output value
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

old = b;
fa = feval(f,a);
fb = feval(f,b);
for i = 1 : Nmax
new = b - fb * ( b - a ) / ( fb - fa );
fnew = feval(f,new);

if ( nargout == 0 )
disp(sprintf('\t\t %3d \t (%.10f,%.10f) \t %.10f \n', i, a, b, new))
end

if ( abs(new-old) < TOL )

if ( nargout == 1 )
y = new;
end

30
 return
elseif ( fa * fnew < 0 )
b = new;
fb = fnew;
else
a = new;
fa = fnew;
end
old = new;

end

disp('Maximum number of iterations exceeded')

if ( nargout == 1 ) y = new;
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function y = false_pos_aitken ( f, a, b, TOL, Nmax )

%FALSE_POS_AITKEN uses method of false position (regula falsi) with Aitken

% extrapolation to locate the root of a nonlinear function
%
% calling sequences:
% y = false_pos_aitken ( 'f', a, b, TOL, Nmax )
% false_pos_aitken ( 'f', a, b, TOL, Nmax )
%
% inputs:
% f string containing name of m-file defining function
% whose root is to be located
% a,b left and right endpoints, respectively, of interval
% known to contain root of f
% TOL absolute error convergence tolerance
% (applied to extrapolated sequence of approximations)
% Nmax maximum number of iterations to be performed
%
% output:
% y approximate value of root
%
% NOTE:
% if FALSE_POS_AITKEN is called with no output arguments,
% the iteration number, the current false position
% approximation to the root and the current extrapolated
% approximation to the root are displayed.
%
% if the maximum number of iterations is exceeded, a message
% to this effect will be displayed and the most recent
% approximation will be returned in the output value
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

old = b; phatold = b;
fa = feval(f,a);
fb = feval(f,b);
for i = 1 : Nmax
new = b - fb * ( b - a ) / ( fb - fa );
fnew = feval(f,new);

if ( i == 1 | i == 2 )

31
 if ( nargout == 0 )
disp ( sprintf ( '\t\t %3d \t %.10f \n', i, new ) )
end
else
phat = new - ( new - old ) ^ 2 / ( new - 2 * old + older );

if ( nargout == 0 )
disp ( sprintf ( '\t\t %3d \t %.10f \t %.10f \n', i, new, phat ) )
end

if ( abs(phat-phatold) < TOL )

if ( nargout == 1 )
y = phat;
end
return
else
phatold = phat;
end
end

if ( sign(fa) * sign(fnew) < 0 )

b = new; fb = fnew;
else
a = new; fa = fnew;
end
older = old; old = new;

end

disp('Maximum number of iterations exceeded')

if ( nargout == 1 )
y = phat;
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function y = fixed_point ( f, x0, TOL, Nmax )

%FIXED_POINT locate the fixed point of an arbitrary function using general

% fixed_point iteration
%
% calling sequences:
% y = fixed_point ( 'g', x0, TOL, Nmax )
% fixed_point ( 'g', x0, TOL, Nmax )
%
% inputs:
% g string containing name of m-file defining the
% iteration function
% x0 initial approximation to location of fixed point
% TOL absolute error convergence tolerance
% Nmax maximum number of iterations to be performed
%
% output:
% y approximate value of fixed point
%
% NOTE:
% if FIXED_POINT is called with no output arguments, the
% iteration number and the current approximation to the
% fixed point are displayed
%
% if the maximum number of iterations is exceeded, a message

32
% to this effect will be displayed and the most recent
% approximation will be returned in the output value
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

old = x0
for i = 1 : Nmax
new = feval(f,old);

if ( nargout == 0 )
disp ( sprintf ( '\t\t %3d \t %.10f \n', i, new ) )
end

if ( abs(new-old) < TOL )

if ( nargout == 1 )
y = new;
end
return
else
old = new;
end

end

disp('Maximum number of iterations exceeded')

if ( nargout == 1 )
y = new;
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function y = fixed_point_aitken ( f, x0, TOL, Nmax )

%FIXED_POINT_AITKEN uses general fixed point iteration with Aitken

% extrapolation to locate the root of a nonlinear function
%
% calling sequences:
% y = fixed_point_aitken ( 'g', x0, TOL, Nmax )
% fixed_point_aitken ( 'g', x0, TOL, Nmax )
%
% inputs:
% g string containing name of m-file defining the
% iteration function
% x0 initial approximation to location of fixed point
% TOL absolute error convergence tolerance
% (applied to extrapolated sequence of approximations)
% Nmax maximum number of iterations to be performed
%
% output:
% y approximate value of fixed point
%
% NOTE:
% if FIXED_POINT_AITKEN is called with no output arguments,
% the iteration number, the current functional iteration
% approximation to the fixed point and the current extrapolated
% approximation to the fixed point are displayed.
%
% if the maximum number of iterations is exceeded, a message
% to this effect will be displayed and the most recent

33
% approximation will be returned in the output value
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

old = x0; phatold = x0;

for i = 1 : Nmax
new = feval(f,old);

if ( i == 1 | i == 2 )
if ( nargout == 0 )
disp ( sprintf ( '\t\t %3d \t %.10f \n', i, new ) )
end
else
phat = new - ( new - old ) ^ 2 / ( new - 2 * old + older );

if ( nargout == 0 )
disp ( sprintf ( '\t\t %3d \t %.10f \t %.10f \n', i, new, phat ) )
end

if ( abs(phat-phatold) < TOL )

if ( nargout == 1 )
y = phat;
end
return
else
phatold = phat;
end
end

older = old;
old = new;

end

disp('Maximum number of iterations exceeded')

if ( nargout == 1 )
y = phat;
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function r = polyroots ( poly, Nmax, TOL )

%POLYROOTS locate the roots of an arbitrary polynomial using Laguerre's

% method
%
% calling sequences:
% r = polyroots ( poly, Nmax, TOL )
% polyroots ( poly, Nmax, TOL )
%
% inputs:
% poly vector containing the coefficients of the polynomial
% whose roots are to be computed. the first entry in
% the vector must be the leading coefficient of the
% polynomial, and the final entry must be the constant
% term. zero coefficients must be explicitly provided.
% Nmax maximum number of iterations to be performed to
% compute each root

34
% TOL absolute error convergence tolerance
%
% output:
% r vector containing the roots of the polynomial
%
% NOTE:
% if POLYROOTS is called with no output arguments, the
% roots of the polynomial and the number of iterations
% required to compute each root are displayed
%
% if the maximum number of iterations is exceeded, a message
% to this effect will be displayed and the roots for which
% convergence had been achieved will be returned
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

q = poly;
degree = length(poly)-1;
r = zeros ( 1, degree );
its = zeros ( 1, degree );
found = 0;

while ( found < degree - 2 )

[new, done] = laguerre ( q, 0.0, TOL, Nmax );

if ( done == 0 )
disp('Maximum number of iterations exceeded')
return
end;

if ( abs(imag(new)) == TOL )
found = found + 1;
r(found) = real(new);
its(found) = done;
q = deconv ( q, [ 1 -real(new) ] );
else
r(found+1) = new;
r(found+2) = conj(new);
its(found+1) = done;
its(found+2) = done;
found = found + 2;
q = deconv ( q, conv ( [1 -new], [1 -conj(new)] ) );
end;
q = real(q);
end;

if ( found == degree - 2 )
if ( q(2) == 0 )
r(degree-1) = sqrt ( -q(3)/q(1) );
r(degree) = -r(degree-1);
else
r(degree-1)=2*q(3)/(-q(2)-sign(q(2))*sqrt(q(2)*q(2)-4*q(1)*q(3)));
r(degree) = q(3) / ( r(degree-1) * q(1) );
end;
else
r(degree) = -q(2)/q(1);
end;

if ( nargout == 0 )

35
 disp ( [r' its'] )
end;

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [P,iter]= seidel(G,P,delta, max1)

%Input - G is the nonlinear fixed-point system

% saved as an M-file function
% - P is the initial guess at the solution
% - delta is the error bound
% - max1 is the number of iterations
%Output - P is the seidel approximation to the solution
% - iter is the number of iterations required

%Use the @ notation call [P,iter]=seidel(@G, P, delta, max1).

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

N=length(P);

for k=1:max1
X=P;

% X is the kth approximation to the solution

for j=1:N
A=G(X);
% Update the terms of X as they are calculated
X(j)=A(j);
end

err=abs(norm(X-P));
relerr=err/(norm(X)+eps);
P=X;
iter=k;
if (err<delta)|(relerr<delta)
break
end
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [P,iter,err]=newdim(F,JF,P,delta,epsilon,max1)

%Input -F is the system saved as the M-file F.m

% -JF is the Jacobian of F saved as the M-file JF.M
% -P is the inital approximation to the solution
% -delta is the tolerance for P
% -epsilon is the tolerance for F(P)
% -max1 is the maximum number of iterations
%Output -P is the approximation to the solution
% -iter is the number of iterations required
% -err is the error estimate for P

%Use the @ notation call

%[P,iter,err]=newdim(@F, @JF, P, delta, epsilon, max1).

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

36
Y=F(P);

for k=1:max1
J=JF(P);
Q=P-(J\Y')';
Z=F(Q);
err=norm(Q-P);
relerr=err/(norm(Q)+eps);
P=Q;
Y=Z;
iter=k;
if (err<delta)|(relerr<delta)|(abs(Y)<epsilon)
break
end
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [y, done] = laguerre ( f, x0, TOL, Nmax )

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

done = 0;
n = length ( f ) - 1;
fp = polyder ( f );
fp2 = polyder ( fp );

for i = 1 : Nmax
fx = polyval ( f, x0 );
fpx = polyval ( fp, x0 );
fp2x = polyval ( fp2, x0 );

if ( abs(fx) < TOL )

y = x0; done = i;
return
end;

gx = fpx / fx;
g2 = gx * gx;
hx = g2 - fp2x / fx;

disc = sqrt ( (n-1) * ( n * hx - g2 ) );

if ( abs(gx-disc) < abs(gx+disc) )
denom = gx+disc;
else
denom = gx-disc;
end
dx = n / denom;
x0 = x0 - dx;

if ( abs(dx) < TOL )

y = x0; done = i;
return
end
end
y = x0;

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

37
function y = steffensen ( f, x0, TOL, Nmax )

%STEFFENSEN locate the fixed point of an arbitrary function using

% Steffensen's method
%
% calling sequences:
% y = steffensen ( 'g', x0, TOL, Nmax )
% steffensen ( 'g', x0, TOL, Nmax )
%
% inputs:
% g string containing name of m-file defining the
% iteration function
% x0 initial approximation to location of fixed point
% TOL absolute error convergence tolerance
% Nmax maximum number of iterations to be performed
%
% output:
% y approximate value of fixed point
%
% NOTE:
% if STEFFENSEN is called with no output arguments, the
% iteration number and the current approximation to the
% fixed point are displayed
%
% if the maximum number of iterations is exceeded, a message
% to this effect will be displayed and the most recent
% approximation will be returned in the output value
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

p0 = x0; phatold = x0;

for i = 1 : Nmax
p1 = feval(f,p0);
p2 = feval(f,p1);

phat = p2 - ( p2 - p1 ) ^ 2 / ( p2 - 2 * p1 + p0 );

if ( nargout == 0 )
disp ( sprintf ( '\t\t %3d \t %.10f \n', i, phat ) )
end

if ( abs(phat-phatold) < TOL )

if ( nargout == 1 )
y = phat;
end
return
else
phatold = phat;
p0 = phat;
end
end
disp('Maximum number of iterations exceeded')
if ( nargout == 1 )
y = phat;
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

38
39
2. Numerical Linear Algebra
function [P,iter]= seidel(G,P,delta, max1)

%Input - G is the nonlinear fixed-point system

% saved as an M-file function
% - P is the initial guess at the solution
% - delta is the error bound
% - max1 is the number of iterations
%Output - P is the seidel approximation to the solution
% - iter is the number of iterations required

%Use the @ notation call [P,iter]=seidel(@G, P, delta, max1).

% Copyright by Reza Abazari 2008-04-23

% Email : abazari-r@[Link]

N=length(P);

for k=1:max1
X=P;

% X is the kth approximation to the solution

for j=1:N
A=G(X);
% Update the terms of X as they are calculated
X(j)=A(j);
end

err=abs(norm(X-P));
relerr=err/(norm(X)+eps);
P=X;
iter=k;
if (err<delta)|(relerr<delta)
break
end
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [p,Q]=steff(f,df,p0,delta,epsilon,max1)

%Input - f is the object function

% - 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 - p is the Steffensen approximation to the zero
% - Q is the matrix containing the Steffensen sequence

%If f and df are defined as M-file functions use the @ notation

% call [p,Q]=steff(@f,@df,p0,delta,epsilon,max1).
%If f and df are defined as anonymous functions use the
% call [p,Q]=steff(f,df,p0,delta,epsilon,max1).

%Initialize the matrix R

40
R=zeros(max1,3);
R(1,1)=p0;

for k=1:max1
for j=2:3

%Denominator in Newton-Raphson method is calculated

nrdenom=df(R(k,j-1));

%Conditional calculates Newton-Raphson approximations

if nrdenom==0
'division by zero in Newton-Raphson method'
break
else
R(k,j)=R(k,j-1)-f(R(k,j-1))/nrdenom;
end

%Denominator in Aitkens Acceleration process is calculated

aadenom=R(k,3)-2*R(k,2)+R(k,1);

%Conditional calculates Aitkens Acceleration approximations

if aadenom==0
'division by zero in Aitkens Acceleration'
break
else
R(k+1,1)=R(k,1)-(R(k,2)-R(k,1))^2/aadenom;
end
end

%Breaks out and ends program if division by zero occured

if (nrdenom==0)|(aadenom==0)
break
end

%Stopping criteria are evaluated; p and the matrix Q are determined

err=abs(R(k,1)-R(k+1,1));
relerr=err/(abs(R(k+1,1))+delta);
y=f(R(k+1,1));
if (err<delta)|(relerr<delta)|(y<epsilon)
p=R(k+1,1);
Q=R(1:k+1,:);
break
end
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [P,iter,err]=newdim(F,JF,P,delta,epsilon,max1)

%Input -F is the system saved as the M-file F.m

% -JF is the Jacobian of F saved as the M-file JF.M
% -P is the inital approximation to the solution
% -delta is the tolerance for P
% -epsilon is the tolerance for F(P)
% -max1 is the maximum number of iterations
%Output -P is the approximation to the solution
% -iter is the number of iterations required
% -err is the error estimate for P

%Use the @ notation call

%[P,iter,err]=newdim(@F, @JF, P, delta, epsilon, max1).

41
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

Y=F(P);

for k=1:max1
J=JF(P);
Q=P-(J\Y')';
Z=F(Q);
err=norm(Q-P);
relerr=err/(norm(Q)+eps);
P=Q;
Y=Z;
iter=k;
if (err<delta)|(relerr<delta)|(abs(Y)<epsilon)
break
end
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function x = gauss_elim ( A, b )

%GAUSS_ELIM solve the linear system Ax = b using Gaussian elimination

% with back substitution
%
% calling sequences:
% x = gauss_elim ( A, b )
% gauss_elim ( A, b )
%
% inputs:
% A coefficient matrix for linear system
% (matrix must be square)
% b right-hand side vector
%
% output:
% x solution vector (i.e., vector for which Ax = b)
%
% NOTE:
% this is intended as a demonstration routine - no pivoting
% strategy is implemented to reduce the effects of roundoff
% error
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

[nrow ncol] = size ( A );

if ( nrow ~= ncol )
disp ( 'gauss_elim error: Square coefficient matrix required' );
return;
end;
nb = length ( b );
if ( nrow ~= nb )
disp('gauss_elim error: Size of b-vector not compatible with matrix dimension')
return;
end;
x = zeros ( 1, nrow );

% Gaussian elimination

42
for i = 1 : nrow - 1
if ( A(i,i) == 0 )
t = min ( find ( A(i+1:nrow,i) ~= 0 ) + i );
if ( isempty(t) )
disp ( 'gauss_elim error: A matrix is singular' );
return
end;
temp = A(i,:); tb = b(i);
A(i,:) = A(t,:); b(i) = b(t);
A(t,:) = temp; b(t) = tb;
end;
for j = i+1 : nrow
m = -A(j,i) / A(i,i);
A(j,i) = 0;
A(j, i+1:nrow) = A(j, i+1:nrow) + m * A(i, i+1:nrow);
b(j) = b(j) + m * b(i);
end;
end;

%
% back substitution
%

x(nrow) = b(nrow) / A(nrow, nrow);

for i = nrow - 1 : -1 : 1
x(i) = ( b(i) - sum ( x(i+1:nrow) .* A(i, i+1:nrow) ) ) / A(i,i);
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function x = gauss_seidel ( A, b, xold, TOL, Nmax )

%GAUSS_SEIDEL approximate the solution of the linear system Ax = b by applying

% the Gauss-Seidel method (successive relaxation)
%
% calling sequences:
% x = gauss_seidel ( A, b, xold, TOL, Nmax )
% gauss_seidel ( A, b, xold, TOL, Nmax )
%
% inputs:
% A coefficient matrix for linear system - must be a
% square matrix
% b right-hand side vector for linear system
% xold vector containing initial guess for solution of
% linear system
% TOL convergence tolerance - applied to maximum norm of
% difference between successive approximations
% NMax maximum number of iterations to be performed
%
% output:
% x approximate solution of linear system
%
% NOTE:
% if GAUSS_SEIDEL is called with no output arguments, the
% iteration number and the current approximation to the
% solution are displayed
%
% if the maximum number of iterations is exceeded, a meesage
% to this effect will be displayed and the current approximation
% will be returned in the output value
%

43
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

n = length ( b );
[r c] = size ( A );
if ( c ~= n )
disp ( 'gauss_seidel error: matrix dimensions and vector dimension not
compatible' )
return
end;
xnew = zeros ( 1, n );

if ( nargout == 0 )
s = sprintf ( '%3d \t %10f ', 0, xold(1) );
for j = 2 : n
s = sprintf ( '%s%10f ', s, xold(j) );
end;
disp ( s );
end;

for its = 1 : Nmax

xnew(1) = ( b(1) - sum ( A(1,2:n) .* xold(2:n) ) ) / A(1,1);
for i = 2 : n-1
xnew(i) = ( b(i) - sum( A(i,1:i-1) .* xnew(1:i-1)) - ...
sum(A(i,i+1:n).*xold(i+1:n))) / A(i,i);
end;
xnew(n) = ( b(n) - sum ( A(n,1:n-1) .* xnew(1:n-1) ) ) / A(n,n);

if ( nargout == 0 )
s = sprintf ( '%3d \t %10f ', its, xnew(1) );
for j = 2 : n
s = sprintf ( '%s%10f ', s, xnew(j) );
end;
disp ( s );
end;

conv = max ( abs ( xnew - xold ) );

if ( conv < TOL )
x = xnew;
return
else
xold = xnew;
end;
end;
disp ( 'gauss_seidel error: maximum number of iterations exceeded' );
if ( nargout == 1 )
x = xnew;
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function y = broyden ( F, J, x0, TOL, Nmax )

%BROYDEN solve the system of nonlinear equations F(x) = 0 using

% Broyden's method
%
% calling sequences:
% y = broyden ( F, J, x0, TOL, Nmax )
% broyden ( F, J, x0, TOL, Nmax )
%
% inputs:

44
% F vector-valued function of a vector argument which
% defines the system of equations to be solved
% J matrix-valued function which computes the Jacobian
% associated with the function F
% x0 vector containing initial guess for solution of
% nonlinear system
% TOL convergence tolerance - applied to maximum norm of
% difference between successive approximations
% NMax maximum number of iterations to be performed
%
% output:
% y approximate solution of nonlinear system
%
%
% NOTE:
% if BROYDEN is called with no output arguments, each
% approximation to the solution is displayed
%
% if the maximum number of iterations is exceeded, a meesage
% to this effect will be displayed and the current approximation
% will be returned in the output value
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

Fold = feval(F,x0)';
Jold = feval(J,x0);
A0 = inv ( Jold );
dx = -A0 * Fold;
x0 = x0 + dx;
if ( nargout == 0 )
disp ( x0' )
end

for i = 2 : Nmax
Fnew = feval(F,x0)';
dy = Fnew - Fold;
u = A0 * dy;
v = dx' * A0;
denom = dx' * u;
A0 = A0 + ( dx - u ) * v / denom;
dx = -A0 * Fnew;
x0 = x0 + dx;

if ( nargout == 0 )
disp ( x0' )
end

if ( max(abs(dx)) < TOL )

if ( nargout == 1 )
y = x0;
end
return
else
Fold = Fnew;
end

end

disp('broyden error: Maximum number of iterations exceeded');

45
if ( nargout == 1 )
y = x0;
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function x = conj_grad ( A, b, x, TOL, Nmax )

%CONJ_GRAD approximate the solution of the linear system Ax = b by applying

% the conjugate gradient method
%
% calling sequences:
% x = conj_grad ( A, b, x, TOL, Nmax )
% conj_grad ( A, b, x, TOL, Nmax )
%
% inputs:
% A coefficient matrix for linear system - must be a
% square matrix
% b right-hand side vector for linear system - must be
% a column vector
% x column vector containing initial guess for solution of
% linear system
% TOL convergence tolerance - applied to Euclidean norm of
% residual vector
% NMax maximum number of iterations to be performed
%
% output:
% x approximate solution of linear system
%
% NOTE:
% if CONJ_GRAD is called with no output arguments, the
% iteration number and the current approximation to the
% solution are displayed
%
% if the maximum number of iterations is exceeded, a meesage
% to this effect will be displayed and the current approximation
% will be returned in the output value
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

n = length ( b );
[r c] = size ( A );
if ( c ~= n )
disp('conj_grad error: matrix dimensions and vector dimension not compatible')
return
end;

if ( nargout == 0 )
s = sprintf ( '%3d \t %10f ', 0, x(1) );
for j = 2 : n
s = sprintf ( '%s%10f ', s, x(j) );
end;
disp ( s );
end;
r = A * x - b;
delta0 = r' * r;
d = -r;

for its = 1 : Nmax

46
 h = A * d;
lambda = delta0 / ( d' * h );
x = x + lambda * d;
r = r + lambda * h;
delta1 = r' * r;

if ( nargout == 0 )
s = sprintf ( '%3d \t %10f ', its, x(1) );
for j = 2 : n
s = sprintf ( '%s%10f ', s, x(j) );
end;
disp ( s );
end;

if ( sqrt(delta1) < TOL )

return
else
alpha = delta1 / delta0;
delta0 = delta1;
d = -r + alpha * d;
end;
end;
disp ( 'conj_grad error: maximum number of iterations exceeded' );

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function x = jacobi ( A, b, xold, TOL, Nmax )

%JACOBI approximate the solution of the linear system Ax = b by applying

% the Jacobi method (simultaneous relaxation)
%
% calling sequences:
% x = jacobi ( A, b, xold, TOL, Nmax )
% jacobi ( A, b, xold, TOL, Nmax )
%
% inputs:
% A coefficient matrix for linear system - must be a
% square matrix
% b right-hand side vector for linear system
% xold vector containing initial guess for solution of
% linear system
% TOL convergence tolerance - applied to maximum norm of
% difference between successive approximations
% NMax maximum number of iterations to be performed
%
% output:
% x approximate solution of linear system
%
% NOTE:
% if JACOBI is called with no output arguments, the
% iteration number and the current approximation to the
% solution are displayed
%
% if the maximum number of iterations is exceeded, a meesage
% to this effect will be displayed and the current approximation
% will be returned in the output value
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

n = length ( b );

47
[r c] = size ( A );
if ( c ~= n )
disp('jacobi error: matrix dimensions and vector dimension not compatible')
return
end;
xnew = zeros ( 1, n );

if ( nargout == 0 )
s = sprintf ( '%3d \t %10f ', 0, xold(1) );
for j = 2 : n
s = sprintf ( '%s%10f ', s, xold(j) );
end;
disp ( s );
end;

for its = 1 : Nmax

xnew(1) = ( b(1) - sum ( A(1,2:n) .* xold(2:n) ) ) / A(1,1);
for i = 2 : n-1
xnew(i) = ( b(i) - sum ( A(i,1:i-1) .* xold(1:i-1) ) -...
sum ( A(i,i+1:n) .* xold(i+1:n) ) ) / A(i,i);
end;
xnew(n) = ( b(n) - sum ( A(n,1:n-1) .* xold(1:n-1) ) ) / A(n,n);

if ( nargout == 0 )
s = sprintf ( '%3d \t %10f ', its, xnew(1) );
for j = 2 : n
s = sprintf ( '%s%10f ', s, xnew(j) );
end;
disp ( s );
end;

conv = max ( abs ( xnew - xold ) );

if ( conv < TOL )
x = xnew;
return
else
xold = xnew;
end;
end;
disp ( 'jacobi error: maximum number of iterations exceeded' );
if ( nargout == 1 )
x = xnew;
end.

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [lu, pvt] = LUfactor ( A )

%LUFACTOR compute the LU decomposition for the matrix A

%
% calling sequence:
% [lu, pvt] = LUfactor ( A )
% input:
% A coefficient matrix for linear system
% (matrix must be square)
% outputs:
% lu matrix containing LU decomposition of input matrix
% A - the L matrix in the decomposition consists of
% 1's along the main diagonal together with the
% strictly lower triangular portion of the matrix lu;
% the U matrix in the decomposition is the upper

48
% triangular portion of the matrix lu
% pvt vector which indicates the permutation of the rows
% performed during the decomposition process
% NOTE:
% this routine performs partial pivoting during the
% decomposition process
%
% the system Ax = b can be solved by first applying LUfactor
% to the coefficient matrix A and then applying the companion
% routine, LUsolve, for each right-hand side vector b
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

[nrow ncol] = size ( A );

if ( nrow ~= ncol )
disp ( 'LUfactor error: Square coefficient matrix required' );
return;
end

%
% initialize row pointers
%

for i=1:nrow
pvt(i) = i;
end

for i = 1 : nrow - 1

%
% partial pivoting
%
t =min(find(abs(A(pvt(i:nrow),i)) == max(abs(A(pvt(i:nrow),i))))+ i-1 );
if ( t ~= i )
temp = pvt(i);
pvt(i) = pvt(t);
pvt(t) = temp;
end
%
% terminate if matrix is singular
%

if ( A(pvt(i),i) == 0 )
disp ( 'LUfactor error: coefficient matrix is singular' );
lu = A;
return
end
%
% elimination steps
%
for j = i+1 : nrow
m = -A(pvt(j),i) / A(pvt(i),i);
A(pvt(j),i) = -m;
A(pvt(j), i+1:nrow) = A(pvt(j), i+1:nrow) + m * A(pvt(i), i+1:nrow);
end
end

lu = A;

49
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function x = LUsolve ( lu, b, pvt )

%LUSOLVE perform forward and backward substitution to obtain the

% solution to the linear system Ax = b, where the LU
% decomposition of the coefficient matrix has previously
% been determined
%
% calling sequence:
% x = LUsolve ( lu, b, pvt )
% LUsolve ( lu, b, pvt )
%
% inputs:
% lu matrix containing LU decomposition of coefficient
% matrix for the linear system - the L matrix in the
% decomposition must consists of 1's along the main
% diagonal together with the strictly lower triangular
% portion of the matrix lu; the U matrix in the
% decomposition must be the upper triangular portion
% of the matrix lu
% b right-hand side vector for linear system
% pvt vector which indicates the permutation of the rows
% performed during the LU decomposition of the
% coefficient matrix
%
% output:
% x solution vector
%
% NOTE:
% the system Ax = b can be solved by first applying LUfactor
% to the coefficient matrix A and then applying LUsolve once
% for each right-hand side vector b
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

[nrow ncol] = size ( lu );

%
% forward substitution
%

z(1) = b(pvt(1));
for i = 2 : nrow
z(i) = b(pvt(i)) - sum ( z(1:i-1) .* lu(pvt(i), 1:i-1) );
end

%
% back substitution
%

x(nrow) = z(nrow) / lu(pvt(nrow), nrow);

for i = nrow - 1 : -1 : 1
x(i) = ( z(i) - sum ( x(i+1:nrow) .* lu(pvt(i), i+1:nrow) ) ) / lu(pvt(i),i);
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function y = newton_sys ( F, J, x0, TOL, Nmax )

50
%NEWTON_SYS solve the system of nonlinear equations F(x) = 0 using
% Newton's method
%
% calling sequences:
% y = newton_sys ( F, J, x0, TOL, Nmax )
% newton_sys ( F, J, x0, TOL, Nmax )
%
% inputs:
% F vector-valued function of a vector argument which
% defines the system of equations to be solved
% J matrix-valued function which computes the Jacobian
% associated with the function F
% x0 vector containing initial guess for solution of
% nonlinear system
% TOL convergence tolerance - applied to maximum norm of
% difference between successive approximations
% NMax maximum number of iterations to be performed
%
% output:
% y approximate solution of nonlinear system
%
% dependencies:
% this routine uses both LUfactor and LUsolve
%
% NOTE:
% if NEWTON_SYS is called with no output arguments, each
% approximation to the solution is displayed
%
% if the maximum number of iterations is exceeded, a meesage
% to this effect will be displayed and the current approximation
% will be returned in the output value
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

old = x0;
for i = 1 : Nmax
Fold = feval(F,old);
Jold = feval(J,old);
[lu pvt] = LUfactor ( Jold );
dx = LUsolve ( lu, -Fold, pvt );
new = old + dx;

if ( nargout == 0 )
disp ( new )
end

if ( max(abs(dx)) < TOL )

if ( nargout == 1 )
y = new;
end
return
else
old = new;
end
end

disp('newton_sys error: Maximum number of iterations exceeded');

if ( nargout == 1 )

51
 y = new;
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function y = tridiagonal ( c, a, b, r )

%TRIDIAGONAL solve a linear system with a tridiagonal coefficient matrix

%
% calling sequence:
% x = tridiagonal ( c, a, b, r )
% tridiagonal ( c, a, b, r )
%
% inputs:
% c vector containing the entries along lower diagonal
% of the coefficient matrix
% a vector containing the entries along main diagonal
% of the coefficient matrix
% b vector containing the entries along upper diagonal
% of the coefficient matrix
% r right-hand side vector
%
% output:
% x solution vector
%
% NOTE:
% the entries in the vectors c, a and b are assumed to be
% numbered as follows:
%
% | a(1) b(1) |
% | c(1) a(2) b(2) |
% | c(2) a(3) b(3) |
% | . . . |
% | . . . |
% | . . b(n-1) |
% | c(n-1) a(n) |
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

n = length ( a );

%
% factorization step
%

for i = 1 : n-1
b(i) = b(i) / a(i);
a(i+1) = a(i+1) - c(i) * b(i);
end

%
% forward substitution
%

r(1) = r(1) / a(1);

for i = 2 : n
r(i) = ( r(i) - c(i-1) * r(i-1) ) / a(i);
end

52
%
% back substitution
%

for i = n-1 : -1 : 1
r(i) = r(i) - r(i+1) * b(i);
end

if ( nargout == 0 )
disp ( r )
else
y = r;
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function x = sor ( A, b, xold, omega, TOL, Nmax )

%SOR approximate the solution of the linear system Ax = b by applying

% the SOR method (successive over-relaxation)
%
% calling sequences:
% x = sor ( A, b, xold, omega, TOL, Nmax )
% sor ( A, b, xold, omega, TOL, Nmax )
%
% inputs:
% A coefficient matrix for linear system - must be a
% square matrix
% b right-hand side vector for linear system
% xold vector containing initial guess for solution of
% linear system
% omega relaxation parameter
% TOL convergence tolerance - applied to maximum norm of
% difference between successive approximations
% NMax maximum number of iterations to be performed
%
% output:
% x approximate solution of linear system
%
% NOTE:
% if SOR is called with no output arguments, the
% iteration number and the current approximation to the
% solution are displayed
%
% if the maximum number of iterations is exceeded, a meesage
% to this effect will be displayed and the current approximation
% will be returned in the output value
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

n = length ( b );
[r c] = size ( A );
if ( c ~= n )
disp ( 'sor error: matrix dimensions and vector dimension not compatible' )
return
end
xnew = zeros ( 1, n );

if ( nargout == 0 )

53
 s = sprintf ( '%3d \t %10f ', 0, xold(1) );
for j = 2 : n
s = sprintf ( '%s%10f ', s, xold(j) );
end
disp ( s );
end

for its = 1 : Nmax

xnew(1) = ( 1 - omega ) * xold(1) + omega *...
( b(1) - sum ( A(1,2:n) .* xold(2:n) ) ) / A(1,1);
for i = 2 : n-1
xnew(i) = ( 1 - omega ) * xold(i) + omega * ( b(i) -...
sum ( A(i,1:i-1) .* xnew(1:i-1) ) -...
sum ( A(i,i+1:n) .* xold(i+1:n) ) ) / A(i,i);
end
xnew(n) = ( 1 - omega ) * xold(n) + omega *...
( b(n) - sum ( A(n,1:n-1) .* xnew(1:n-1) ) ) / A(n,n);

if ( nargout == 0 )
s = sprintf ( '%3d \t %10f ', its, xnew(1) );
for j = 2 : n
s = sprintf ( '%s%10f ', s, xnew(j) );
end
disp ( s );
end

conv = max ( abs ( xnew - xold ) );

if ( conv < TOL )
x = xnew;
return
else
xold = xnew;
end
end
disp ( 'sor error: maximum number of iterations exceeded' );
if ( nargout == 1 )
x = xnew;
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [lambda, v] = sympower ( A, x, TOL, Nmax )

%SYMPOWER approximate the dominant eigenvalue and an associated

% eigenvector for a symmetric matrix using the power
% method
%
% calling sequences:
% [lambda, v] = sympower ( A, x, TOL, Nmax )
% lambda = sympower ( A, x, TOL, Nmax )
% sympower ( A, x, TOL, Nmax )
%
% inputs:
% A square symmetric matrix whose dominant eigenvalue
% is to be approximated
% x initial approximation to eigenvector corresponding
% to the dominant eigenvalue
% TOL absolute error convergence tolerance
% (convergence is measured in terms of the Euclidean
% norm of the difference between successive terms
% in the eigenvector seqeunce)

54
% Nmax maximum number of iterations to be performed
%
% outputs:
% lambda approximation to dominant eigenvalue of A
% v an eigenvector of A corresponding to the eigenvalue
% lambda - vector will be normalized to unit length
% in the Euclidean norm
%
% NOTE:
% if SYMPOWER is called with no output arguments, the
% iteration number, the current eigenvector approximation,
% the current eigenvalue approximation and an estimate of
% the rate of convergence of the eigenvalue sequence will
% be displayed
%
% if the maximum number of iterations is exceeded, a message
% to this effect will be displayed and the most recent
% approximations to the dominant eigenvalue and its corresponding
% eigenvector will be returned in the output values
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

[r c] = size ( A );
[rx rc] = size ( x );
if ( rx == 1 ) x = x'; rx = rc; end;

if ( r ~= c )
disp ( 'sympower error: matrix must be square' );
return;
elseif ( r ~= rx )
disp ( 'sympower error: dimensions of matrix and vector are not compatible' );
return;
end

x = x / sqrt ( sum ( x .* x ) );
mu_old = 0;

if ( nargout == 0 )
s = sprintf ( '%3d \t %10f ', 0, x(1) );
for j = 2 : rx
s = sprintf ( '%s%10f ', s, x(j) );
end
disp ( s );
end

for i = 1 : Nmax
xnew = A * x;
mu = sum ( x .* xnew );
xnew = xnew / sqrt ( sum ( xnew .* xnew ) );

if ( nargout == 0 )
s = sprintf ( '%3d \t %10f ', i, xnew(1) );
for j = 2 : rx
s = sprintf ( '%s%10f ', s, xnew(j) );
end
s = sprintf ( '%s \t %10f', s, mu );
if ( i >= 3 )
s = sprintf ( '%s \t \t %10f', s, abs((mu-mu_old)/(mu_old-mu_older)) );

55
 end
disp ( s );
end

delta = x - sign(mu) * xnew;

err = sqrt ( sum ( delta .* delta ) );
if ( err < TOL )
if ( nargout >= 1 ) lambda = mu; end;
if ( nargout >= 2 ) v = xnew; end;
return;
else
x = xnew;
mu_older = mu_old;
mu_old = mu;
end

end

disp ( 'sympower error: Maximum number of iteration exceeded' );

if ( nargout >= 1 )
lambda = mu;
end
if ( nargout >= 2 )
v = xnew;
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [lambda, vects] = hotelling ( A, num_pairs, TOL, Nmax )

%HOTELLING apply Hotelling deflation to a symmetric matrix to

% approximate a specified number of eigenpairs associated
% with the largest eigenvalues of the matrix - the routine
% SYMPOWER is used to approximate each eigenpair
%
% calling sequences:
% [lambda, vects] = hotelling ( A, num_pairs, TOL, Nmax )
% lambda = hotelling ( A, num_pairs, TOL, Nmax )
% hotelling ( A, num_pairs, TOL, Nmax )
%
% inputs:
% A square symmetric matrix whose eigenpairs are to be
% approximated
% num_pairs
% number of eigenpairs to approximate
% TOL absolute error convergence tolerance applied during
% each call to SYMPOWER
% (convergence is measured in terms of the Euclidean
% norm of the difference between successive terms
% in the eigenvector seqeunce)
% Nmax maximum number of iterations to be performed during
% each call to SYMPOWER
%
% outputs:
% lambda vector containing the largest 'num_pairs' eigenvalues
% of the matrix A
% vects matrix containing eigenvectors corresponding to the
% entries in the output vector 'lambda'
% - the i-th column of this matrix is an eigenvector,
% normalized to unit length in the Euclidean norm,
% corresponding to the i-th entry in the vector

56
% 'lambda'
%
% dependencies:
% this routine makes use of the routine SYMPOWER
%
% NOTE:
% if the maximum number of iterations is exceeded during any
% call to SYMPOWER, a message to this effect will be displayed
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

[r c] = size ( A );
if ( r ~= c )
disp ( 'hotelling error: matrix must be square' );
return;
end
n = r;

l = zeros ( 1, num_pairs );
v = zeros ( n, num_pairs );

[l(1) v(:,1)] = sympower ( A, rand(n,1), TOL, Nmax );

for i = 2:num_pairs
A = hd ( A, l(i-1), v(:,i-1) );
[l(i) v(:,i)] = sympower ( A, rand(n,1), TOL, Nmax );
end

lambda = l;
vects = v;

return;
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function B = hd ( A, l, v )
B = A - ( l / dot ( v, v ) ) * v * v';
return;

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [lambda, v] = inv_power ( A, q, x, TOL, Nmax )

%INV_POWER approximate the eigenvalue nearest to the number q,

% and an associated eigenvector, for an arbitrary matrix
% using the inverse power method
%
% calling sequences:
% [lambda, v] = inv_power ( A, q, x, TOL, Nmax )
% lambda = inv_power ( A, q, x, TOL, Nmax )
% inv_power ( A, q, x, TOL, Nmax )
%
% inputs:
% A square matrix whose eigenvalue nearest to the value
% q is to be approximated
% q approximation to an eigenvalue of A
% x initial approximation to eigenvector corresponding
% to the eigenvalue nearest to q
% TOL absolute error convergence tolerance
% (convergence is measured in terms of the maximum
% norm of the difference between successive terms
% in the eigenvector seqeunce)

57
% Nmax maximum number of iterations to be performed
%
% outputs:
% lambda approximation to dominant eigenvalue of A
% v an eigenvector of A corresponding to the eigenvalue
% lambda - vector will be normalized to unit length
% in the maximum norm
%
% dependencies:
% this routine makes use of both LUfactor and LUsolve from
% "Systems of Equations" library
%
% NOTE:
% if INV_POWER is called with no output arguments, the
% iteration number, the current eigenvector approximation,
% the current eigenvalue approximation and an estimate of
% the rate of convergence of the eigenvalue sequence will
% be displayed
%
% if the maximum number of iterations is exceeded, a message
% to this effect will be displayed and the most recent
% approximations to the eigenvalue nearest to q and its
% corresponding eigenvector will be returned in the output values
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

[r c] = size ( A );
[rx rc] = size ( x );
if ( rc == 1 ) x = x'; rc = rx; end;

if ( r ~= c )
disp ( 'inv_power error: matrix must be square' );
return;
elseif ( r ~= rc )
disp ( 'inv_power error: dimensions of matrix and vector are not compatible' );
return;
end;

p = min ( find ( abs(x) == max(abs(x)) ) );

x = x / x(p);
mu_old = 0;

if ( nargout == 0 )
s = sprintf ( '%3d \t %10f ', 0, x(1) );
for j = 2 : rc
s = sprintf ( '%s%10f ', s, x(j) );
end;
disp ( s );
end;

[lu pvt] = LUfactor ( A - q*eye(rc) );

for i = 1 : Nmax
xnew = LUsolve ( lu, x, pvt );
mu = xnew(p);
p = min ( find ( abs(xnew) == max(abs(xnew)) ) );
xnew = xnew / xnew(p);

58
 if ( nargout == 0 )
s = sprintf ( '%3d \t %10f ', i, xnew(1) );
for j = 2 : rc
s = sprintf ( '%s%10f ', s, xnew(j) );
end
s = sprintf ( '%s \t %10f', s, 1/mu+q );
if ( i >= 2 )
s = sprintf ( '%s \t \t %10f', s, abs((mu-mu_old)/(mu_old-mu_older)) );
end
disp ( s );
end

err = max ( abs ( x - xnew ) );

if ( err < TOL )
if ( nargout >= 1 ) lambda = 1/mu+q; end;
if ( nargout >= 2 ) v = xnew'; end;
return;
else
x = xnew;
mu_older = mu_old;
mu_old = mu;
end

end

disp ( 'inv_power error: Maximum number of iteration exceeded' );

if ( nargout >= 1 )
lambda = 1/mu+q;
end
if ( nargout >= 2 )
v = xnew;
end

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [lambda, vects] = wielandt ( A, num_pairs, TOL, Nmax )

%WIELANDT apply Wielandt deflation to an arbitrary matrix to

% approximate a specified number of eigenpairs associated
% with the either largest or the smallest eigenvalues of
% the matrix - the routine POWER_METHOD is used to approximate
% eigenpairs associated with the largest eigenvalues, the
% routine INV_POWER is used to approximate eigenpairs
% associated with the smallest eigenvalues
%
% calling sequences:
% [lambda, vects] = wielandt ( A, num_pairs, TOL, Nmax )
% lambda = wielandt ( A, num_pairs, TOL, Nmax )
% wielandt ( A, num_pairs, TOL, Nmax )
%
% inputs:
% A square symmetric matrix whose eigenpairs are to be
% approximated
% num_pairs
% number of eigenpairs to approximate
% if num_pairs > 0 then approximate largest eigenvalues
% if num_pairs < 0 then approximate smallest eigenvalues
% TOL absolute error convergence tolerance applied during
% each call to either POWER_METHOD or INV_POWER
% (convergence is measured in terms of the maximum
% norm of the difference between successive terms

59
% in the eigenvector seqeunce)
% Nmax maximum number of iterations to be performed during
% each call to either POWER_METHOD or INV_POWER
%
% outputs:
% lambda vector containing the largest/smallest 'num_pairs'
% eigenvalues of the matrix A
% vects matrix containing eigenvectors corresponding to the
% entries in the output vector 'lambda'
% - the i-th column of this matrix is an eigenvector,
% normalized to unit length in the maximum norm,
% corresponding to the i-th entry in the vector
% 'lambda'
%
% dependencies:
% this routine makes use of the routines POWER_METHOD and
% INV_POWER
%
% NOTE:
% if the maximum number of iterations is exceeded during any
% call to either POWER_METHOD or INV_POWER, a message to this
% effect will be displayed
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

[r c] = size ( A );
if ( r ~= c )
disp ( 'wielandt error: matrix must be square' );
return;
end;
n = r;

if ( num_pairs < 0 )
num_pairs = abs(num_pairs);
small = 1;
else
small = 0;
end;

l = zeros ( 1, num_pairs );
v = zeros ( n, num_pairs );
j = zeros ( 1, num_pairs );
x = zeros ( num_pairs, n );

if ( small == 0 )
[l(1) v(:,1)] = power_method ( A, rand(n,1), TOL, Nmax );
else
[l(1) v(:,1)] = inv_power ( A, 0, rand(1,n), TOL, Nmax );
end;

for i = 2:num_pairs
j(i) = min ( find ( abs(v(:,i-1)) == max(abs(v(:,i-1))) ) );
x(i,1:n+2-i) = A(j(i),:);
A = wd ( A, n+2-i, j(i), v(:,i-1) );
if ( small == 0 )
[l(i) v(1:n-i+1,i)] = power_method ( A, rand(n-i+1,1), TOL, Nmax );
else
[l(i) v(1:n-i+1,i)] = inv_power ( A, 0, rand(1,n-i+1), TOL, Nmax );
end;

60
end;

for i=num_pairs:-1:2
for k=i:-1:2
temp = [v(1:j(k)-1, i); 0; v(j(k):n-k+1, i)];
v(1:n-k+2,i) = (l(i) - l(k-1))*temp +(x(k,1:n+2-k)*temp)*v(1:n-k+2,k-
1)/v(j(k),k-1);
end;
p = min ( find ( abs(v(:,i)) == max(abs(v(:,i))) ) );
v(:,i) = v(:,i) / v(p,i);
end;

lambda = l;
vects = v;

return;
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function B = wd ( A, n, j, v )

temp = ( v / v(j) ) * A(j,:);

B = [A(1:j-1,1:j-1)-temp(1:j-1,1:j-1) A(1:j-1,j+1:n) - temp(1:j-1,j+1:n); ...
A(j+1:n,1:j-1)-temp(j+1:n,1:j-1) A(j+1:n,j+1:n) - temp(j+1:n,j+1:n)];
return;

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [lambda, v] = power_method ( A, x, TOL, Nmax )

%POWER_METHOD approximate the dominant eigenvalue and an associated

% eigenvector for an arbitrary matrix using the power
% method
%
% calling sequences:
% [lambda, v] = power_method ( A, x, TOL, Nmax )
% lambda = power_method ( A, x, TOL, Nmax )
% power_method ( A, x, TOL, Nmax )
%
% inputs:
% A square matrix whose dominant eigenvalue is to be
% approximated
% x initial approximation to eigenvector corresponding
% to the dominant eigenvalue
% TOL absolute error convergence tolerance
% (convergence is measured in terms of the maximum
% norm of the difference between successive terms
% in the eigenvector seqeunce)
% Nmax maximum number of iterations to be performed
%
% outputs:
% lambda approximation to dominant eigenvalue of A
% v an eigenvector of A corresponding to the eigenvalue
% lambda - vector will be normalized to unit length
% in the maximum norm
%
% NOTE:
% if POWER_METHOD is called with no output arguments, the
% iteration number, the current eigenvector approximation,
% the current eigenvalue approximation and an estimate of
% the rate of convergence of the eigenvalue sequence will
% be displayed
%

61
% if the maximum number of iterations is exceeded, a message
% to this effect will be displayed and the most recent
% approximations to the dominant eigenvalue and its corresponding
% eigenvector will be returned in the output values
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

[r c] = size ( A );
[rx rc] = size ( x );
if ( rx == 1 ) x = x'; rx = rc; end;

if ( r ~= c )
disp ( 'power_method error: matrix must be square' );
return;
elseif ( r ~= rx )
disp ( 'power_method error: dimensions of matrix and vector are not compatible'
);
return;
end;

p = min ( find ( abs(x) == max(abs(x)) ) );

x = x / x(p);
mu_old = 0;

if ( nargout == 0 )
s = sprintf ( '%3d \t %10f ', 0, x(1) );
for j = 2 : rx
s = sprintf ( '%s%10f ', s, x(j) );
end;
disp ( s );
end;

for i = 1 : Nmax
xnew = A * x;
mu = xnew(p);
p = min ( find ( abs(xnew) == max(abs(xnew)) ) );
xnew = xnew / xnew(p);

if ( nargout == 0 )
s = sprintf ( '%3d \t %10f ', i, xnew(1) );
for j = 2 : rx
s = sprintf ( '%s%10f ', s, xnew(j) );
end;
s = sprintf ( '%s \t %10f', s, mu );
if ( i >= 2 )
s = sprintf ( '%s \t \t %10f', s, abs((mu-mu_old)/(mu_old-mu_older)) );
end;
disp ( s );
end;

err = max ( abs ( x - xnew ) );

if ( err < TOL )
if ( nargout >= 1 ) lambda = mu; end;
if ( nargout >= 2 ) v = xnew; end;
return;
else
x = xnew;
mu_older = mu_old;
mu_old = mu;

62
 end;

end;
disp ( 'power_method error: Maximum number of iteration exceeded' );
if ( nargout >= 1 ) lambda = mu; end;
if ( nargout >= 2 ) v = xnew;
end.

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [lambda, v] = qrst ( a, b, TOL, Nmax, vects )

%QRST determine all of the eigenvalues (and optionally all of the

% eigenvectors) of a symmetric tridiagonal matrix using the
% QR algorithm with Wilkinson shift
%
% calling sequences:
% [lambda, v] = qrst ( a, b, TOL, Nmax, vects )
% [lambda, v] = qrst ( a, b, TOL, Nmax )
% lambda = qrst ( a, b, TOL, Nmax )
% qrst ( a, b, TOL, Nmax, vects )
% qrst ( a, b, TOL, Nmax )
%
% inputs:
% a vector containing elements along the main diagonal
% of the symmetric tridiagonal matrix whose eigenvalues
% are to be determined
% b vector containing elements along the off diagonal
% of the symmetric tridiagonal matrix whose eigenvalues
% are to be determined
% TOL convergence tolerance
% Nmax maximum number of iterations
% vects optional input argument
% matrix containing eigenvector information produced
% during the reduction of the original symmetric
% matrix to symmetric tridiagonal form
% - this input is needed only if computation of the
% eigenvectors is requested (by including the second
% output argument) and the original matrix was not in
% symmetric tridiagonal form
%
% output:
% lambda vector containing the eigenvalues of the symmetric
% tridiagonal matrix determined by the vectors a and b
% v optional output argument
% matrix containing the eigenvectors of the symmetric
% tridiagonal matrix determined by the vectors a and b
% - the i-th column of this matrix is an eigenvector
% which corresponds to the i-th eigenvalue in the
% vector lambda
% - eigenvectors will be not computed if this second
% output argument is omitted
%
% NOTE:
% if the maximum number of iterations is exceeded, a message
% to this effect will be displayed, along with the number of
% eigenvalues which had been determined - these eigenvalues
% will be returned in the last entries of the output vector
% lambda
%
% Copyright by Reza Abazari 2008-04-23

63
% Email : abazari-r@[Link]

n = length(a);
if ( length(b) == n-1 ) b(2:n) = b(1:n-1); end;

c = zeros ( 1, n );
s = zeros ( 1, n );
shift = 0;
togo = n;

if ( nargout >= 2 )
if ( nargin >= 5 )
v = vects;
else
v = eye(n);
end;
end;

for its = 1 : Nmax

if ( togo == 1 )
lambda(1) = a(1) + shift;
disp ( its );
return;
end;

trace = a(togo-1) + a(togo);

det = a(togo-1)*a(togo) - b(togo)*b(togo);
disc = sqrt ( trace*trace - 4*det );
mu1 = (1/2) * ( trace + disc );
mu2 = (1/2) * ( trace - disc );
if ( abs ( mu1 - a(togo) ) < abs ( mu2 - a(togo) ) )
s = mu1;
else
s = mu2;
end;

shift = shift + s;
for i = 1:togo
a(i) = a(i) - s;
end;

oldb = b(2);
for i = 2:togo
j = i-1;
r = sqrt ( a(j)^2 + oldb^2 );
c(i) = a(j) / r;
s(i) = oldb / r;
a(j) = r;
temp1 = c(i)*b(i) + s(i)*a(i);
temp2 = -s(i)*b(i) + c(i)*a(i);
b(i) = temp1;
a(i) = temp2;
if ( i ~= togo ) oldb = b(i+1); b(i+1) = c(i)*b(i+1); end;
end;

a(1) = c(2)*a(1) + s(2)*b(2);

b(2) = s(2)*a(2);
for i = 2:togo-1
a(i) = s(i+1)*b(i+1) + c(i)*c(i+1)*a(i);

64
 b(i+1) = s(i+1)*a(i+1);
end;
a(togo) = c(togo)*a(togo);

if ( nargout >= 2 )
for i = 2 : togo
col1 = v(:,i-1) * c(i) + v(:,i) * s(i);
v(:,i) = -s(i) * v(:,i-1) + c(i) * v(:,i);
v(:,i-1) = col1;
end;
end;

if ( abs(b(togo)) < TOL )

lambda(togo) = a(togo) + shift;
disp([lambda(togo) its]);
togo = togo - 1;
end;

end;

disp ( 'qrst error: Maximum number of iterations exceeded' );

disp ( sprintf ( '%d eigenvalues determined \n', n-togo ) );

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [a, b, V] = red2st ( A )

%RED2ST perform similarity transformations to reduce the symmetric

% matrix A to symmetric tridiagonal form
%
% calling sequences:
% [a, b, V] = red2st ( A )
% [a, b] = red2st ( A )
% red2st ( A )
%
% input:
% A square symmetric matrix to be reduced to symmetric
% tridiagonal form
%
% outputs:
% a vector containing elements along the main diagonal
% of the symmetric tridiagonal from of A
% b vector containing elements along the off diagonal
% of the symmetric tridiagonal from of A
% V optional output argument
% matrix containing eigenvector information for the
% matrix A
%
% Copyright by Reza Abazari 2008-04-23
% Email : abazari-r@[Link]

[nrow ncol] = size ( A );

if ( nrow ~= ncol )
disp ( 'red2st error: square matrix required' );
return;
end;
n = nrow;

if ( nargout >= 3 ) V = eye(n); end;

65
for i = 1 : n-2
w = zeros ( n, 1 );
x = A(:,n-i+1);
alpha = - sign(x(n-i)) * norm ( x(1:n-i) );
if ( alpha ~= 0 )
w(n-i) = sqrt ( (1/2) * ( 1 - x(n-i)/alpha ) );
w(1:n-i-1) = -(1/2) * x(1:n-i-1) / ( alpha * w(n-i) );

u = A * w;
K = dot ( w, u );
q = u - K * w;
A = A - 2*w*q' - 2*q*w';

if ( nargout >= 3 ) V = V - 2*V*w*w'; V; end;

end;
end;

a = diag(A);
b = zeros ( n, 1 );
for i = 2:n b(i) = A(i,i-1);
end;

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

3. Polynomial Interpolation

66