Derivative.m    Polygons.m      deriv_err.m     plots.m         sorte.m*
Operations      RunSort.m       derr.m          polygons.m
PadeArray.m     TangentPlot.m   floats.m        powers.m
PadeCoeff.m     cells.m         pade_driver.m   sort.m
%-----------------------------------------------------------------------
%Derivative.m
function [d,err] = Derivative(fname,a,delta,M2)
%  [d,err] = Derivative(fname,a,delta,M2)
%  fname is 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 is supplied; assume the 2nd derivative bound is 1.
   M2 = 1
end
if nargin == 2
   % no function evaluation errors supplied, so
   % set delta to eps.
   delta = eps;
end
%  Plot error (and it's derivative) vs. h to exhibit optimal step
x = linspace(1000*eps,10000000000*eps,250);
   y = (M2/2)*x+2*delta./x;
plot(x,y)
%  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);
%-----------------------------------------------------------------------
  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,n);
for i=1:m
   for j=1:n
      P{i,j} = PadeCoeff(i,j);
   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); 
%-----------------------------------------------------------------------
  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.
p = p-1; q = q-1;
   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); 
%-----------------------------------------------------------------------
% 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 equal
end  
%-----------------------------------------------------------------------
% Script File: RunSort
while(input('Another Example? (1=yes, 0=no)'))
   diary SortOut
   sorte
   diary off
   if (input('Keep output? (1=yes, 0=no)')~=1)
      delete SortOut
   end
end 
%-----------------------------------------------------------------------
% Script File: TangentPlot
% Displays increasingly smooth plots of Tangent(2*pi*x)
% in the range -pi/2 <= x <= 9pi/2
close all
% n=40
for n =[4 8 12 16 20 50 100 200 400 ]
% for n =4:4:36
  x=linspace(-pi/2,pi/2,n); y = tan(x);
   z=cos(x); plot(x,y,x,z)
  ymax = 4;
  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)
%   plot(xnew,y,'*')
    z = cos(xnew);
%   plot(xnew,z)
    plot(xnew,y,'*',xnew,z,'o')
  title(' The cos-Tangent functions')
  end
 hold off
  pause(1)
end
%-----------------------------------------------------------------------
% Script: cells
NYC_Lat  = struct('d',40,'m',45,'s',27);
NYC_Long = struct('d',75,'m',12,'s',32);
C1       = struct('name','New York','lat',NYC_Lat,'long',NYC_Long);
C1.name
C1.long.m
C=cell(2,2);
C{1,1}=[1 2;3 4];
C{1,2}=[5;6];
C{2,1}=[7 8];
C{2,2}=9;
M=[C{1,1} C{1,2};C{2,1} C{2,2}] 
%-----------------------------------------------------------------------
% script: deriv_err
%------------------
%Problem 1 (1.4.10)
%------------------
%  Let h be a small number. The derivative of a function f at x_0 can be
% approximated by a forward divided difference:
%                         f(x_0 + h) - f(x_0)
%                f'(x) = ---------------------
%                                  h
% and by a central divided difference:
%                         f(x_0 + h) - f(x_0 - h)
%                f'(x) = ------------------------- .
%                                   2h
%  For the function f(x) = sin(x), plot the error when these approximations are
% used to estimate f'(1) = cos(1) for h = 10^-1, 10^-2,..., 10^-16. Repeat
% with x_0 = 1 + 10^6 pi. Use loglog for each of the four plots and display
% them all in the same window using subplot. Be sure to title each plot and to
% label the axes appropriately. Your script should not have any loops.
%------------------
%Theory:
% the truncation error for the forward difference approximation is found in
%the text (p. 56): 
%              |D_h - f'(x_0)|<= M_2 h/2 + 2 d/h
% and the truncation error for the central difference approximation is given in
%problem 1.5.1 (p. 57):
%              |C_h - f'(x_0)|<= M_3 h^2/6 + 2 d/h
%              |C_h - f'(x_0)|<= M_3 h^2/6 + 2 d/h
%------------------
% Begining of the script
h = [.1e0 .1e-1 .1e-2 .1e-3 .1e-4 .1e-5 .1e-6 .1e-7 .1e-8 .1e-9 .1e-10 ...
      .1e-11 .1e-12 .1e-13 .1e-14 .1e-15];
a=1+10^6 * pi;
x = [1;a];
x0 = x*ones(1,16);
ho = ones(2,1)*h;
xp = x0+ho;
xm = x0-ho;
D  = (sin(xp)-sin(x0))./ho;
C  = (sin(xp)-sin(xm))./(2*ho);
exact = cos(x0);
DRelErr = abs(exact - D)./exact;
CRelErr = abs(exact - C)./exact;
for i=1:2   
      subplot(2,2,i)
   loglog(h,DRelErr(i,:))
   title(sprintf('Forward difference; x = %2.0f',x(i)))
end
for i=1:2
   subplot(2,2,i+2)
   loglog(h,CRelErr(i,:))
   title(sprintf('Centered difference; x = %2.0f',x(i)))
end
%  subplot(2,2,1)
%  loglog(h,DRelErr(1,:))
%  title(sprintf('Forward difference; x = %2.0f',x(1)))
%  subplot(2,2,2)
%  loglog(h,DRelErr(2,:))
%  title(sprintf('Forward difference; x = %2.0f',x(2)))
%  subplot(2,2,3)
%  loglog(h,CRelErr(1,:))
%  title(sprintf('Centered difference; x = %2.0f',x(1)))
%  subplot(2,2,4)
%  loglog(h,CRelErr(2,:))
%  title(sprintf('Centered difference; x = %2.0f',x(2)))   
%-----------------------------------------------------------------------
% Script File: derr
%  Plot error (and it's derivative) vs. h to exhibit optimal step
clear
fname = 'sin';
fname1= 'cos';
delta = eps;
hmin  = 10^6*eps;
hmax  = 10^10*eps;
xmax  = 1;
x=pi/4;
M2=1;
[d,err,hopt] = Derivative(fname,x,delta,M2);
yopt = (M2/2)*hopt+2*delta/hopt;
   h = linspace(hmin,hmax,1000);
   y = (M2/2)*h+2*delta./h;
   y1=   (M2/2)*h;
   y2= 2*delta./h;
   derivh =(feval(fname,x+h) - feval(fname,x))./h;
   abserr=abs(feval(fname1,x)-derivh);
   loglog(h,y,h,y1,h,y2,h,abserr,hopt,yopt,'bo');
   xlabel(sprintf('hopt = %3.2d, d%s(%3.2d)/dx=%3.2d +/-%3.2d',hopt,...
      fname ,x,d,err))
   ylabel('abserr, maxerr')
   title(sprintf('derivative error as function of h for %s(x)',fname))
   axis([hmin hmax 10^(-12) 10^(-4)])
%-----------------------------------------------------------------------
% Script: floats
% find machine epsilon
x = 1; p = 0; y = 1; z = x + y;
while x ~= z, y = y/2; p = p+1; z = x+y; end
disp(sprintf('Machine epsilon = %22.15e and exponent = %1.0f',2*y,p))
% find underflow exponent
x = 1; q = 0;
while x > 0, x = x/2; q = q+1; z = x+y; end
disp(sprintf('underflow exponent: base 2 = %1.0f, base 10 = %1.0f',-q/2,-q/2/log(10)))
x = 1; r = 0;
while x ~= inf, x = x*2; r = r+1; z = x+y; end
disp(sprintf('overflow exponent: base 2 = %1.0f, base 10 = %1.0f',r/2,r/2/log(10)))  
%-----------------------------------------------------------------------
% Script File: pade_driver
clear
exp(1)
x = 1;
P = PadeArray(8,8);
for p=1:8
   for q=1:8
      numer = 0;
      for i=1:p
        numer = numer*x+P{p,q}.num(p-i+1);
     end 
     denom = 0;
      for j=1:q
         denom = denom*x+P{p,q}.den(q-j+1);
      end
      expr(p,q) = numer/denom;
   end
end
expr-exp(1)
%-----------------------------------------------------------------------
% Script File: plots
m = 50;
x = linspace(-1,1,50);
y = zeros(1,m);
exact = exp(x);
k = 0;
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(3,2,k)
   plot(x,RelErr)
   title(sprintf('n = %2.0f',n))
end
RelErr = abs(exact - MyExp1(x')')./exact;
   subplot(3,2,5)
   plot(x,RelErr)
   title(sprintf('MyExp1 '))
RelErr = abs(exact - MyExp2(x))./exact;
   subplot(3,2,6)
   plot(x,RelErr)
   title(sprintf('MyExp2 '))
%-----------------------------------------------------------------------
% Script File: Polygons
% Plots selected regular polygons
close all
theta = linspace(0,2*pi,360);
c=cos(theta); s = sin(theta);
k=0;
for sides = [3 4 5 6 7 10 12 18 24]
   stride = 360/sides;
   k = k+1;
   subplot(3,3,k)
   plot(c(1:stride:360),s(1:stride:360))
   axis equal
end
%-----------------------------------------------------------------------
% Script File: powers
% Problem 1.2.5, p.28
% Superimposes plots of x^k, k = 1:10
% in the range 0 <= x <= 1
close all
m=10
x=linspace(0,1,100)'; y = x; plot(x,y)
axis([0 1 0 1])
  title(' The Power function')
  xlabel(' x ')
  ylabel(' x^k ; k=1,...,10 ')
hold on
for k = 2:m
   y=y.*x; plot(x,y)
end
hold off 
% Script: sort
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);
[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));
if title('Local Maxima'),I+1,x(I+1),I+1,x(I+1),'*')
end % Script: sorte
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);
[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));
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 
