Numerical Mathematics With Matlab

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 : [email protected] 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 : [email protected] 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 : [email protected] 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 : [email protected] 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 : [email protected] 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
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
7

% Copyright by Reza Abazari 2008-04-23 % Email : [email protected] 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 (y0big | 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
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 : [email protected] % 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 % % % % %Output % %

-

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 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 : [email protected] %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))
- 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 - 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 : [email protected] ya=f(a); yb=f(b); if ya*yb rel="nofollow">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 % % % % %Output % % %

-

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 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 : [email protected] 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 : [email protected] 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 % % % % %Output % % % % % %Call

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

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)

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

% Copyright by Reza Abazari 2008-04-23 % Email : [email protected]

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
20

yC2=F(C2); if (yC2=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 technique.this 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 : [email protected] 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 : [email protected]

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 : [email protected] 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 M.file used findzero.m % findzero.m

and err.m

that’s given the end of this m.file.

% Copyright by Reza Abazari 2008-12-09 % Email : [email protected] 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];

N= feval(fun,y);

A=M\N;

% the matrix of the LHS of the % linear system of equations % the column vector of the RHS % of the linear system of equations

% % % %

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
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 : [email protected] % 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 m.file: %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); polynomial E=A(end);

% the 3 coefficients of the second order % 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 : [email protected] % 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 : [email protected] % % % % % % % % % % % % % %

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); % f = feval(fun,x); %

Vector of potential bracket limits 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 : [email protected] 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 : [email protected]

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 : [email protected]

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 % % % calling % % % % inputs: % % % % % % % output: % % % NOTE: % % % % %

locate the fixed point of an arbitrary function using general fixed_point iteration sequences: y = fixed_point ( 'g', x0, TOL, Nmax ) fixed_point ( 'g', x0, TOL, Nmax )

g x0 TOL Nmax

string containing name of m-file defining the iteration function initial approximation to location of fixed point absolute error convergence tolerance maximum number of iterations to be performed

y

approximate value of fixed point

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 : [email protected]

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 : [email protected]

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 % % % calling % % % % inputs: % % % % % % %

locate the roots of an arbitrary polynomial using Laguerre's method sequences: r = polyroots ( poly, Nmax, TOL ) polyroots ( poly, Nmax, TOL )

poly

Nmax

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. 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 : [email protected]

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 : [email protected] 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 : [email protected]

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 : [email protected]

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 : [email protected] 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 : [email protected] 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 % % % % % %Output %

- 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 - 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 : [email protected] 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 : [email protected] [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 % % % calling % % % % inputs: % % % % % % % % % % output: % % % NOTE: % % % % % % % %

approximate the solution of the linear system Ax = b by applying the Gauss-Seidel method (successive relaxation) sequences: x = gauss_seidel ( A, b, xold, TOL, Nmax ) gauss_seidel ( A, b, xold, TOL, Nmax )

A

NMax

coefficient matrix for linear system - must be a square matrix right-hand side vector for linear system vector containing initial guess for solution of linear system convergence tolerance - applied to maximum norm of difference between successive approximations maximum number of iterations to be performed

x

approximate solution of linear system

b xold TOL

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 : [email protected] 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 % % % calling % % % % inputs:

solve the system of nonlinear equations F(x) = 0 using Broyden's method sequences: y = broyden ( F, J, x0, TOL, Nmax ) broyden ( F, J, x0, TOL, Nmax )

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 : [email protected]

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 : [email protected]

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 : [email protected] 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 : [email protected]

[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 : [email protected]

[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 : [email protected]

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 : [email protected]

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 : [email protected]

n = length ( b [r c] = size ( if ( c ~= n ) disp ( 'sor return end xnew = zeros (

); A ); error: matrix dimensions and vector dimension not compatible' )

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 % % % % calling % % % % % inputs: % % % % % % % %

approximate the dominant eigenvalue and an associated eigenvector for a symmetric matrix using the power method sequences: [lambda, v] = sympower ( A, x, TOL, Nmax ) lambda = sympower ( A, x, TOL, Nmax ) sympower ( A, x, TOL, Nmax )

A x TOL

square symmetric matrix whose dominant eigenvalue is to be approximated initial approximation to eigenvector corresponding to the dominant eigenvalue 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 : [email protected]

[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 : [email protected]

[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 % % % % calling % % % % % inputs: % % % % % % % % %

approximate the eigenvalue nearest to the number q, and an associated eigenvector, for an arbitrary matrix using the inverse power method 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 )

A q x TOL

square matrix whose eigenvalue nearest to the value q is to be approximated approximation to an eigenvalue of A initial approximation to eigenvector corresponding to the eigenvalue nearest to q 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 : [email protected]

[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) x = x / x(p); mu_old = 0;

== max(abs(x)) ) );

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 % % % % % % % % calling % % % % % inputs: % % % % % % % % % %

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 sequences: [lambda, vects] = wielandt ( A, num_pairs, TOL, Nmax ) lambda = wielandt ( A, num_pairs, TOL, Nmax ) wielandt ( A, num_pairs, TOL, Nmax )

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 : [email protected]

[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 v j x

= = = =

zeros zeros zeros zeros

( ( ( (

1, num_pairs n, num_pairs 1, num_pairs 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,k1)/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 : [email protected]

[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) x = x / x(p); mu_old = 0;

== max(abs(x)) ) );

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) xnew = xnew / xnew(p);

== max(abs(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 >= 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 : [email protected]

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 : [email protected]

[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 = K q A

A = = =

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

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

V; end;

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

3. Polynomial Interpolation 66