\documentclass{article}

\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphics}
\usepackage{listings}
\usepackage{tcolorbox}
\usepackage{verbatim}

\DeclareMathOperator*{\argmin}{arg\,min}

\marginparwidth 0pt
\oddsidemargin  0pt
\evensidemargin  0pt
\marginparsep 0pt
\topmargin   0pt
\textwidth   6.5in
\textheight  8.5in
\setlength{\parindent}{0em}

\newenvironment{packed_itemize}{
\begin{itemize}
  \setlength{\itemsep}{0pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
}{\end{itemize}}

\newenvironment{packed_enumerate}{
\begin{enumerate}
  \setlength{\itemsep}{1pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
}{\end{enumerate}}

\definecolor{codeback}{RGB}{235,255,255}

\lstset{ %
language=Matlab,                % choose the language of the code
basicstyle=\footnotesize,       % the size of the fonts that are used for the code
numbers=left,                   % where to put the line-numbers
numberstyle=\footnotesize,      % the size of the fonts that are used for the line-numbers
stepnumber=1,                   % the step between two line-numbers. If it is 1 each line will be numbered
numbersep=10pt,                  % how far the line-numbers are from the code
backgroundcolor=\color{codeback},  % choose the background color. You must add \usepackage{color}
showspaces=false,               % show spaces adding particular underscores
showstringspaces=false,         % underline spaces within strings
showtabs=false,                 % show tabs within strings adding particular underscores
frame=single,           % adds a frame around the code
tabsize=2,          % sets default tabsize to 2 spaces
captionpos=b,           % sets the caption-position to bottom
breaklines=true,        % sets automatic line breaking
breakatwhitespace=false,    % sets if automatic breaks should only happen at whitespace
escapeinside={\%*}{*)}          % if you want to add a comment within your code
}

\begin{document}

\begin{center}
  \Huge{MATLAB Demonstration}\\
  \Large{MATH1070: Numerical Mathematical Analysis}
\end{center}
Location: http://people.sc.fsu.edu/$\sim$jburkardt/classes/math1070\_2019/matlab\_demo/matlab\_demo.pdf
\vskip 0.1in
\hrule
\vskip 0.1in

\begin{center}
  \includegraphics[scale=0.30]{two_flies.png}
  \vskip 0.1in
  {\it{Two flies land in the unit square.  On average, how far apart are they?}}
\end{center}

MATLAB lets you implement numerical algorithms, and interactively experiment
with them.   It has many numerical algorithms as built-in commands, 
and uses a linear-algebra notation that easily expresses vector and array operations.
It includes extensive graphics.  MATLAB ``toolboxes'' supply additional packages of
very specialized algorithms for symbolic mathematics, machine learning, signal analysis,
image processing, financial mathematics, and bioinformatics.
\vskip 0.1in
\begin{tcolorbox}[colback=red!5,title=MATLAB Demo]
{\it{
  How can we use a scientific programming language like MATLAB to:
  \begin{packed_itemize}
    \item{create a mathematical model of a real-world problem;}
    \item{devise an algorithm that can simulate the problem;}
    \item{write a computer program to implement the algorithm;}
    \item{experimentally search for solutions;}
    \item{graphically illustrate the properties of the solution;}
  \end{packed_itemize}
}}
\end{tcolorbox}
\vskip 0.1in
This is a demonstration of MATLAB, and so we will concentrate on solving the problem,
and hope that the MATLAB commands we use will make sense to you.
\vskip 0.1in
This particular fly-paper problem can be solved analytically, but we are interested in how
MATLAB can simulate the problem, estimate the answer, and produce convincing tables
and plots.  We are studying a random process, and trying to say something (the average)
about a huge number of cases.  We begin by simulating a single case.

\section{One simulation}

Being mathematicians, we assume the paper is a $1 \times 1$ square.  To simulate
the randomly chosen position of the first fly, $p_1=(x_1,y_1)$, we need to choose
uniform random values between 0 and 1.  Then we can place the second fly.  Once
we have their locations, we simply use the Pythagorean theorem to get the distance.
\begin{packed_enumerate}
  \item{First fly: {\tt{p1 = rand(2,1)}} }
  \item{Second fly: {\tt{p2 = rand(2,1)}} }
  \item{Distance: {\tt{d = norm ( p1 - p2 )}} }
  \item{Plot: 
  \begin{lstlisting}
    hold ( 'on' )
    plot ( p1(1), p1(2), 'r*' )
    plot ( p2(1), p2(2), 'b*' )
    plot ( [p1(1),p2(1)], [p1(2),p2(2)], 'k-' )
    axis ( [0,1,0,1] )
    hold ( 'off')
  \end{lstlisting}
  }
\end{packed_enumerate}

\begin{center}
%\begin{figure}
% \centering
  \includegraphics[scale=0.50]{simulate_one.png}
% \caption{{\bf{A few simulations}}}
%\end{figure}
\end{center}

\section{Many Simulations, Averaged}

The average distance between random flies can be estimated by simulating many examples.
In MATLAB, the {\it{for}} loop lets us repeat an operation many times.  We can
initialize a variable {\tt{average}} to zero, loop over {\tt{n}} simulations, adding
each distance to {\tt{average}}, and then dividing by {\tt{n}}:
\begin{lstlisting}
  average = 0.0;
  n = 10;
  for test = 1 : n
    p1 = rand(2,1);
    p2 = rand(2,1);
    d = vecnorm ( p1 - p2 );
    average = average + d;
  end
  average = average / n
\end{lstlisting}
How much variation do we see in our answers?  Let's try {\tt{n = 100}}
and compare answers.  What about {\tt{n = 1000}}?
\vskip 0.1in
{\bf{About how large a value of $n$ do we need to use before we are pretty confident
that the first two digits after the decimal point are correct?}}

\section{Using Vector Notation}

\begin{lstlisting}
  n = 10;
  p1 = rand(2,n);
  p2 = rand(2,n);
  d = vecnorm ( p1 - p2 );
  average = mean ( d )
\end{lstlisting}

Notice that, when moving to the vector version of the calculation, each of
the three ways of computing the distance {\tt{d}} had to be adjusted in some way.

The advantages of writing a calculation in the vector version are that
the code is somewhat more compact, it is easier to see that the computation
could be done in parallel, and in particular, MATLAB is able to compute many
vector operations very efficiently.

\section{Timing Comparison}

Let's compare the scalar and vector codes when {\tt{n}} is large.

First, let's rewrite our scalar calculation as an {\it{M-file}}, which we
will call {\it{scalar\_code.m}}, with the value of {\tt{n}} an input parameter.

\begin{lstlisting}
function average = scalar_code ( n )
  average = 0.0;
  for test = 1 : n
    p1 = rand(2,1);
    p2 = rand(2,1);
    d = vecnorm ( p1 - p2 );
    average = average + d;
  end
  average = average / n
  return
end
\end{lstlisting}

Doing this means we have created our own MATLAB command, which can be executed
with a single line:
\begin{lstlisting}
  average = scalar_code ( n )
\end{lstlisting}
We can similarly create {\it{vector\_code.m}}.  Now we could compare the speed
of the two methods by simply waiting for each to complete, or using a stopwatch,
but instead, we will use the MATLAB commands {\tt{tic}} and {\tt{toc}}:

\begin{lstlisting}
  for logn = 1 : 6
    n = 10^logn;
    tic;
    scalar_code ( n );
    t1 = toc;
    tic;
    vector_code ( n );
    t2 = toc;
    fprintf ( 1, '%d  %g  %g  %g\n', n, t1, t2, t1/t2 );
  end
\end{lstlisting}

\section{A convergence table}

The average is an average over an infinite number of cases.  We are estimating
that average using {\tt{n}} cases.  How do we convince ourselves that the estimate
is a good approximation?  One way is to repeat the calculation for larger values
of {\tt{n}} and make a table.  If the estimates seem to converge, then this 
suggests that our many choices of {\tt{n}} seem to point towards a common answer
{\it{although we can't be sure this answer is right!}}.

\begin{lstlisting}
  average_old = 0.0;
  for logn = 1 : 6
    n = 10^logn;
    average = vector_code ( n );
    fprintf ( 1, '%7d  %14.6f  %14.6f\n', n, average, average - average_old );
    average_old = average;
  end
\end{lstlisting}
Our results show that even a value of {\tt{n=1,000,000}} only gives us 
confidence in the first two decimals of our estimate.
\begin{lstlisting}
>> converge
     10        0.476492        0.476492
    100        0.535507        0.059016
   1000        0.518386       -0.017121
  10000        0.523612        0.005226
 100000        0.521069       -0.002543
1000000        0.521247        0.000178
\end{lstlisting}

\section{A histogram of the results}

In the problem we are thinking about, possible distances $d$ range between 0 and 
$\sqrt(2)$.  Some distances are more likely than others, as described by a 
probability density function.  If we choose {\tt{n}} large enough, we can make
a histogram of our results that will suggest what the true PDF looks like.

The basic command is simply
\begin{lstlisting}
  histogram ( d )
\end{lstlisting}
but I want to specify the number of bins to be 20, and I want the vertical
$y$ axis to be scaled like a PDF, rather than a frequency count.
\begin{lstlisting}
  n = 10000;
  p1 = rand(2,n);
  p2 = rand(2,n);
  d = vecnorm ( p1 - p2 );
  bins = 20;
  histogram ( d, bins, 'Normalization', 'pdf' )
\end{lstlisting}

\begin{center}
  \includegraphics[scale=0.50]{distance_histogram.png}
\end{center}

\section{Histogram versus PDF}

For this particular problem, there happens to be an exact formula for the PDF.
We can plot it together with the histogram, to see how we did.  The formula
for the exact PDF is a little messy:
\begin{displaymath}
  \text{pdf}(d) =
\begin{cases}
  2 d ( d^2 - 4 d + \pi ) & \text{if } d <= 1.0 \\
  2 d (4 \sqrt { d^2 - 1.0 } - ( d^2 + 2 - \pi )
      - 4 \arctan ( \sqrt { d^2 - 1 } ) ) & \text{if } 1.0 < d
\end{cases}
\end{displaymath}
\begin{center}
  \includegraphics[scale=0.50]{distance_histogram_compare.png}
\end{center}

\section{Challenge:}

Consider a new problem:
{\it{Two flies land inside the unit {\bf{disk}}.  On average, how far apart are they?}}
\vskip 0.1in
To pick a random point $(r,\theta)$ in the unit circle:
\begin{packed_enumerate}
  \item{pick a random number $r_1$, and set $\theta = 2 \pi r_1$; }
  \item{pick a random number $r_2$, and set $r = \sqrt{r_2}$;}
\end{packed_enumerate}
\vskip 0.1in
If you are surprised by the computation of $r$, plot 1,000 random points this way,
and than compare with 1,000 random points computed by $r=r_2$.
\vskip 0.1in
Generate a lot of pairs (the right way!), save the distances in a vector $d$, 
and then histogram the results.  Compare your histogram to a plot of the exact 
PDF for the distance $0 \leq d \leq 2$:
\begin{displaymath}
  pdf(d) = \frac{d}{\pi} (4\arccos(\frac{d}{2}) - d \sqrt{4-d^2})
\end{displaymath}

\end{document}
