\documentclass[12pt]{article}
\usepackage{fullpage,mathpazo,amsmath,amssymb,amsfonts,nicefrac,microtype,graphicx}

\usepackage{hyperref,color,textcomp}
\usepackage[skip=-2pt,font=footnotesize]{caption}

\definecolor{webgreen}{rgb}{0,.35,0}
\definecolor{webbrown}{rgb}{.6,0,0}
\definecolor{RoyalBlue}{rgb}{0,0,0.9}
\hypersetup{
   colorlinks=true, linktocpage=true, pdfstartpage=3, pdfstartview=FitV,
   breaklinks=true, pdfpagemode=UseNone, pageanchor=true, pdfpagemode=UseOutlines,
   plainpages=false, bookmarksnumbered, bookmarksopen=true, bookmarksopenlevel=1,
   hypertexnames=true, pdfhighlight=/O,
   urlcolor=webbrown, linkcolor=RoyalBlue, citecolor=webgreen,
   pdfauthor={Chris H. Rycroft, Yue Sun},
   pdfsubject={UW--Madison Math/CS 714 (Fall 2025)},
   pdfkeywords={},
   pdfcreator={pdfLaTeX},
   pdfproducer={LaTeX with hyperref}
}
\hypersetup{pdftitle={Math/CS 714: Homework 2}}

\newcommand{\N}{\mathbb{N}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\p}{\partial}
\renewcommand{\vec}[1]{\mathbf{#1}}
\newcommand{\vx}{\vec{x}}
\newcommand{\vy}{\vec{y}}
\newcommand{\Lag}{\mathcal{L}}
\newcommand{\sep}{\,:\,}
\setlength{\unitlength}{0.5mm}
\renewcommand{\labelitemii}{$\diamond$}

\graphicspath{{./}{gr/}}

\begin{document}
\section*{Math/CS 714: Assignment 2}
\begin{enumerate}
  \item \textbf{Iterative methods (2 points).}
        \begin{enumerate}
          \item Consider the $2\times 2$ matrix
                \begin{equation}
                  A = \left( \begin{array}{cc} 1 & \rho \\ - \rho & 1 \end{array}\right).
                \end{equation}
                Under what conditions will the Jacobi and Gauss--Seidel methods converge?
          \item Consider the $n \times n$ matrix
                \begin{equation}
                  \footnotesize
                  C = \left(
                  \begin{array}{cccccc}
                    3  & -1 &    &        &    &    \\
                    -1 & 3  & -1 &        &    &    \\
                       & -1 & 3  & -1     &    &    \\
                       &    &    & \ddots &    &    \\
                       &    &    & -1     & 3  & -1 \\
                       &    &    &        & -1 & 3
                  \end{array}
                  \right).
                \end{equation}
                Starting from $u_0\in\R^n$, for which values of $\omega \in \R$ does the iteration
                \begin{equation}
                  u_{k+1} = u_k + \omega (b - C u_k)
                \end{equation}
                for $k=0,1,2,\ldots$ converge to a solution of $Cu=b$? What iterative
                method from the lectures does this most closely resemble? How is it
                different?
        \end{enumerate}

  \item \textbf{Triangular domain revisited with conjugate gradient (8 points).}
        Question 4 of the first homework assignment looked at solving the Poisson equation
        \begin{equation}
          \nabla^2 u = f \label{eq:poiss}
        \end{equation}
        on the triangular domain $T$ with vertices at $(0,0)$, $(1,0)$, and
        $(s,\tfrac12)$ where $s=\tfrac{\sqrt{3}}{2}$. Dirichlet boundary conditions $u=0$ are
        applied on the boundary $\p T$. For $n\in \N$ and $h=1/n$, the domain is
        discretized with points
        \begin{equation}
          \vx_{i,j} = (h(i+\tfrac12 j), hsj))
        \end{equation}
        for $0\le i\le n$, $0\le j \le n-i$. The points where $i=0$, $j=0$, or
        $i+j=n$ are boundary points, and all others are interior. Let $u_{i,j}$
        be the numerical approximation for $u(\vx_{i,j})$. The Laplacian
        is discretized using
        \begin{equation}
          \nabla^2_6u_{i,j} = \frac{2(-6 u_{i,j} + u_{i+1,j} + u_{i,j-1} + u_{i-1,j+1} +u_{i-1,j} + u_{i,j+1} + u_{i+1,j-1})}{3h^2}.
        \end{equation}
        As in Homework 1, consider solving Eq.~\eqref{eq:poiss} using the $f$ that creates the exact
        solution
        \begin{equation}
          u^{\text{ex}}(x,y) = \left( (2y-\sqrt{3})^2- 3(2x-1)^2 \right) \sin y.
        \end{equation}
        The problem can be expressed as a linear system $Au=f$ where $u\in \R^N$ is
        the solution vector at the $N=(n-1)(n-2)/2$ interior grid points. Define
        an appropriate 2-norm as
        \begin{equation}
          \|r\|_2 = \sqrt{\frac{s}{2n^2} \sum_{j=1}^{n-1} \sum_{i=1}^{n-j-1} r_{i,j}^2}
        \end{equation}
        for a vector $r$ describing a field on the grid.
        \begin{enumerate}
          \item For a range of grid sizes from $n=10$ up to at least $n=160$
                measure the wall-clock time $T_n$ to solve the system using your
                code.\footnote{If you had difficulty with completing this question,
                  you can alternatively use the code \texttt{poisson\_\,tri.py} from
                  the homework 1 solutions.} By making a log--log plot of $T_n$ versus
                $n$, fit the timing data to
                \begin{equation}
                  T_n = a n^b, \label{eq:time_model}
                \end{equation}
                and comment on the exponent $b$.
          \item Write a code to implement the conjugate gradient (CG) algorithm to
                solve $Au=f$. Your code should not explicitly build $A$ as a dense
                matrix. It should either represent $A$ as a sparse matrix, or contain
                a function that can directly compute $Aq$ for a given vector $q$. The
                CG algorithm should terminate when $\|r\|_2 < 10^{-10}$, where $r$ is
                the residual vector. For the case when $n=40$ test that your program
                gives the same results as the original version.
          \item For $n=10,20,40,80,160$, calculate the number of iterations $k$
                required in order to reach the termination criterion.
          \item Measure the wall-clock time of the CG algorithm to solve the
                linear system for a range of grid sizes from $n=10$ up to at least
                $n=160$. Fit the data to Eq.~\eqref{eq:time_model} and comment on
                how the exponent compares to your result from part (a).
          \item Consider the block Jacobi preconditioner $M$ with block sizes of
                $w=\lfloor \sqrt{N} \rfloor$. In general, $w$ will not evenly divide
                $N$, and therefore let the final block of $M$ be smaller. Implement
                the preconditioned CG algorithm, and repeat part (c) to measure the
                number of iterations required to reach the termination criterion.
          \item Repeat part (d), fitting the timing data to
                Eq.~\eqref{eq:time_model}
          \item \textbf{Optional.} Would this question be possible if
                $\nabla_3^2$ was used instead of $\nabla_6^2$?
          \item \textbf{Optional.} With the block Jacobi preconditioner, the
                ordering of the points within the linear system will affect the
                performance. In particular, having a large number of matrix entries
                within the Jacobi blocks will likely make the preconditioner more
                effective. Consider approaches to reorder the points to achieve
                better performance.
        \end{enumerate}

  \newpage
  \item \textbf{ODE integration methods (7 points).}
        \begin{enumerate}
          \item Consider solving the differential equation $y'=f(t,y)$ at
                timepoints $t_k$ with corresponding numerical solutions $y_k$. The
                multi-step Nystr\"om numerical method is based upon the integral
                relation
                \begin{equation}
                  y(t_{k+1}) = y(t_{k-1}) + \int_{t_{k-1}}^{t_{k+1}} f(t,y) dt.
                \end{equation}
                Derive an implicit multi-step numerical method by approximating the
                integrand $f(t,y)$ with the polynomial interpolant using the function
                values at $t_{k-2}$, $t_{k-1}$, $t_k$, and $t_{k+1}$. Your method
                should have the form
                \begin{equation}
                  y_{k+1} = y_{k-1} + h(\alpha f_{k-2}+\beta f_{k-1} + \gamma f_k+\eta
                  f_{k+1}) \label{eq:nys}
                \end{equation}
                where $\alpha$, $\beta$, $\gamma$, and $\eta$ are constants to be
                determined, $h$ is the timestep interval size, and $f_l=f(t_l,y_l)$ for
                an arbitrary $l$.
          \item Find the exact solution to the second-order differential equation
                \begin{equation}
                  y''(t)+2y'(t) + 26y(t)=0. \label{eq:ode2}
                \end{equation}
                subject to the initial conditions $y(0)=1$, $y'(0)=0$.
          \item Write Eq.~\eqref{eq:ode2} as a coupled system of two first-order
                differential equations for $\vy=(y,v)=(y,y')$. Solve the system
                over the interval $0\le t\le 3$ with a timestep of $h=0.02$ using your
                multi-step method from part (a).\footnote{Even though Eq.~\ref{eq:nys}
                  is an implicit method and would be hard to solve in general, the
                  linearity of the system in Eq.~\ref{eq:ode2} means that for this case,
                  the update equation for $\vy_{k+1}$ can be derived analytically.}

                Before Eq.~\eqref{eq:nys} can be applied, $\vy_1$ and $\vy_2$
                must be calculated accurately. Use one of the following two approaches:
                \begin{enumerate}
                  \item set them based on the exact solution from part (b),
                  \item calculate them using the classical fourth-order Runge--Kutta
                        method.
                \end{enumerate}
                Plot the exact and numerical solutions over the range $0\le t\le 3$.

                Make a log--log plot of the absolute error between the numerical and
                exact values of $y$ at $t=3$ as a function of $h$, over the range from
                $h=10^{-3}$ to $h=10^{-1}$. Show that your method is fourth-order
                accurate.
          \item Suppose that instead of setting $\vy_1$ and $\vy_2$ accurately,
                you instead make use of forward Euler steps. Create a log--log plot
                of the absolute error of $y$ at $t=3$ as a function of $h$. Determine
                the order of accuracy, and discuss why this is the case.
        \end{enumerate}

  \item \textbf{Linear difference equation (3 points).}
        \begin{enumerate}
          \item Find the general solution of the linear difference equation
                \begin{equation}
                  U_{n+3} + 2U_{n+2} - 4U_{n+1} - 8U_n =0. \label{eq:dfnce}
                \end{equation}
          \item Determine the particular solution with initial data $U_0=4$, $U_1=-2$, and $U_2=8$.
          \item Consider the iteration
                \begin{equation}
                  \left(\begin{array}{c} U_{n+1} \\ U_{n+2} \\ U_{n+3} \end{array} \right)
                  =
                  \left( \begin{array}{ccc} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 8 & 4 & -2 \end{array} \right)
                  \left(\begin{array}{c} U_{n} \\ U_{n+1} \\ U_{n+2} \end{array} \right),
                \end{equation}
                which we can define as $\mathbf{U}_{n+1} = A \mathbf{U}_n$. The matrix
                $A$ is called the companion matrix for the difference equation in
                Eq.~\eqref{eq:dfnce}. A general solution of the difference
                equation is given by $\mathbf{U}_{n} = A^n \mathbf{U}_0$. If
                $A=RJR^{-1}$ is the Jordan canonical form for $A$, then $A^n = RJ^n
                  R^{-1}$. Determine the eigenvalues and Jordan canonical form for this
                matrix and show how this is related to the general solution found in
                part (a).
        \end{enumerate}
\end{enumerate}
\end{document}
