Introduction to MATLAB Exercise 3

1. The formula \begin{equation} x_\pm = \frac{a\pm\sqrt{a^2-4}}{2} \end{equation} for the two roots of the special quadratic function \(x^2 - ax+1=0\) is well-known. When the absolute value of \(a\) is large, one of the root may not be accurate as expected. For example, if a=1e8, the smaller root computed from (a-sqrt(a^2-4))/2 is 7.450580596923828e-09 (using the default double precision in MATLAB), which is wrong in the first digit. In stead, we should use the equivalent formula \begin{equation} x_- = \frac{2}{a+\sqrt{a^2-4}} \end{equation} to avoid cancellation in the significant digits between \(a\) and \(\sqrt{a^2-4}\). Similarly, if \(a\) is negative, the two roots are computed from \begin{equation} x_+ = \frac{2}{a-\sqrt{a^2-4}},\qquad x_- = \frac{a-\sqrt{a^2-4}}{2}. \end{equation} Create a function quadroot by continuing the following block of code (save it as a matlab file "quadroot.m"). You can assume that the input argument \(a\) satisfies \(a^2>4\), and there is no need to discuss the case when \(a^2\leq 4\).

   function y=quadroot(a)
   %QUADROOT    Roots of the quadratic function x^2-a*x+1=0
   %   quadroot(a) returns the two roots of the quadratic function x^2-a*x+1=0
   %   such that both roots are accurate even if the absolute value of a is large
   

2. The number \(\sqrt{2}\) is irrational, but can be approximated through various iterations using elementary operations. Compute \(x_n\) and \(y_n\) which are defined by \begin{equation} x_{n} = \frac{x_{n-1}^2+2}{2x_{n-1}}, \qquad x_0=1 \end{equation} and \begin{equation} y_{n} = \frac{y_{n-1}+2}{y_{n-1}+1}, \qquad y_0=1. \end{equation} Find the number of iterations needed such that \(|x_n-\sqrt{2}|\) or \(|y_n-\sqrt{2}|\) is less than \(10^{-12}\). You can just use the same variable during the iteration, like x = (x^2+2)/2/x;, instead of saving all the history. Which iteration is faster (you can counter the number of steps, by adding the counter by one inside each iteration)?
3. During the lecture, the following function is used to determine whether the input is a leap year:

   function isleap(year)
   %% for the input year, determine whether this is a leap year
   
   if mod(year,4)~=0
       disp('This is not a leap year');
   elseif mod(year,100)~=0
       disp('This is a leap year');
   elseif mod(year,400)==0
       disp('This is a leap year');
   else
       disp('This is not a leap year');
   end
   

   This function can be written in alternative ways to make it clear, using the following figure (the shaded set indicates leap year).

   

   • We can use nested if statements as follow. Replace % insert the display here by either disp('This is not a leap year'); or disp('This is a leap year');.

     function isleap1(year)
     %% for the input year, determine whether this is a leap year
     % using nested "if" statements
     % save this file as "isleap1.m"
     if (mod(year,4)==0)
         if (mod(year,100)==0)
             if (mod(year,400)==0)
                 % insert the display here
             else
                 % insert the display here
             end
         else
             % insert the display here
         end
     else
         % insert the display here
     end
     

   • We can also use logical "and" & and/or logical "or" |. For example, there are two white regions representing non-leap years in the figure. The outer one is characterized by mod(year,4)~=0, and the inner one by (mod(year,100)==0 & mod(year,400)~=0). Therefore, the corresponding function becomes

     function isleap2(year)
     %% for the input year, determine whether this is a leap year
     % using logical "and" and "or"
     % save this file as "isleap2.m"
     
     if ((mod(year,4)~=0) | (mod(year,100)==0 & mod(year,400)~=0))
         disp('This is not a leap year');
     else
         disp('This is a leap year');
     end	
     

     Now write a function using some logical statements that first characterize the two shaded region in the figure. You can just modify the following code (add some statements inside the bracket following if:

     function isleap3(year)
     %% for the input year, determine whether this is a leap year
     % using logical "and" and "or"
     % save this file as "isleap3.m"
     
     if ()
         disp('This is a leap year');
     else
         disp('This is not a leap year');
     end	
     

   You can test your function by trying the following leap years (2000 and 2016) and non-leap years (1900 and 2015).
4. Approximate \(\pi\) by tiling a quarter disk. There are many ways to approximate the irrational number \(\pi\), for instance by counting the number of tiles (with unit size) of the following quarter disk with radius \(n\) (assumed to be an integer).
   • Write a function, for given integer \(n\) and row \(k\) (another integer not larger than \(n\)), such that number_of_tiles(n,k) returns the number of tiles in k-th row. In the following figure, \(n=8\) and you should get, for example, \(7\) for number_of_tiles(8,2) and \(5\) for number_of_tiles(8,6). How is the output related to \(n\) and \(k\)? (Hint: use floor the integer part of a decimal number, of course there are other related commands like ceil and round)
     
     

     

     function t = number_of_tiles(n,k)
     %% for the input integer radius n, output the number of tiles in k-th row
     % output a vector, if the input k is also a vector (we would like to call "number_of_tiles(n,1:n)" ),
     % so use vector operator if possible
     

     In fact, the function is so short, you can write it more conveniently as an anonymous function

     number_of_tiles = @(n,k)  % complete this line
     

   • As \(n\) increase, the area total number of tiles (the sum of number_of_tiles(n,1:n)) is about the area of the quarter circle \( \pi n^2/4\). Therefore, \(\pi\) can be approximated as

     format long;
     n = 100; [4*sum(number_of_tiles(n,1:n))/n^2; pi]
     

     Increase the value of \(n\), say \(n=1000\) and \(n=10000\), to see how the approximation approaches \(\pi\).
5. Collatz Conjecture. Consider the following operation on an arbitrary positive integer from \(n\) to \(f(n)\) \[ n \mapsto f(n)=\begin{cases} 3n+1, &\mbox{if } n \mbox{ is odd},\cr n/2, & \mbox{if } n \mbox{ is even}. \end{cases} \] That is, if the number is even, divide it by two; if the number is odd, triple it and add one. For example, by applying this operation successively, we get \[ 3\ \mapsto\ 10\ \mapsto\ 5\ \mapsto\ 16\ \mapsto\ 8\ \mapsto\ 4\ \mapsto\ 2\ \mapsto\ 1\ \mapsto\ 4\ \mapsto\ 2\ \mapsto\ 1\ \mapsto\ \cdots. \] It was conjectured that, starting with any positive integer, this operation eventually ends up with the periodic sequence \[ 1\ \mapsto\ 4\ \mapsto\ 2\ \mapsto\ 1\ \mapsto\ 4\ \mapsto\ 2\ \mapsto\ 1 \cdots. \] Write a function to determine the number of steps such that the number \(1\) first appears in the sequence.

   function counter = collatz(N)
   %% For input positive integer N, return of the number of steps, such that "1" first appears
   counter = 0; % the initial value
   n = N;
   while (n ~= 1)
       % complete this part, perform one step of the iteration on n, and increase the counter by one
   end
   

   You can test your function, for the following values: (1) the output of collatz(1) should be 0; (2) the output of collatz(2) should be 1; (3) the output of collatz(3) should be 7. Some small numbers may take many iterations to terminate. For example, the output for collatz(27) should be larger than 100.