ExpPlot.m       MyExp2.m        SinePlot.m      gasket.m        taylor2.m
ExpTaylor.m     MyExp3.m        SineTable.m     gasket_driver.m taylorFlops.m
ExpTaylorerr.m* MyExp4.m        Stirling.m      lorenzde.m      taylors.m
Lorenz.m*       MyExpF.m        TangentPlot.m   mandel.m        taylorv.m
MyExp1.m        MyExpW.m        fourier.m       polygons.m      experr.m
                                                                Zoom.m
%-----------------------------------------------------------------------
% Script File: ExpPlot
% Examine the function f(x)=((1+x/24)(1-x/12+x^2/384))^8
% as an approximation to exp(x) across [0,1]
close all
x = linspace(0,50,1000);
num = 1+x/24;
denom = 1 - x/12 + (x/384).*x;
quot = num./denom;
y = quot.^8;
plot(x,y',x,exp(x),'o')
  axis([0 50 0 1000000])
%-----------------------------------------------------------------------
% 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]
x = input('Enter x: ');
   figure
   term = 1; s = 1; f = exp(x)*ones(nTerms,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)
%  semilogy(1:nTerms,err)
   ylabel('Relative Error in Partial Sum.')
   xlabel('Order of Partial Sum.')
   title(sprintf('x = %5.2f',x))
   figure
   semilogy(1:nTerms,err)
%end
%-----------------------------------------------------------------------
% Script File: ExpTaylorerr.m
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 - MyExp4(x))./exact;
   subplot(3,2,6)
   plot(x,RelErr)
   title(sprintf('MyExp4 '))
%-----------------------------------------------------------------------
% Script File: Lorenz
%LRUN ODE solving example: Lorenz
tspan = [0 50];                % solve for 0 <= t <= 50
yzero = [0;1;0];               % Initial conditions
[t,y] = ode45('lorenzde',tspan,yzero);
plot(y(:,1),y(:,3))
xlabel('y_1','FontSize',14)
ylabel('y_3','FontSize',14,'Rotation',0,'HorizontalAlignment','right')
title('Lorenz equations','FontSize',16
%-----------------------------------------------------------------------
  function y = MyExp1(x)
% y = MyExp1(x)
% x is a column vector of length p, n is a positive integer
% and y is a column vector with the property that y(i) = n-th order Taylor 
% approximation to exp(x(i)) for i = 1:p
n = 17; p = length(x);
y = ones(p,1);
for i = 1:p
   y(i) = MyExpF(x(i),n);
end
%-----------------------------------------------------------------------
  function y = MyExp2(x)
% y = MyExp2(x)
% x is a p-vector, n is a positive integer
% and y is a p-vector of the same shape as x and with the property 
% that y(i) = nTerms-th order Taylor 
% approximation to exp(x(i)) for i = 1:p
% p = size(x)
nTerms = 17;
y    = ones(size(x));
term = ones(size(x));
for k = 1:nTerms
   term = x.*term/k;
   y    = y +term;
end
%-----------------------------------------------------------------------
 function y = MyExp3(x)
% y = MyExp3(x)
% x is a column vector of length p,
% and y is a column vector with the property that y(i) is a Taylor 
% approximation to exp(x(i)) for i = 1:p
p = length(x);
y = ones(p,1);
for i = 1:p
   y(i) = MyExpW(x(i));
end
%-----------------------------------------------------------------------
  function y = MyExp4(x)
% y = MyExp4(x)
% x is a p-vector
% y is a p-vector of the same shape as x
% and with the property that y(i) is the
% Taylor approximation to exp(x(i)) for 
% i = 1:p; p = size(x)
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
%-----------------------------------------------------------------------
  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)
term = 1; y = 1;
for k = 1:n
   term = x*term/k; y = y + term;
end
%-----------------------------------------------------------------------
  function y = MyExpW(x)
% y = MyExpW(x)
% x is a scalar  and y is Taylor approximation to exp(x)
term = 1; y = 0;
k = 0;
while abs(term) > eps*abs(y)
   k = k+1;
   y = y+term;
   term = x*term/k;
end
%-----------------------------------------------------------------------
% 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(' The function y = sin(2*pi*x) ')
  xlabel(' x (in rads) ')
  ylabel(' y ')
  pause(1)
end
%-----------------------------------------------------------------------
% Script: Sine Table
% Prints a short table of sine values
clc
n=30;
x=linspace(0,1,n);
y=sin(2*pi*x);
disp(' ')
   disp(        '  k        x(k)   sin(x(k))')
   disp(        '---------------------------------------')
for k=1:n
   degrees = (k-1)*360/(n-1);
   disp(sprintf(' %2.0f        %3.0f   %8.5f    ',k,degrees,y(k)))
end
disp(' ')
disp(' x(k) is given in degrees ')
disp(sprintf(' One Degree = %7.5e  radians ',pi/180))
%-----------------------------------------------------------------------
% Script File: Strirling
% Prints a table showing error in Stirling's formula for n!
clc
disp('                          Strirling       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         
%-----------------------------------------------------------------------
% 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 fourier
% 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 ;3 ;7 ;5 ];
y = A*[2 8 1;3 2 0;7 6 0;5 9 0];
plot(x,y)
%-----------------------------------------------------------------------
function gasket(Pa,Pb,Pc,level)
%GASKET Recursively generated Sierpinski gasket
%       GASKET(PA,PB,PC,LEVEL) generates an approximation to
%       the Sierpinski gasket, where the 2-vectors PA, PB and PC
%       define the triangle vertices.
%       LEVEL is the level of recursion
if level == 0
   % Fill the triangle with vertices Pa, Pb, Pc
   fill([Pa(1),Pb(1),Pc(1)],[Pa(2),Pb(2),Pc(2)],[0.5,0.5,0.5]);
   hold on
else
   % Recursive calls for the three subtriangles
   gasket(Pa,(Pa+Pb)/2,(Pa+Pc)/2,level-1)
   gasket(Pb,(Pb+Pa)/2,(Pb+Pc)/2,level-1)
   gasket(Pc,(Pc+Pa)/2,(Pc+Pb)/2,level-1)
end
%-----------------------------------------------------------------------
% Script File: gasket_driver
level = 8;
Pa=[0;0];
Pb=[1;0];
Pc=[0.5;sqrt(3)/2];
gasket(Pa,Pb,Pc,level)
hold off
title(['Gasket level =',num2str(level)],'FontSize',16)
axis('equal','off')
%-----------------------------------------------------------------------
function yprime = lorenzde(t,y)
%LORENZDE Lorenz Equation
% YPRIME = LORENZDE(T,Y)
yprime =[10*(y(2)-y(1));28*y(1)-y(2)-y(1)*y(3);y(1)*y(2)-8*y(3)/3];
%-----------------------------------------------------------------------
%script MANDEL   Mandelbrot Set
clc
close
h = waitbar(0,'Computing...');
x = linspace(-2.1,0.6,301);
y = linspace(-1.1,1.1,301);
[X,Y] = meshgrid(x,y);
%C = complex(X,Y);
C1 = X; C2 = Y;

Z_max = 1e6; it_max = 50;
Z1 = X; Z2 = Y;
for k = 1:it_max
   %Z = Z.^2 + C;
   Z1 = Z1.*Z1-Z2.*Z2 + X;
   Z2 = 2*Z1.*Z2 + Y;
   waitbar(k/it_max)
end
close(h)
Z = sqrt(Z1.^2+Z2.^2);
contourf(x,y,Z<Z_max,1)
title('Mandelbrot Set','FontSize',16)
%-----------------------------------------------------------------------
% 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 bench34
nRepeat = 100;
disp('T3 = MyExp3 benchmark.')
disp('T4 = MyExp4 benchmark.')
disp(' Length(x)     T4/T3')
disp('----------------------------------------------------')
for L = 1000:100:2000
   xL = linspace(-1,1,L);
   tic, for k=1:nRepeat, y = MyExp3(xL); end, T3 = toc;
   tic, for k=1:nRepeat, y = MyExp4(xL); end, T4 = toc;
   disp(sprintf('%6.0f  %13.6f  ',L,T4/T3))
end
disp('----------------------------------------------------')
for L = 1000:100:2000
% T3 = MyExp3 benchmark.
% T4 = MyExp4 benchmark.
%  Length(x)     T4/T3
% ----------------------------------------------------
%   1000       0.022783  
%   1100       0.022148  
%   1200       0.022003  
%   1300       0.021454  
%   1400       0.021335  
%   1500       0.021154  
%   1600       0.020791  
%   1700       0.020608  
%   1800       0.020592  
%   1900       0.020197  
%   2000       0.019923
% ----------------------------------------------------
%-----------------------------------------------------------------------
% script TaylorFlops
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 taylors
m = 50;
x = linspace(-1,1,m);
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(2,2,k)
   plot(x,RelErr)
   title(sprintf('n = %2.0f',n))
end
%-----------------------------------------------------------------------
% script taylorv
nRepeat = 100;
disp('T1 = MyExp1 benchmark.')
disp('T2 = MyExp2 benchmark.')
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
% T1 = MyExp1 benchmark.
% T2 = MyExp2 benchmark.
%  Length(x)     T2/T1
% ----------------------------------------------------
%   1000       0.012505  
%   1100       0.012224  
%   1200       0.012098  
%   1300       0.011663  
%   1400       0.011508  
%   1500       0.011173  
%   1600       0.011123  
%   1700       0.010729  
%   1800       0.009902  
%   1900       0.010098  
%   2000       0.010312  
%-----------------------------------------------------------------------
%script file experr
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 
z=MyExp1(x');
RelErr = abs(exact - MyExp1(x')')./exact;
   subplot(3,2,5)
   plot(x,RelErr)
   title(sprintf('MyExp1 '))
RelErr = abs(exact - MyExp4(x))./exact;
   subplot(3,2,6)
   plot(x,RelErr)
   title(sprintf('MyExp4 '))
%-----------------------------------------------------------------------
% Script File: Zoom
% Plots (x-1)^6 near x=1 with increasingly refined scale.
% Evaluation via binomial expansion leads to severe cancellation
% (x-1)^6 = x^6 - 6*x^5 + 15*x^4 - 20*x^3 + 15*x^2 -6*x + 1
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);
   z = ones(n,1) + x.*(-6 +x.*(15 +x.*(-20 +x.*(15 +x.*(-6+x)))));
   w = (x-1).^6;
   k = k+1; subplot(6,3,k); plot(x,y,x,zeros(1,n))
   axis([1-delta 1+delta -max(abs(y)) max(abs(y))])
   subplot(6,3,k+6); plot(x,z,x,zeros(1,n))
   axis([1-delta 1+delta -max(abs(z)) max(abs(z))])
   subplot(6,3,k+12); plot(x,w,x,zeros(1,n))
   axis([1-delta 1+delta -max(abs(w)) max(abs(w))])
end
