\documentclass[]{article}
\usepackage{lmodern}
\usepackage{amssymb,amsmath}
\usepackage{ifxetex,ifluatex}
\usepackage{fixltx2e} % provides \textsubscript
\ifnum 0\ifxetex 1\fi\ifluatex 1\fi=0 % if pdftex
  \usepackage[T1]{fontenc}
  \usepackage[utf8]{inputenc}
\else % if luatex or xelatex
  \ifxetex
    \usepackage{mathspec}
    \makeatletter % undo the wrong changes made by mathspec
    \let\RequirePackage\original@RequirePackage
    \let\usepackage\RequirePackage
    \makeatother
  \else
    \usepackage{fontspec}
  \fi
  \defaultfontfeatures{Ligatures=TeX,Scale=MatchLowercase}
\fi
% use upquote if available, for straight quotes in verbatim environments
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
% use microtype if available
\IfFileExists{microtype.sty}{%
\usepackage{microtype}
\UseMicrotypeSet[protrusion]{basicmath} % disable protrusion for tt fonts
}{}
\usepackage{hyperref}
\hypersetup{unicode=true,
            pdftitle={Homework 4},
            pdfborder={0 0 0},
            breaklinks=true}
\urlstyle{same}  % don't use monospace font for urls
\IfFileExists{parskip.sty}{%
\usepackage{parskip}
}{% else
\setlength{\parindent}{0pt}
\setlength{\parskip}{6pt plus 2pt minus 1pt}
}
\setlength{\emergencystretch}{3em}  % prevent overfull lines
\providecommand{\tightlist}{%
  \setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}
\setcounter{secnumdepth}{0}
% Redefines (sub)paragraphs to behave more like sections
\ifx\paragraph\undefined\else
\let\oldparagraph\paragraph
\renewcommand{\paragraph}[1]{\oldparagraph{#1}\mbox{}}
\fi
\ifx\subparagraph\undefined\else
\let\oldsubparagraph\subparagraph
\renewcommand{\subparagraph}[1]{\oldsubparagraph{#1}\mbox{}}
\fi

\title{Homework 4}

\input{preamble.tex}
\title{Homework}



\begin{document}
\maketitle

Please answer the following questions in complete sentences in a typed
manuscript and submit the solution on Gradescope by April 5th around 5am
like our usual deadline.

\hypertarget{problem-0-homework-checklist}{%
\subsection{Problem 0: Homework
checklist}\label{problem-0-homework-checklist}}

\begin{itemize}
\item
  Please identify anyone, whether or not they are in the class, with
  whom you discussed your homework. This problem is worth 1 point, but
  on a multiplicative scale.
\item
  Make sure you have included your source-code and prepared your
  solution according to the most recent Piazza note on homework
  submissions.
\end{itemize}

\hypertarget{problem-1-gautschi-exercise-3.29-more}{%
\subsection{Problem 1: Gautschi Exercise 3.29 +
more}\label{problem-1-gautschi-exercise-3.29-more}}

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  Derive the 2-point Gauss-Hermite quadrature rule.
\item
  Example application (from
  \url{http://ice.uchicago.edu/2012_presentations/Faculty/Judd/Quadrature_ICE11.pdf}).
  An investor holds one bond that will be worth 1 in the future and
  equity whose value is \(Z\) where \(\log Z \sim N(\mu,\sigma^2)\). (So
  this means that the log of the value of the expected utility is a
  normally distributed random variable.) The expected utility is the
  random number \(U = f(1 + Z)\). where \(f\) is a utility function,
  we'll use \(f(x) = x^{1+\gamma} / (1+\gamma)\), where \(\gamma < 0\).
  (This is a concave utility function because having more money doesn't
  give you all that much more utility.) We'll use \(\gamma = -0.5\).
  Suppose also that \(\mu = 0.15\) and \(\sigma = 0.25\). We want to
  find the expected utility to the investor! This involves evaluating
  the integral
  \[ E[U] = \int_{-\infty}^\infty f(1+e^x) e^{-(x-\mu)^2/(2 \sigma^2)} \, dx. \]
  Write a compute program to use Gauss-Hermite quadrature to approximate
  the value of this integral. You need to justify the number of points
  you use for the approximation.
\end{enumerate}

\hypertarget{problem-2-multivariate-quadrature}{%
\subsection{Problem 2: Multivariate
quadrature}\label{problem-2-multivariate-quadrature}}

In this problem, we will investigate multivariate quadrature for
integrals such as \[ I = \int_{-1}^1 \int_{-1}^1 f(x,y) \, dx dy \]
using \emph{tensor product} rules. Let \(t_i\) and \(w_i\) be the nodes
and weights of an \(n\)-point 1-dimensional Gauss-Legendre rule. Then
the multidimensional quadrature rule is:
\[ I \approx \sum_{i=1}^n \sum_{j=1}^n f(t_i, t_j) w_i w_j. \] We can
derive this as follows:
\[ I \approx \int_{-1}^1 \sum_{i=1}^n f(t_i,y) w_i dy \approx \sum_{i=1}^n \sum_{j=1}^n f(t_i, t_j) w_i w_j. \]
Of course, this makes it clear we don't have to use the same number of
points to integrate in \(x\) and \(y\)! So in general, let
\(t_i^{(x)}, w_i^{(x)}\) be an \(N_x\)-point Gauss-Legendre quadrature
rule for the \(x\)-variable and \(t_i^{(y)}, w_i^{(y)}\) be an
\(N_y\)-point Gaussian quadrature for the \(y\) variable.

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  Implement a computer code to perform two-dimensional Gauss-Legendre
  quadrature using this type of construction. Your code should allow a
  user to input \(N_x\) and \(N_y\) to determine the number of points in
  each variable.
\item
  Use your code to estimate the integrals using \(10\) points in each
  dimension.

  \begin{itemize}
  \tightlist
  \item
    \(f(x,y) = x^2 + y^2\)
  \item
    \(f(x,y) = x^2 y^2\)
  \item
    \(f(x,y) = exp(x^2 + y^2)\)
  \item
    \(f(x,y) = (1-x^2) + 100*(y - x^2)^2\)
  \end{itemize}
\item
  \emph{This is an open ended question that requires you to
  investigate.} You can also find the answer in many textbooks, but if
  you do so, make sure you document your sources and demonstrate the
  effect that is claimed. We saw in class that an \(n\)-point Gaussian
  quadrature rule exactly integrated polynomials of up to degree
  \(2n-1\). In this problem, I want you to generate a conjecture about
  the degree of exactness of multidimensional Gaussian quadrature. You
  should use your code from part 1, along with carefully constructed
  examples, to support a statement such as:

\begin{verbatim}
 My evidence suggests that 2d Gauss-Legendre quadrature will exactly 
 integrate two-dimensional functions $f(x,y)$ when ...
\end{verbatim}

  Here are some helpful ideas that may play a role.

  \begin{itemize}
  \item
    The total degree of a multidimensional polynomial is the maximum of
    the sum of degrees of each term. So \(f(x,y) = x^2 y^2\) has total
    degree \(4\).
  \item
    Another type of degree is the largest degree in each variable, so
    \(f(x,y) = x^2 y^2\) involves polynomials of degree \(2\) only.
  \end{itemize}
\end{enumerate}

\hypertarget{problem-3-monte-carlo-gaussian-legendre-and-clenshaw-curtis-quadrature}{%
\subsection{Problem 3: Monte Carlo, Gaussian-Legendre, and
Clenshaw-Curtis
Quadrature}\label{problem-3-monte-carlo-gaussian-legendre-and-clenshaw-curtis-quadrature}}

In this problem, we'll study how different quadrature methods converge
on a variety of problems. For a technical paper on this idea, see
\href{https://people.maths.ox.ac.uk/trefethen/publication/PDF/2008_127.pdf}{Trefethen,
Is Gaussian Quadrature better than Clenshaw-Curtis? SIAM Review, 2008}.
In this problem, we'll be studying and reproducing Fig. 2 from that
paper.

The functions are all defined on the region \([-1,1]\) and are:

\begin{itemize}
\tightlist
\item
  \(f_1(x) = x^{20}\)
\item
  \(f_2(x) = e^{-x^2}\)
\item
  \(f_3(x) = e^{-1/(x^2)}\)
\item
  \(f_4(x) = e^x\)
\item
  \(f_5(x) = 1/(1+16x^2)\)
\item
  \(f_6(x) = |x|^3\)
\end{itemize}

The quadrature methods are:

\begin{itemize}
\item
  \textbf{Monte Carlo quadrature.} Monte Carlo quadrature is a
  randomized method. We simply guess \(n\) points between \([-1,1]\)
  uniformly at random and then take the average of all the function
  values. For instance, the following Julia code evaluates a Monte Carlo
  approximation

\begin{verbatim}
  function montecarlo(f::Function, n::Int, a::Float64, b::Float64)
      assert(a < b)
      xi = (b-a)*rand(n)+a
      fi = mean(f(xi))
  end
\end{verbatim}
\item
  \textbf{Clenshaw-Curtis quadrature.} This quadrature uses Chebyshev
  points of the second kind to build an interpolatory quadrature formula
  instead of uniformly spaced points (as is common in Newton-Cotes
  quadrature). It just so happens that there is an incredibly elegant
  method to compute the weights associated with this quadrature based on
  the Fast-Fourier transform. See Trefethen's paper above for a 6-line
  Matlab code that implements Clenshaw-Curtis quadrature.
\item
  \textbf{Gauss-Legendre quadrature.} In Gauss-Legendre quadrature, we
  pick the quadrature nodes and weights together. This gives even more
  accuracy. To find these nodes and weights, we must evaluate the
  eigenvalues and one component of each eigenvector of the Jacobi matrix
  associated with the Legendre orthogonal polynomials. The Jacobi matrix
  for these polynomials is easy:
  \[ \mJ = \bmat{ 0 & \sqrt{1/(4-1/(1^2))} & 0 & 0 & 0 & \ldots \\
                \sqrt{1/(4-1/(1^2))} & 0 & \sqrt{1/(4-1/(2^2))} & 0 & 0 & \ldots \\
                0 & \sqrt{1/(4-1/(2^2))} & 0 & \sqrt{1/(4-1/(3^2))} & 0 & \ldots \\
                \vdots & \vdots & \ddots & \ddots & \ddots & \ldots \\
                \vdots & \vdots & \vdots & \vdots & \vdots & \ldots \\ } \]
  The size of the matrix should be \(n \times n\) if you want an
  \(n\)-point formula. To get the nodes, we just look at the eigenvalues
  of the matrix. To get the weights, we need to get the first component
  of each eigenvector, and square it. (Hint: see Trefethen's paper for a
  simple Matlab code.)
\end{itemize}

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  By whatever method you want, determine the exact values of these 6
  integrals. (Hint: Wolfram Alpha is okay!)
\item
  Write a program to create the Jacobi matrix for Gauss-Legendre
  quadrature and show the eigenvalues of the matrix for \(n=11\).
\item
  Implement a computer program for Clenshaw-Curtis quadrature. (Hint:
  you can copy Trefethen's routines, with attribution) but you must
  explain the steps. Julia implementations are advisable too.
\item
  Implement a computer program for Gauss-Legendre quadrature.
\item
  Note that the Monte Carlo method is a randomized algorithm, so the
  result will change if you do it again. Each run is called a trial or
  sample. Use a computer implementation of Monte Carlo integration to
  how much the values in a Monte Carlo integration can vary between one
  trial and the next. Compute the variance for \(n=100\) and \(1000\)
  samples for each of the above functions.
\item
  Prepare 6 figures like Trefethen had in his paper for the 6 functions.
  Except use Monte Carlo integration instead of Newton-Cotes.
\item
  Empirically determine how computing the Gaussian Quadrature nodes and
  weights scales as a function of \(n\), the number of points, in terms
  of CPU time. (Hint, consider \(n\) between \(100\) and \(1000\).) You
  should be able to justify your final scaling constant in terms of the
  runtime of a known algorithm.
\item
  (Optional) There are some recent developments in \emph{fast
  quadrature} methods that produce the nodes and weights much more
  quickly. In Julia, these are implemented in the
  \texttt{FastGaussQuadrature} package and the \texttt{gausslegendre}
  function and in Matlab, the function \texttt{legpts.m} from Chebfun
  implements them (an old simple version is here:
  \url{http://www.mathworks.com/matlabcentral/fileexchange/23972-chebfun-v4-old-version--please-download-current-version-instead/content/chebfun/legpts.m}
  Determine the simplest equation that describes the empirical scaling
  for these routines to find the quadrature points and weights. (Hint,
  consider \(n\) between \(1000\) and \(100000\).)
\end{enumerate}

\hypertarget{problem-4-gautschi-exercise-4.1}{%
\subsection{Problem 4: Gautschi Exercise
4.1}\label{problem-4-gautschi-exercise-4.1}}

The following sequences converge to \(0\) as \(n\to\infty\):

\begin{itemize}
\tightlist
\item
  \(v_n = n^{-10}\)
\item
  \(w_n = 10^{-n}\)
\item
  \(x_n = 10^{-n^2}\)
\item
  \(y_n = n^{10} 3^{-n}\)
\item
  \(z_n = 10^{-3 \cdot 2^n}\)
\end{itemize}

Indicate the type of convergence for each sequence in terms of

\begin{itemize}
\tightlist
\item
  Sublinear
\item
  Linear
\item
  Superlinear
\item
  Quadratic
\item
  Cubic
\item
  None of the above
\end{itemize}

\hypertarget{problem-5-implementing-a-core-routine}{%
\subsection{Problem 5: Implementing a core
routine}\label{problem-5-implementing-a-core-routine}}

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\item
  Using any method we've seen to solve a scalar nonlinear equation
  (bisection, false position, secant), develop a routine to compute
  \(\sqrt{x}\) using only addition, subtraction, multiplication, and
  division (and basic control structures) to numerical precision. (Use
  double-precision.)
\item
  Compare the results of your method to the Matlab/Julia/Python function
  \texttt{sqrt}. Comment on any differences that surprise you.
\end{enumerate}

\end{document}
