Solutions to Prelim 2 Review Questions
--------------------------------------

%-------------------------------------
% Insight Ch4 Sec1 #5
% Solicits input for values a, b,c and prints the length of a:c:b and the
% distance from the last vector component to b.

a = input('a = ');
b = input('b = ');
c = input('c = ');
v = a:c:b;

fprintf('The length of a:c:b is %d.\n', length(v))
fprintf('The distance from v(end) to b is %f.\n', v(end)-b)

% Note:  When used as an index, the keyword  end  is the same as the 
% length of the vector.  So the distance from the last vector component to 
% b can be written as   v(length(v)) - b


%-------------------------------------
% Insight Ch4 Sec1 #6

x = 1:10:100;
while length(x) < 1000
  x = [x x(1) x];
end

fprintf('length(x) = %d.\n',length(x));

% Note: If the length of x in iteration k is N, then after iteration k+1
% the length of x is 2*N + 1. At the start length(x) = 10 so as the loop
% progresses the length of x is
% 10, 21, 43, 87, 175, 351, 703, 1407. 


%-------------------------------------
% Insight Ch4 Sec1 #11
% Generate a length-10 vector according to a recursion formula.

n = input('n = ');

f = zeros(1,10);
for k=1:10
  if k==1                  % Starting value
    f(k) = n;
  elseif rem(f(k-1),2)==1  % f(k-1) is odd
    f(k) = 3*f(k-1) + 1;
  else                     % f(k-1) is even
    f(k) = f(k-1)/2;
  end
end

plot(1:10,f,'*')           % Plot the result using the star-marker


%-------------------------------------
% Insight Ch4 Sec3 #3
% Evaluate error
%   e_h(x) = abs( (sin(x+h)-sin(x))/h - cos(x) ) = big_O(h)

clc
x = input('Enter a real number in [0,2pi]: ');

n= 16;
h= 1./logspace(1,n,n);
e_h= zeros(1,n);
e_best= Inf;

disp('  h          error ')
for k= 1:n
    e_h(k)= abs( (sin(x+h(k))-sin(x))/h(k) - cos(x) );
    fprintf('%.0e    %.4e \n', h(k), e_h(k))
    if e_h(k)<e_best
        e_best= e_h(k);
        h_best= h(k);
    end
end

fprintf('\n Optimal h: %.0e \n', h_best)
fprintf(' Error: %.4e \n', e_best)


%-------------------------------------
% Insight Ch5 Sec1 #9

function Snm = Subset(n, m)
% Snm is the number ways a set of n objects can be partitioned into m nonempty subsets.
% n and m are positive integers.

total= 0;
sgn= 1;        % sign of the term
for j= m:-1:1  % count down since the mth term has positive sign (sgn is 1)
    total= total + sgn * j^n / factorial(m-j) / factorial(j);
    sgn= -sgn;
end
Snm= total;

% Note:  Due to round-off errors, Snm may not be an integer.  To avoid this
% problem, one can calculate the summation of   Snm * m!  instead, i.e., 
% every term in the summumation is multiplied by m!.
% This means you actually calculate the summation from j=1 to j=m of
%    (-1)^m-j * j^n * m! / ((m-j)! * j!)
%                     ^^^^^^^^^^^^^^^^^^
%                     This part is the binomial coefficient B(m,j)
%                     which can be calculated as   bchoosek(m,j).
%                     Refer to problems P5.1.6 and P5.1.7.
%
% Below is the alternative solution:
%
% total = 0;
% sgn = 1;    % sign of the term
% for j=m:-1:1
%     total = total + sgn * nchoosek(m, j) * j^n;
%     sgn = -sgn;
% end
% Snm = total / factorial(m);


%-------------------------------------
% Insight Ch5 Sec1 #9, cont'd
% Print the table of Sigma(n,m)

% print the header line
disp('  n    S_n^(1)    S_n^(2)    S_n^(3)    S_n^(4)    S_n^(5)    S_n^(6)')
disp('---------------------------------------------------------------------')

% print the values
for n=1:10
    % string of the line to be printed
    line = sprintf('%3d  ', n);
    for m=1:6
        val= Subset(n,m);  % val may be non-integer due to round-off error 
        val= round(val);
        line = sprintf('%s%6d     ', line, val);
    end
    disp(line)
end


%-------------------------------------
% Insight Ch6 Sec1 #15
% Simulate selecting two points on the unit circle (uniformly random)
% and calculate the probability that the connecting chord is longer than 1.

N = 10000; % number of trials
num = 0;   % number of times that the the connecting chord is longer than 1
for k=1:N
    % generate the random points on the unit circle
    theta1 = rand(1) * 2*pi;
    theta2 = rand(1) * 2*pi;
    cx1 = cos(theta1);
    cy1 = sin(theta1);
    cx2 = cos(theta2);
    cy2 = sin(theta2);
    
    % test if the connecting chord is longer than 1
    if (cx1 - cx2)^2 + (cy1 - cy2)^2 > 1
        num = num + 1;
    end
end
fprintf('Probability that connecting chord is longer than 1 is %f\n', num/N)


%-------------------------------------
% Insight Ch6 Sec1 #19a

function p = ProbG(L,R)
% p is the area under the Gaussian function (bell curve) between L and R

% initialize...
N= 10000;   % Number of trials
count = 0;  % Number of points under the curve

% count the number of random points that lie beneath the curve
for i = 1:N
    % Random point in the rectangle bounded by x=L, x=R, y=0, and y=1
    x= rand(1)*(R-L)+L;
    y= rand(1);

    % If the point is under the curve, i.e. f(x), then count as a "hit."
    if y <= (2*pi)^(-1/2) * exp(-x^2/2)
        count=count+1;
    end
end

% Monte Carlo estimate of area under the curve
p= (count/N)*(R-L);


%-------------------------------------
% Insight Ch6 Sec2 #8

function [x, y, z] = RandomWalk3D(n)
% Random walk in 3D.  x,y,z are vectors representing the path.
% The kth step has the coordinates (x(k),y(k),z(k)).

%initialize coordinates and counter
xc=0; yc=0; zc=0; k=0;

%while not on the boundary, take a random step
while abs(xc)<n && abs(yc)<n && abs(zc)<n
    r=rand(1);
    if r <= 1/6
        xc = xc + 1;
    elseif r <= 2/6
        xc = xc - 1;
    elseif r <= 3/6
        yc = yc + 1;
    elseif r <= 4/6
        yc = yc - 1;
    elseif r <= 5/6
        zc = zc + 1;
    else
        zc = zc - 1;
    end

    %update position arrays
    k=k+1; x(k)=xc; y(k)=yc; z(k)=zc;
end


%-------------------------------------
% Insight Ch6 Sec2 #8, cont'd
% Estimate the average number of hops required for the robot to
% reach the boundary in a 3D random walk.

clc
disp('   n    Average #steps to Boundary')
disp('----------------------------------')
nTrials = 100;
for n = 5:5:50
    s = 0;
    for k=1:nTrials
        [x,y,z] = RandomWalk3D(n);
        s = s + length(x);
    end
    ave = s/nTrials;
    fprintf(' %3d             %8.3f\n',n,ave)
end
fprintf('\n\n(Results based on %d trials)\n',nTrials)


%-------------------------------------
% Review Question 1

function  seconds = hms2s(h, m, s)
% Convert a time expressed in hours, minutes, seconds to a time in seconds.
%   h2ms( h, m, s )  returns  3600*h + 60*m + s

seconds = 3600*h+60*m+s;


%-------------------------------------
% Review Question 1, cont'd

function [h, m, s] = s2hms( seconds )
% Convert a time in seconds to a time in hrs, minutes, seconds.
% seconds = total number of seconds
% h = maximum number of hours in seconds
% m = maximum number of minutes in seconds, after max hours are removed
% s = number of seconds remaining after max hours and max minutes removed

h = floor( seconds/3600 );
m = floor( rem( seconds, 3600 )/60 );
s = rem( rem(seconds, 3600), 60 );


%-------------------------------------
% Review Question 1, cont'd
% Script to demonstrate functions hms2s and s2hms

sec = hms2s( 7, 45, 13 )
[h, m, s ] = s2hms( 3682 )


%-------------------------------------
% Review Question 2
% Script to plot f(x)=0.2*(x^3-x)/(1.1+cos(x)) for -2pi <= x <= 2pi

x = linspace(-2*pi, 2*pi, 200);
y = 0.2*(x.^3-x)./(1.1+cos(x));
plot(x,y)