Matlab Lecture Examples From Week 2

% eg1
theta = linspace(0,2*pi);
c = cos(theta); s = sin(theta);
plot(c,s)

% eg2
theta = linspace(0,2*pi);
c = cos(theta); s = sin(theta);
plot(c,s)
axis equal

% eg3
theta = linspace(0,2*pi);
c = cos(theta); s = sin(theta);
plot(c,s)
axis equal
axis([-1.2 1.2 -1.2 1.2])

% eg4
theta = linspace(0,2*pi);
c = cos(theta); s = sin(theta);
plot(c,s)
axis equal
hold on
plot(c+1,s+2)
plot(c+2,s-1)
hold off

% eg5
theta = linspace(0,2*pi);
c = cos(theta); s = sin(theta);
plot(c,s,c+1,s+2,c+2,s-1)
axis equal

% eg6
x = linspace(0,2*pi)';
y1 = sin(x); y2 = sin(2*x); y3 = sin(3*x); y4 = sin(4*x);
A = [y1 y2 y3 y4];
for k=1:4
   subplot(2,2,k)
   plot(x,A(:,k))
   axis([0,2*pi,-1.1 1.1])
   title(sprintf('sin(%1dx)',k))
end

% eg7
x = linspace(0,2*pi)';
y1 = sin(x); y2 = sin(2*x); y3 = sin(3*x); y4 = sin(4*x);
A = [y1 y2 y3 y4];
for k=1:4
   subplot(2,2,k)
   plot(x,A(:,k))
   axis([0,2*pi,-1.1 1.1])
   title(sprintf('sin(%1dx)',k))
end
figure
c = [1;-2;-3;2];
y = A*c;
plot(x,y)
axis([0 2*pi -6 6])
title('f(x) = sin(x) - 2sin(2x) -3sin(3x) + 2sin(4x)')


% eg8
h = logspace(-1,-16,16);
x = pi/3;
dsin = (sin(x+h) - sin(x))./h;
error = abs(dsin - cos(x));
loglog(h,error,h,error,'o')
xlabel('h')
ylabel('|dsin(x) - cos(x)|')


% Smallest positive integer p so 1 + 1/2^p = 1 in floating point.
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 = %2d',p))

% Smallest positive integer q so 1/2^q = 0 in floating point.
x = 1;
q = 0;
while (x>0)
   x = x/2;
   q = q+1;
end
disp(sprintf('q = %2d',q))


% Smallest positive integer r so 2^r = inf in floating point.
x = 1;
r = 0;
while (x~=inf)
   x = 2*x;
   r = r+1;
end
disp(sprintf('r = %2d',r))


% Warm-Up Problem.
% What is the output...
h = 1./2.;
x = 2./3. - h;                 % = 1/6
y = 3./5. - h;                 % = 1/10
e = (x+x+x) - h;               % = 3/6 - 1/2  = 0
f = (y+y+y+y+y)-h;             % = 5/10 - 1/2 = 0
z = f/e


% Error associated with MySin0
n = 200;
x = linspace(0,4*pi,n);
y = zeros(1,n);
for k=1:n
   y(k) = MySin0(x(k));
end
error = abs(y-sin(x));
semilogy(x,error)
title('|MySin0(x)  -  sin(x)|')

function y = MySin0(x)
% y = sin(x) where x is a scalar
s = 0; 
term = x; 
k = 0;
z = -x^2;
while(s+term~=s)
   s = s+term;
   term = z*term/((2*k+2)*(2*k+3));
   k = k+1;
end
y = s;


% Error associated with MySin1
n = 200;
x = linspace(0,4*pi,n);
y = MySin1(x);
error = abs(y-sin(x));
semilogy(x,error)
title('|MySin1(x)  -  sin(x)|')

function y = MySin1(x)
% y = sin(x) where x is a vector
s = zeros(size(x)); 
term = x; 
k = 0;
z = -x.^2;
while(any(s+term~=s))
   s = s+term;
   term = z.*term/((2*k+2)*(2*k+3));
   k = k+1;
end
y = s;


% Error associated with MySin2
n = 200;
x = linspace(0,4*pi,n);
y = MySin2(x);
error = abs(y-sin(x));
semilogy(x,error)
title('|MySin2(x)  -  sin(x)|')

function y = MySin2(x)
% y = sin(x) where x is a vector
y = MySin1(rem(x,2*pi));


function SineTerms(x)
% Plots the the terms in the sine series used by MySin0
close all 
figure
s = 0; term = x; k = 0; z = -x^2;
tau = [x];
while(s+term~=s)
   s = s+term;
   term = z*term/((2*k+2)*(2*k+3));
   k = k+1;
   tau = [tau term];
end
subplot(2,1,1)
plot(tau)
title(' (-1)^k x^{2k+1}/(2k+1)!')
xlabel('k')
subplot(2,1,2)
semilogy(abs(tau))


% Benchmarking script
n = 3000; 
x = linspace(0,2*pi,n);
tic 
y = zeros(1,n);
for k=1:n
   y(k) = MySin0(x(k));
end
tMySin0 = toc;

nTrials = 100;
tic
for mu = 1:nTrials
   y = MySin2(x);
end
tMySin2 = toc/nTrials;

nTrials = 100;
tic
for mu = 1:nTrials
   y = sin(x);
end
tSin = toc/nTrials;