Solutions to Prelim 1 Review Questions QUESTION 2 ------------------------------------------------------ "M Questions" from Insight. You can check these on your own using Matlab. Ask a course staff member if you have questions. QUESTION 3 ------------------------------------------------------ amt = input('Enter the Income: '); if (amt > 60000) tax = (25/100)*30000+(30/100)*30000+(40/100)*(amt-60000); elseif (amt > 30000) tax = (25/100)*30000+(30/100)*(amt-30000); else tax = (25/100)*(amt); end fprintf('The tax to be paid for an income of %d is %d',amt,tax); QUESTION 4 ----------------------------------------------------- % Sec 2.1 Problem 8 % The Euler constant disp(' n E_n'); disp('-------------------------'); kMax = 100; step = 100; % Initialization total = 0; % Keep track of the summation for n=1:kMax*step % Calculate the sum total = total + 1/n; if rem(n,step) == 0 % Calculate E_n E_n = total - log(n); fprintf(' %5d %.8f\n', n, E_n) end end ----------------------------------------------------- % Sec 2.2 Problem 8 % Two sequences that approximate pi n = 0; a = 6/sqrt(3)*(-1)^n/3^n/(2*n+1); fprintf('a%d = %.7f\n',0, a) while abs(a - pi) > .000001 n = n + 1; a = a + 6/sqrt(3)*(-1)^n/3^n/(2*n+1); fprintf('a%d = %.7f\n',n, a) end n = 0; b = 16 * (-1)^n / 5^(2*n+1) / (2*n+1) - 4 * (-1)^n / 239^(2*n + 1); fprintf('\nb%d = %.7f\n',0, b) while abs(b - pi) > .000001 n = n + 1; b = b + 16 * (-1)^n / 5^(2*n+1) / (2*n+1) - 4 * (-1)^n / 239^(2*n + 1); fprintf('b%d = %.7f\n',n, b) end ----------------------------------------------------- % Section 3.1 Problem 8 % Draw ASCII figure % * % * * % * * % * * % * * % * * % * * % * * % * n= input('Input n, n>3: '); % Top half for r= 1:n % A general row has spaces, star, spaces, star, \n for c= 1:n-r fprintf(' '); end fprintf('*'); if (r > 1) for c = 1:2*(r-1)-1 fprintf(' '); end fprintf('*'); end fprintf('\n') end % Bottom half for r= n-1:-1:1 for c= 1:n-r fprintf(' '); end fprintf('*'); if (r > 1) for c = 1:2*(r-1)-1 fprintf(' '); end fprintf('*'); end fprintf('\n') end ----------------------------------------------------- %Section 3.1 Problem 10 %Print a sequence a(2), a(3), ..., a(n) where n is smallest integer %such that abs(a(n-1)-a(n))<=.01. % Initialization tol= 0.01; % Stopping tolerance n= 1; a_current= 0.5; % a(n), i.e., a(1) a_previous= 0; % a(n-1) while abs(a_previous - a_current) > tol % Update n= n+1; a_previous= a_current; % Calculate current term a(n), which is a summation a_current = 0; for j = 1:n^2 a_current = a_current + (n/(n^2+j^2)); end fprintf('%4d %10.6f\n', n, a_current) end ----------------------------------------------------- % Section 3.2 Problem 7 % Print t_1 to t_26 disp(' n t_n') disp('--------------------') n = 1; while n<=26 % Initialization current = 1; k = n; % Calculate from inside to outside. % current is the value of the kth squre root while k>1 || (n==1 && k>0) % Update... current = sqrt(1+k*current); k = k-1; end fprintf(' %2d %10.8f \n',n,current) % Update... n = n+1; end QUESTION 5 ------------------------------------------------------ % Assume variables nBig and nSm are given. % Outer loop to iterate from nBig to nSm (while-loop used here; % you can use for loop). while nBig >= nSm % Check nBig to see if it is prime d=2; % divisor % Iterate until first proper divisor is found while (mod(nBig,d) ~= 0); d = d + 1; end if (d == nBig) fprintf('%d is a prime\n',nBig); else fprintf('%d\n',nBig); end nBig = nBig – 1; end QUESTION 6 ----------------------------------------------------- % Example program for computing the mode of a sequence of integers. disp('Determine mode of a set of nonnegative integers.') disp('Use a negative number to quit.') previous = -1; % Previous number seen count = 1; % Count of current number mode = -1; % Mode seen so far modeCount = 0; % Count for mode so far number = input('Give me a number: '); while number >= 0 % Quit when negative number is encountered if number == previous % same run, so increment count count = count + 1; else % new run, so reset count count = 1; end if count > modeCount mode = number; modeCount = count; end previous = number; number = input('Another number not smaller than the previous: '); end if mode == -1 disp('Mode is undefined') else fprintf('Mode is %d which occurred %d times.\n', mode, modeCount) end QUESTION 7 ----------------------------------------------------- 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; %------------------------------------- % Question 7, 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 ); %------------------------------------- % Question 7, cont'd % Script to demonstrate functions hms2s and s2hms sec = hms2s( 7, 45, 13 ) [h, m, s ] = s2hms( 3682 ) QUESTION 8 ----------------------------------------------------- % 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 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)