Chapter 1 M-Files

M-File Home

Script File Index

 SineTable Prints a short table of sine evaluations. SinePlot Displays a sequence of sin(x) plots. ExpPlot Plots exp(x) and an approximation to exp(x). TangentPlot Plots tan(x). SineAndCosPlot Superimposes plots of sin(x) and cos(x). Polygons Displays nine regular polygons, one per window. SumOfSines Displays the sum of four sine functions. SumOfSines2 Displays a pair of sum-of-sine functions. UpDown Sample core exploratory environment. RunUpDown Framework for running UpDown. Histograms Displays the distribution of rand and randn. Clouds Displays 2-dimensional rand and randn. Dice Histogram of 1000 dice rolls. Darts Monte Carlo computation of pi. Smooth Polygon smoothing. Stirling Relative and absolute error in Stirling formula. ExpTaylor Plots relative error in Taylor approximation to exp(x). Zoom Roundoff in the expansion of (x-1)^ 6. FpFacts Examines precision, overflow, and underflow. TestMyExp Examines MyExp1 , MyExp2 , MyExp3 , and MyExp4. TestDerivative Applies Derivative to Sin10. Euler Three-digit arithmetic sum of 1 + 1/2 + ... + 1/n. ShowPadeArray Tests the function PadeArray. ShowFonts Illustrates how to use fonts. ShowSymbols Shows how to generate math symbols. ShowGreek Shows how to generate Greek letters. ShowText Shows how to align with text. ShowLineWidth Shows how vary line width in a plot. ShowAxes Shows how to set tick marks on axes. ShowLegend Shows how to add a legend to a plot. ShowColor Shows how to use built-in colors and user-defined colors.

Function File Index

 MyExpF For-loop Taylor approximation to exp(x). MyExp1 Vectorized version of MyExpF. MyExp2 Better vectorized version of MyExpF. MyExpW While-loop Taylor approximation to exp(x). MyExp3 Vectorized version of MyExpW. MyExp4 Better vectorized version of MyExpW. Derivative Numerical differentiation. Sin10 sin(10x). Represent Sets up 3-digit arithmetic representation. Convert Converts 3-digit representation to float. Float Simulates 3-digit arithmetic. Pretty Pretty prints a 3-digit representation. PadeArray Builds a cell array of Pade coefficients. PadeCoeff Generates Pade approximant coefficients

% 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


% 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


% 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)],'*')


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


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


% 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


% 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


% 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


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


% 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


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


% 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


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


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


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

clc
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


% 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


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


% 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


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

% 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


  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


  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


  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


  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


  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;


  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


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




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