Chapter 1 M-Files

M-File Home

Script File Index

Function File Index

SineTable

% Script File: SineTable
% Prints a short table of sine evaluations.

clc
n = 21;
x = linspace(0,1,n);
y = sin(2*pi*x);
disp(' ')
disp('  k     x(k)   sin(x(k))')
disp('------------------------')
for k=1:21
   degrees = (k-1)*360/(n-1);
   disp(sprintf(' %2.0f     %3.0f     %6.3f    ',k,degrees,y(k)));
end
disp( ' ');
disp('x(k) is given in degrees.')
disp(sprintf('One Degree = %5.3e Radians',pi/180))

SinePlot

% Script File: SinePlot
% Displays increasingly smooth plots of sin(2*pi*x).

close all
for n = [4 8 12 16 20 50 100 200 400]
   x = linspace(0,1,n);
   y = sin(2*pi*x);
   plot(x,y)
   title(sprintf('Plot of sin(2*pi*x) based upon n = %3.0f points.',n))
   pause(1)
end

ExpPlot

% Script File: ExpPlot
% Examines the function
%
%      f(x) = ((1 + x/24)/(1 - x/24 + x^2/384))^8
%
% as an approximation to exp(z) across [0,1].

close all
x  = linspace(0,1,200);
num   = 1 + x/24;
denom = 1 - x/12 + (x/384).*x;
quot = num./denom;
y = quot.^8;
plot(x,y,x,exp(x))

TangentPlot

% Script File: TangentPlot
% Plots the function tan(x), -pi/2 <= x <= 9pi/2

close all
ymax = 10;
x = linspace(-pi/2,pi/2,40);
y = tan(x);
plot(x,y)
axis([-pi/2 9*pi/2 -ymax ymax])
title('The Tangent Function')
xlabel('x')
ylabel('tan(x)')
hold on
for k=1:4
   xnew = x+ k*pi;
   plot(xnew,y);
end
hold off

SineAndCosPlot

% Script File: SineAndCosPlot
% Plots the functions sin(2*pi*x) and cos(2*pi*x) across [0,1]
% and marks their intersection.

close all
x = linspace(0,1,200);
y1 = sin(2*pi*x);
y2 = cos(2*pi*x);
plot(x,y1,x,y2,'--',[1/8 5/8],[1/sqrt(2) -1/sqrt(2)],'*')

Polygons

% Script File: Polygons
% Plots selected regular polygons.

close all
theta = linspace(0,2*pi,361);
c = cos(theta);
s = sin(theta);
k=0;
for sides = [3 4 5 6 8 10 12 18 24]
   stride = 360/sides;
   k=k+1;
   subplot(3,3,k)
   plot(c(1:stride:361),s(1:stride:361))
   axis([-1.2 1.2 -1.2 1.2])
   axis equal
end

SumOfSines

% Script File: SumOfSines
% Plots f(x) = 2sin(x) + 3sin(2x) + 7sin(3x) + 5sin(4x)
% across the interval [-10,10].

close all
x = linspace(-10,10,200)';
A = [sin(x) sin(2*x) sin(3*x) sin(4*x)];
y = A*[2;3;7;5];
plot(x,y)
title('f(x) = 2sin(x) + 3sin(2x) + 7sin(3x)  + 5sin(4x)')

SumOfSines2

% Script File: SumOfSines2
% Plots the functions
%         f(x) = 2sin(x) + 3sin(2x) + 7sin(3x) + 5sin(4x)
%         g(x) = 8sin(x) + 2sin(2x) + 6sin(3x) + 9sin(4x)
% across the interval [-10,10].
   
close all
n = 200;
x = linspace(-10,10,n)';
A = [sin(x) sin(2*x) sin(3*x) sin(4*x)];
y = A*[2 8;3 2;7 6;5 9];
plot(x,y)

UpDown

% Script File: UpDown
% Generates a column vector x(1:n) of positive integers 
% where x(1) is solicited and
% 
%          x(k+1) = x(k)/2    if x(k) is even.
%          x(k+1) = 3x(k)+1   if x(k) is odd.
%
% The value of n is either 500 or the first index with the
% property that x(n) = 1, whichever comes first.

x = zeros(500,1);
x(1) = input('Enter initial positive integer:');
k = 1;
while ((x(k) ~= 1) & (k < 500))
   if rem(x(k),2) == 0
      x(k+1) = x(k)/2;
   else
      x(k+1) = 3*x(k)+1;
   end
   k = k+1;
end
n = k;
x = x(1:n);

clc
disp(sprintf('x(1:%1.0f) = \n',n))
disp(sprintf('%-8.0f',x))
[xmax,imax] = max(x);
disp(sprintf('\n x(%1.0f) = %1.0f is the max.',imax,xmax))
density = sum(x<=x(1))/x(1);
disp(sprintf(' The density is %5.3f.',density))

close all
figure
plot(x)
title(sprintf('x(1) = %1.0f, n = %1.0f',x(1),n));
figure
plot(-sort(-x))
title('Sequence values sorted.')
I = find(rem(x(1:n-1),2)==1);
if length(I)>1
   figure
   plot((1:n),zeros(1,n),I+1,x(I+1),I+1,x(I+1),'*')
   title('Local Maxima')
end

RunUpDown

% Script File: RunUpDown
% Environment for studying the up/down sequence.
% Stores selected results in file UpDownOutput.

while(input('Another Example? (1=yes, 0=no)'))
   diary UpDownOutput
   UpDown
   diary off
   disp(' ')
   if (input('Keep Output? (1=yes, 0=no)')~=1)
      delete UpDownOutput
   end
end

Histograms

% Script File: Histograms
% Histograms of rand(1000,1) and randn(1000,1).

close all
subplot(2,1,1)
x = rand(1000,1);
hist(x,30)
axis([-1 2 0 60])
title('Distribution of Values in rand(1000,1)')
xlabel(sprintf('Mean = %5.3f. Median = %5.3f.',mean(x),median(x)))
subplot(2,1,2)
x = randn(1000,1);
hist(x,linspace(-2.9,2.9,100));
title('Distribution of Values in randn(1000,1)')
xlabel(sprintf('Mean = %5.3f. Standard Deviation = %5.3f',mean(x),std(x)))

Clouds

% Script File: Clouds
% 2-dimensional pictures of the uniform and normal distributions.

close all
Points = rand(1000,2);
subplot(1,2,1)
plot(Points(:,1),Points(:,2),'.')     
title('Uniform Distribution.')
axis([0 1 0 1])
axis square
Points = randn(1000,2);
subplot(1,2,2)
plot(Points(:,1),Points(:,2),'.')    
title('Normal Distribution.')
axis([-3 3 -3 3])
axis square

Dice

% Script File: Dice
% Simulates 1000 rollings of a pair of dice.

close all
First  = 1 + floor(6*rand(1000,1));
Second = 1 + floor(6*rand(1000,1));
Throws = First + Second;
hist(Throws, linspace(2,12,11));
title('Outcome of 1000 Dice Rolls.')

Darts

% Script File: Darts
% Estimates pi using random dart throws.

close all
rand('seed',.123456) 
NumberInside = 0;
PiEstimate = zeros(500,1);
for k=1:500
   x = -1+2*rand(100,1);
   y = -1+2*rand(100,1);
   NumberInside = NumberInside + sum(x.^2 + y.^2 <= 1);
   PiEstimate(k) = (NumberInside/(k*100))*4;    
end
plot(PiEstimate)
title(sprintf('Monte Carlo Estimate of Pi = %5.3f',PiEstimate(500)));
xlabel('Hundreds of Trials')

Smooth

% Script File: Smooth
% Solicits n, draws an n-gon, and then smooths it.

close all
n = input('Enter the number of edges:');
figure
axis([0 1 0 1])
axis square
hold on
x = zeros(n,1);
y = zeros(n,1);
for k=1:n
   title(sprintf('Click in %2.0f more points.',n-k+1))
   [x(k) y(k)] = ginput(1);
   plot(x(1:k),y(1:k), x(1:k),y(1:k),'*')
end
x0 = x; y0 = y;
x = [x;x(1)];
y = [y;y(1)];
plot(x,y,x,y,'*')
title('The Original Polygon') 
hold off
k=0;
xlabel('Click inside window to smooth, outside window to quit.')  
[a,b] = ginput(1);
v = axis;
while (v(1)<=a) & (a<=v(2)) & (v(3)<=b) & (b<=v(4));
   k = k+1;
   x = [(x(1:n)+x(2:n+1))/2;(x(1)+x(2))/2]; 
   y = [(y(1:n)+y(2:n+1))/2;(y(1)+y(2))/2]; 
   m = max(abs([x;y]));
   x = x/m;
   y = y/m;
   figure
   plot(x,y,x,y,'*')

   axis square
   title(sprintf('Number of Smoothings = %1.0f',k)) 
   xlabel('Click inside window to smooth, outside window to quit.')   
   v = axis;
   [a,b] = ginput(1);
end

Stirling

% Script File: Stirling
% Prints a table showing error in Stirling's formula for n!

close all
clc
disp( '  ')
disp('                       Stirling          Absolute    Relative')
disp('   n          n!     Approximation         Error       Error')
disp('----------------------------------------------------------------')
e=exp(1);
nfact=1;
for n=1:13
   nfact = n*nfact;
   s = sqrt(2*pi*n)*((n/e)^n);
   abserror = abs(nfact - s);
   relerror = abserror/nfact;
   s1 = sprintf('  %2.0f   %10.0f   %13.2f',n,nfact,s);
   s2 = sprintf('   %13.2f    %5.2e',abserror,relerror);
   disp([s1 s2])
end

ExpTaylor

  % Script File: ExpTaylor
  % Plots, as a function of n, the relative error in the 
  % Taylor approximation
  %
  %           1 + x + x^2/2! +...+ x^n/n! 
  %
  % to exp(x).

  close all
  nTerms = 50;
  for x=[10 5 1 -1 -5 -10]
     figure
     f = exp(x)*ones(nTerms,1);
     s = 1;
     term = 1;
     for k=1:50
        term = x.*term/k;
        s = s+ term;
        err(k) = abs(f(k) - s);
     end
     relerr = err/exp(x);
     semilogy(1:nTerms,relerr)
     ylabel('Relative Error in Partial Sum.')
     xlabel('Order of Partial Sum.')
     title(sprintf('x = %5.2f',x))	  
  end

Zoom

% Script File: Zoom
% Plots (x-1)^6 near x=1 with increasingly refined scale.
% Evaluation via x^6 - 6x^5 + 15x^4 - 20x^3 + 15x^2 - 6x +1
% leads to severe cancellation.

close all
k=0;
n=100;
for delta = [.1 .01 .008 .007 .005  .003 ]
   x = linspace(1-delta,1+delta,n)';
   y = x.^6 - 6*x.^5 + 15*x.^4 - 20*x.^3 + 15*x.^2  - 6*x + ones(n,1);
   k = k+1;
   subplot(2,3,k)
   plot(x,y,x,zeros(1,n))
   axis([1-delta 1+delta -max(abs(y)) max(abs(y))])
end

FpFacts

% Script File: FpFacts
% Prints some facts about the underlying floating point system.

close all
clc

% p = smallest positive integer so 1+1/2^p = 1.
x=1; p=0; y=1; z=x+y;
while x~=z  
   y=y/2; 
   p=p+1; 
   z=x+y;
end
disp(sprintf('p = %2.0f  is the smallest positive integer so 1+1/2^p = 1.',p))

% q = smallest positive integer so 1/2^q = 0.
x=1; q=0;
while x>0
   x=x/2; 
   q=q+1;
end;
disp(sprintf('q = %2.0f is the smallest positive integer so 1/2^q = 0.',q))


% r = smallest positive integer so 2^r = inf.
x=1; r=0;
while x~=inf  
   x=2*x; 
   r=r+1; 
end
disp(sprintf('r = %2.0f is the smallest positive integer so 2^r = inf.',r))
   
% The Kahan example. In exact arithmetic f/e = 0/0:

disp(' ')
echo on
h = 1./2.;
x = 2./3. - h;
y = 3./5. - h;
e = (x+x+x) - h;
f = (y+y+y+y+y)-h;
z = f/e;
echo off
disp(sprintf('\nz = %10.5f',z))

TestMyExp

% Script File: TestMyExp
% Checks MyExpF, MyExp1, MyExp2, MyExp3, and MyExp4.

close all
m = 50;
x = linspace(-1,1,m);
exact = exp(x);
figure
k=0;
y = zeros(1,m);
for n = [4 8 16 20]
   for i=1:m
      y(i) = MyExpF(x(i),n);
   end
   RelErr = abs(exact - y)./exact;
   k=k+1;
   subplot(2,2,k)
   plot(x,RelErr)
   title(sprintf('n = %2.0f',n))
end
  
clc
nRepeat = 1; 

% Note: The following fragment may take a long time to execute
% with this value of nRepeat.

disp('T1 = MyExp1 benchmark.')
disp('T2 = MyExp2 benchmark.')
disp(' ')
disp(' Length(x)     T2/T1')
disp('-----------------------')
for L = 1000:100:2000
   xL = linspace(-1,1,L);
   tic
   for k=1:nRepeat
      y = MyExp1(xL);
   end
   T1 = toc;
   tic
   for k=1:nRepeat
      y = MyExp2(xL);
   end
   T2 = toc;	 	 
   disp(sprintf('%6.0f  %13.6f  ',L,T2/T1))
end
  
disp(' ')
disp(' ')
disp(' Length(x)     MyExp3(x)      MyExp4(x) ')
disp('                 Flops          Flops')
disp('---------------------------------------')
for L = 50:50:300
   xL = linspace(-1,1,L);
   flops(0);
   y = MyExp3(xL);
   f3 = flops;
   flops(0);
   y = MyExp4(xL);
   f4 = flops;	 	 
   disp(sprintf('%6.0f  %13.0f  %13.0f  ',L,f3,f4))
end

TestDerivative

% Script File: TestDerivative
% Numerical differentiation of f(x) = sin(10x)

close all
a = 10;
M2=100; 

% Look at errors across [0,1]. 
% err1 = error with no derivative information.
% err2 = error with 2nd derivative information.
m = 100;
x=linspace(0,1,m);
for i=1:m
   exactDer = a*cos(a*x(i));
   err1(i) = abs(exactDer - Derivative('sin10',x(i)));
   err2(i) = abs(exactDer - Derivative('sin10',x(i),eps,M2));     
end

% Plot err1(1:m)
figure
plot(x,err1);
title('Derivative(''sin10'',a)')
xlabel('a')
ylabel('Der Error')
    
% Plot err2(1:m)
figure
plot(x,err2);
title('Derivative(''sin10'',a,eps,100)')
xlabel('a')
ylabel('Der Error')

Euler

% Script File: Euler
% Sums the series 1 + 1/2 + 1/3 + .. in 3-digit floating point arithmetic.
% Terminates when the addition of the next term does not change
% the value of the running sum.

oldsum = represent(0);
one = represent(1);
sum = one;
k = 1;
while convert(sum) ~= convert(oldsum)
   k = k+1;
   kay  = represent(k);
   term = float(one,kay,'/');
   oldsum = sum;
   sum  = float(sum,term,'+');
end   
clc
disp(['The sum for ' num2str(k) ' or more terms is ' pretty(sum)])

ShowPadeArray

% Script File: ShowPadeArray
% Illustrates the function PadeArray that creates a cell array.

  clc
  P = PadeArray(5,5);
  for i=1:5
     for j=1:5
        disp(sprintf('\n (%1d,%1d) Pade Coefficients:\n',i-1,j-1))
        disp(['    ' sprintf('%7.4f  ',P{i,j}.num)])
        disp(['    ' sprintf('%7.4f  ',P{i,j}.den)])
     end
  end

ShowFonts

% Script File: ShowFonts

close all
HA = 'HorizontalAlign';

fonts = {'Times-Roman' 'Helvetica' 'AvantGarde' 'Bookman' 'Palatino'...
         'ZapfChancery' 'Courier' 'NewCenturySchlbk' 'Helvetica-Narrow'};
for k=1:9
   figure
   plot([-10 100 100 -10 -10],[0  0  60  60  0])
   axis([-10 100 0 60])
   axis off
   v=38;
   F = fonts{k};  
   text(45,55,F,'FontName',F,'FontSize',24,HA,'center')
   text(10,47,'Plain','FontName',F,'FontSize',22,HA,'center')
   text(45,47,'Bold','FontName',F,'Fontweight','bold','FontSize',22,HA,'center')
   text(82,47,'Oblique','FontName',F,'FontAngle','oblique','FontSize',22,HA,'center')
   for size=[22 18 14 12 11 10 9]
      text(10,v,'Matlab','FontName',F,'FontSize',size,HA,'center')
      text(45,v,'Matlab','FontName',F,'FontSize',size,HA,'center','FontWeight','bold')
      text(82,v,'Matlab','FontName',F,'FontSize',size,HA,'center','FontAngle','oblique')
      v = v-6;
   end
end

ShowSymbols

% Script File: ShowSymbols

close all
plot([0 12 12 0 0],[0 0 12 12 0])
axis off
text(6,10.5,'Math Symbols','FontSize',18,'HorizontalAlign','center')
x = 1; x1 = x+.7;
y = 4.6; y1 = y+.7;
z = 9; z1 = z+.7;
n = 12;

text(y,9,'\leftarrow','Fontsize',n);   
text(y,8,'\rightarrow','Fontsize',n);
text(y,7,'\uparrow','Fontsize',n)
text(y,6,'\downarrow','Fontsize',n)
text(y,5,'\Leftarrow','FontSize',n)
text(y,4,'\Rightarrow','FontSize',n)
text(y,3,'\Leftrightarrow','FontSize',n)
text(y,2,'\partial','FontSize',n)


text(x,9,'\neq','FontSize',n)
text(x,8,'\geq','FontSize',n)
text(x,7,'\approx','FontSize',n)
text(x,6,'\equiv','Fontsize',n)
text(x,5,'\cong','Fontsize',n)
text(x,4,'\pm','Fontsize',n)
text(x,3,'\nabla','FontSize',n)
text(x,2,'\angle','FontSize',n)


text(z,9,'\in','FontSize',n)
text(z,8,'\subset','Fontsize',n)
text(z,7,'\cup','Fontsize',n)
text(z,6,'\cap','Fontsize',n)
text(z,5,'\perp','FontSize',n)
text(z,4,'\infty','FontSize',n)
text(z,3,'\int','Fontsize',n)
text(z,2,'\times','FontSize',n)

text(y1,9,'\\leftarrow','Fontsize',n);   
text(y1,8,'\\rightarrow','Fontsize',n);
text(y1,7,'\\uparrow','Fontsize',n)
text(y1,6,'\\downarrow','Fontsize',n)
text(y1,5,'\\Leftarrow','FontSize',n)
text(y1,4,'\\Rightarrow','FontSize',n)
text(y1,3,'\\Leftrightarrow','FontSize',n)
text(y1,2,'\\partial','FontSize',n)


text(x1,9,'\\neq','FontSize',n)
text(x1,8,'\\geq','FontSize',n)
text(x1,7,'\\approx','FontSize',n)
text(x1,6,'\\equiv','Fontsize',n)
text(x1,5,'\\cong','Fontsize',n)
text(x1,4,'\\pm','Fontsize',n)
text(x1,3,'\\nabla','FontSize',n)
text(x1,2,'\\angle','FontSize',n)


text(z1,9,'\\in','FontSize',n)
text(z1,8,'\\subset','Fontsize',n)
text(z1,7,'\\cup','Fontsize',n)
text(z1,6,'\\cap','Fontsize',n)
text(z1,5,'\\perp','FontSize',n)
text(z1,4,'\\infty','FontSize',n)
text(z1,3,'\\int','Fontsize',n)
text(z1,2,'\\times','FontSize',n)

ShowGreek

% Script File: ShowGreek

close all
plot([-1 12 12 -1 -1],[-1 -1 12 12 -1])
axis off
text(4,10,'Greek Symbols','FontSize',18)
x = 0; x1 = x+.7;
y = 4; y1 = y+.7;
z = 8; z1 = z+.7;
n = 12;

text(x,8,'\alpha','Fontsize',n);   
text(x,7,'\beta','Fontsize',n);
text(x,6,'\gamma','Fontsize',n)
text(x,5,'\delta','Fontsize',n)
text(x,4,'\epsilon','FontSize',n)
text(x,3,'\kappa','FontSize',n)
text(x,2,'\lambda','Fontsize',n)
text(x,1,'\mu','Fontsize',n)
text(x,0,'\nu','Fontsize',n)

text(x1,8,'\\alpha','Fontsize',n);   
text(x1,7,'\\beta','Fontsize',n);
text(x1,6,'\\gamma','Fontsize',n)
text(x1,5,'\\delta','Fontsize',n)
text(x1,4,'\\epsilon','FontSize',n)
text(x1,3,'\\kappa','FontSize',n)
text(x1,2,'\\lambda','Fontsize',n)
text(x1,1,'\\mu','Fontsize',n)
text(x1,0,'\\nu','Fontsize',n)


text(y,8,'\omega','Fontsize',n)
text(y,7,'\phi','Fontsize',n)
text(y,6,'\pi','FontSize',n)
text(y,5,'\chi','Fontsize',n)
text(y,4,'\psi','Fontsize',n)
text(y,3,'\rho','FontSize',n)
text(y,2,'\sigma','Fontsize',n)
text(y,1,'\tau','FontSize',n)
text(y,0,'\upsilon','FontSize',n)

text(y1,8,'\\omega','Fontsize',n)
text(y1,7,'\\phi','Fontsize',n)
text(y1,6,'\\pi','FontSize',n)
text(y1,5,'\\chi','Fontsize',n)
text(y1,4,'\\psi','Fontsize',n)
text(y1,3,'\\rho','FontSize',n)
text(y1,2,'\\sigma','Fontsize',n)
text(y1,1,'\\tau','FontSize',n)
text(y1,0,'\\upsilon','FontSize',n)

text(z,8,'\Sigma','FontSize',n)
text(z,7,'\Pi','FontSize',n)
text(z,6,'\Lambda','FontSize',n)
text(z,5,'\Omega','FontSize',n)
text(z,4,'\Gamma','FontSize',n)




text(z1,8,'\\Sigma','FontSize',n)
text(z1,7,'\\Pi','FontSize',n)
text(z1,6,'\\Lambda','FontSize',n)
text(z1,5,'\\Omega','FontSize',n)
text(z1,4,'\\Gamma','FontSize',n)

ShowText

% Script File: ShowText

close all
r = 1;
t = pi/6 + linspace(0,2*pi,7);
x = r*cos(t);
y = r*sin(t);
plot(x,y);

axis equal off
HA = 'HorizontalAlignment';
VA = 'VerticalAlignment';
text(x(1),y(1),'\leftarrow {\itP}_{1}',  HA,'left','FontSize',14)
text(x(2),y(2),'\downarrow',             HA,'center',VA,'baseline','FontSize',14)
text(x(2),y(2),'{  \itP}_{2}',           HA,'left',VA,'bottom','FontSize',14)
text(x(3),y(3),'{\itP}_{3} \rightarrow', HA,'right','FontSize',14)
text(x(4),y(4),'{\itP}_{4} \rightarrow', HA,'right','FontSize',14)
text(x(5),y(5),'\uparrow',               HA,'center',VA,'top','FontSize',14)
text(x(5),y(5),'{\itP}_{5}  ',           HA,'right',VA,'top','FontSize',14)
text(x(6),y(6),'\leftarrow {\itP}_{6}',  HA,'left','FontSize',14)

text(0,1.4*r,'A Labeled Hexagon^{1}',HA,'center','FontSize',14)
text(0,-1.4*r,'^{1} A hexagon has six sides.',HA,'center','FontSize',10)

ShowLineWidth

% Script File: ShowLineWidth

close all
plot([0 14 14 0 0],[0 0 14 14 0])
text(7,13,'LineWidth','FontSize',18,'HorizontalAlign','center')
axis off
hold on
for width=0:10
   h = plot([3 12],[11-width 11-width]);
   if width>0
      set(h,'LineWidth',width)
      text(1,11-width,sprintf('%3d',width),'FontSize',14)
   else
      text(1,11-width,'default','FontSize',14)

ShowAxes

% Script File: ShowAxes

F = 'Times-Roman'; n = 12;
close all
t = linspace(0,2*pi); c = cos(t); s = sin(t);
plot(c,s)
axis([-1.3 1.3,-1.3 1.3])
axis equal

title('The Circle ({\itx-a})^{2} + ({\ity-b})^{2} = {\itr}^{2}','FontName',F,'FontSize',n)
xlabel('x','FontName',F,'FontSize',n)
ylabel('y','FontName',F,'FontSize',n)
set(gca,'XTick',[-.5 0 .5])
set(gca,'YTick',[-.5 0 .5])
grid on

ShowLegend

% Script File: ShowLegend

close all
t = linspace(0,2);
axis([0 2 -1.5 1.5])
y1 = sin(t*pi);
y2 = cos(t*pi);
plot(t,y1,t,y2,'--',[0 .5 1 1.5 2],[0 0 0 0 0],'o')
set(gca,'XTick',[])
set(gca,'YTick',[0])
grid on
legend('sin(\pi t)','cos(\pi t)','roots',0)

ShowColor

% Script File: ShowColor

close all
plot([-1 16 16 -1 -1],[-3 -3 18 18 -3])
axis off
text(0.5,16.5,'Built-In Colors','FontSize',14)
text(9.5,16.5,'A Gradient','FontSize',14)


% The predefined colors.
x = [3 6 6 3 3];
y = [-.5 -.5 .5 .5 -.5];
hold on
fill(x,y,'c')
text(0,0,'cyan','FontSize',12)
fill(x,y+2,'m')
text(0,2,'magenta','FontSize',12)
fill(x,y+4,'y')
text(0,4,'yellow','FontSize',12)
fill(x,y+6,'r')
text(0,6,'red','FontSize',12)
fill(x,y+8,'g')
text(0,8,'green','FontSize',12)
fill(x,y+10,'b')
text(0,10,'blue','FontSize',12)
fill(x,y+12,'k')
text(0,12,'black','FontSize',12)
fill(x,y+14,'w')
text(0,14,'white','FontSize',12)

% Making your own colors.
a = -.5;
b = 14.5;
L = 11;
R = 15;
n = 100;
k=1;
for f = [.4 .6 .7 .8 .85 .90 .95 1] 
   color = [f 1 1];
   text(x(1)+5,y(1)+(2*(k-1)),sprintf('[ %4.2f , %1d , %1d ]',f,1,1),'VerticalAlign','bottom')
   fill(x+9,y+(2*(k-1)),color)
   k = k+1;
end
hold off

MyExpF

  function y = MyExpF(x,n)
% y = MyExpF(x,n)
% x is a scalar, n is a positive integer
% and y = n-th order Taylor approximation to exp(x).
y = 1;
term = 1;
for k = 1:n
   term = x*term/k;
   y = y + term;
end

MyExp1

  function y = MyExp1(x)
% y = MyExp1(x)
% x is a column n-vector and y is a column n-vector with the property that
% y(i) is a Taylor approximation to exp(x(i)) for i=1:n.

  nTerms = 17;   
  n = length(x);
  y = ones(n,1);
  for i=1:n
     y(i) = MyExpF(x(i),nTerms);
  end

MyExp2

  function y = MyExp2(x)
% y = MyExp2(x)
% x is an n-vector and y is an n-vector with with the same shape
% and the property that y(i) is a Taylor approximation to 
% exp(x(i)) for i=1:n.

  y = ones(size(x));
  nTerms = 17;
  term = ones(size(x));
  for k=1:nTerms
     term = x.*term/k;
     y = y + term;
  end

MyExpW

  function y = MyExpW(x)
% y = MyExpW(x)
% x is a scalar and y is a Taylor approximation to exp(x).

y = 0;
term = 1;
k=0;
while abs(term) > eps*abs(y)
   k = k + 1;
   y = y + term;
   term = x*term/k;
end

MyExp3

  function y = MyExp3(x)
% y = MyExp3(x)
% x is a column n-vector and y is a column n-vector with the property that
% y(i) is a Taylor approximation to exp(x(i)) for i=1:n.

  n = length(x);
  y = ones(n,1);
  for i=1:n
     y(i) = MyExpW(x(i));
  end

MyExp4

  function y = MyExp4(x)
% y = MyExp4(x)
% x is an n-vector and y is an n-vector with  the same shape and the 
% property that y(i) is a Taylor approximation to exp(x(i)) for i=1:n.

  y = zeros(size(x));
  term = ones(size(x));
  k = 0;
  while any(abs(term) > eps*abs(y))
     y = y + term;
     k = k+1;
     term = x.*term/k;
  end

Derivative

  function [d,err] = Derivative(fname,a,delta,M2)
% d = Derivative(fname,a,delta,M2);
% fname  a string that names a function f(x) whose derivative
% at x = a is sought. delta is the absolute error associated with
% an f-evaluation and M2 is an estimate of the second derivative
% magnitude near a. d is an approximation to f'(a) and err is an estimate 
% of its absolute error.
%
% Usage:
%     [d,err] =  Derivative(fname,a)   
%     [d,err] =  Derivative(fname,a,delta) 
%     [d,err] =  Derivative(fname,a,delta,M2)

if nargin &= 3
   % No derivative bound supplied, so assume the
   % second derivative bound is 1.
   M2 = 1;
end
if nargin == 2
   % No function evaluation error supplied, so
   % set delta to eps.
   delta = eps;
end
% Compute optimum h and divided difference
hopt = 2*sqrt(delta/M2);
d   = (feval(fname,a+hopt) - feval(fname,a))/hopt;
err = 2*sqrt(delta*M2);

Sin10

  function y = sin10(x)
y = sin(10*x);

Represent

  function f = Represent(x)
% f = Represent(x)
% Yields a 3-digit mantissa floating point representation of f:
%
%    f.mSignBit   mantissa sign bit (0 if x>=0, 1 otherwise)
%    f.m          mantissa (= f.m(1)/10 + f.m(2)/100 + f.m(3)/1000)
%    f.eSignBit   the exponent sign bit (0 if exponent nonnegative, 1 otherwise)
%    f.e          the exponent (-9<=f.e<=9)
%
% If x is out side of [-.999*10^9,.999*10^9], f.m is set to inf.
% If x is in the range [-.100*10^-9,.100*10^-9] f is the representation of zero
% in which both sign bits are 0, e is zero, and m = [0 0 0].

f = struct('mSignBit',0,'m',[0 0 0],'eSignBit',0,'e',0);

if abs(x)<.100*10^-9 
   % Underflow. Return 0
   return
end

if x>.999*10^9
   % Overflow. Return inf
   f.m = inf;
   return
end
if x<-.999*10^9
   % Overflow. Return -inf
   f.mSignBit = 1;
   f.m = inf;
   return
end
  
% Set the mantissa sign bit
if x>0
   f.mSignBit = 0;
else
   f.mSignBit = 1;
end

% Determine m and e so .1 <= m < 1 and abs(x) = m*10^e 
m = abs(x); e = 0;
while m >= 1,  m = m/10; e = e+1; end
while m < 1/10,m = 10*m; e = e-1; end

% Determine nearest integer to 1000m 
z = round(1000*m);
% Set the mantissa  
f.m(1) = floor(z/100);
f.m(2) = floor((z - f.m(1)*100)/10);
f.m(3) = z - f.m(1)*100 - f.m(2)*10;
% Set the exponent and its sign bit.
if e>=0
   f.eSignBit = 0;
   f.e = e;
else
   f.eSignBit = 1;
   f.e = -e;
end

Convert

  function x = Convert(f)
% x = Convert(f)
% f is a is a representation of a 3-digit floating point number. (For details
% type help represent. x is the value of v.

% Overflow situations
if (f.m == inf) & (f.mSignBit==0)
   x = inf;
   return
end
if (f.m == inf) & (f.mSignBit==1)
   x = -inf;
   return
end

% Mantissa value
mValue = (100*f.m(1) + 10*f.m(2) + f.m(3))/1000;
if f.mSignBit==1
   mValue = -mValue;
end

% Exponent value
eValue = f.e;
if f.eSignBit==1
   eValue = -eValue;
end

x = mValue * 10^eValue;

Float

  function z = Float(x,y,op)
% z = Float(x,y,op)
% x and y are representations of a 3-digit floating point number. (For details
% type help represent.
% op is one of the strings '+', '-', '*', or '/'.
% z is the 3-digit floating point representation of x op y.

sx = num2str(convert(x)); 
sy = num2str(convert(y)); 
z = represent(eval(['(' sx ')' op '(' sy ')' ]));

Pretty

  function s = Pretty(f)
% s = Pretty(f)
% f is a is a representation of a 3-digit floating point number. (For details
% type help represent.
% s is a string so that disp(s) "pretty prints" the value of f.

if (f.m == inf) & (f.mSignBit==0)
   s = 'inf';
elseif (f.m==inf) & (f.mSignBit==1)
   s = '-inf';
else
   % Convert the mantissa.
   m = ['.' num2str(f.m(1)) num2str(f.m(2)) num2str(f.m(3))];
   if f.mSignBit == 1
      m = ['-' m];
   end
   % Convert the exponent.
   e = num2str(f.e);
   if f.eSignBit == 0
      e = ['+' num2str(f.e)];
   else
      e = ['-' num2str(f.e)];
   end
   s = [m ' x 10^' e ];
end

PadeArray

  function P = PadeArray(m,n)
% P = PadeArray(m,n)
% m and n are nonnegative integers.
% P is an (m+1)-by-(n+1) cell array.
%
% P{i,j} represents the (i-1,j-1) Pade approximation N(x)/D(x) to exp(x).

P = cell(m+1,n+1);
for i=1:m+1
   for j=1:n+1
      P{i,j} = PadeCoeff(i-1,j-1);
   end
end
  

PadeCoeff

  function R = PadeCoeff(p,q)
% R = PadeCoeff(p,q)
% p and q are nonnegative integers and R is a representation of the
% (p,q)-Pade approximation N(x)/D(x) to exp(x):
%
%     R.num  is a row (p+1)-vector whose entries are the coefficients of the 
%            p-degree numerator polynomial N(x).
%
%     R.den  is a row (q+1)-vector whose entries are the coefficients of the 
%            q-degree denominator polynomial D(x).
%
% Thus, 
%                  R.num(1) + R.num(2)x + R.num(3)x^2
%                 ------------------------------------
%                          R.den(1) + R.den(2)x
%
% is the (2,1) Pade approximation.

a(1) = 1; for k=1:p, a(k+1) =  a(k)*(p-k+1)/(k*(p+q-k+1)); end
b(1) = 1; for k=1:q, b(k+1) = -b(k)*(q-k+1)/(k*(p+q-k+1)); end
R = struct('num',a,'den',b);