\documentclass[12pt]{article}

\usepackage{graphicx}			% Use this package to include images
\usepackage{caption}            % To include captions of figures inside frames
\usepackage{amsmath, amssymb, amsthm, gensymb, cases} % A library of many standard math expressions
\usepackage[margin=1in]{geometry}% Sets 1in margins. 
\usepackage{fancyhdr}			% Creates headers and footers
\usepackage{enumerate}          % These two packages give custom labels to a list
\usepackage[shortlabels]{enumitem}
\usepackage{matlab-prettifier}  % Import and format MATLAB code/files
\usepackage{hyperref}           % To include urls
\usepackage{xurl}               % To force long urls to respect margins
\usepackage{framed}             % To make frames/boxes


% Creates the header and footer. You can adjust the look and feel of these here.
\pagestyle{fancy}
\fancyhf{}                      % Clear all fields
\fancyhead[L]{CE 2989 \\Numerical Methods in Civil \& Env. Engineering \\ Prof. Taha} % Course number, course name, prof name
\fancyhead[R]{Final Exam \\ Joshua Oliver Abigania \\Due: May 01, 2025} % Due date
\fancyfoot[C]{\thepage} % Page number
\renewcommand{\headrulewidth}{0.2pt} % Creates a horizontal line underneath the header
\setlength{\headheight}{42pt} % Sets enough space for the header

% Lists:
\setlist[enumerate]{font=\bfseries}     % Sets all enumerate to be bold
\setlist[enumerate, 2]{wide, labelwidth=!, labelindent=0in}     % Set all 2nd level enumerate to have these properties

% new commands:
\newcommand{\units}[1]{\left[\textrm{#1}\right]}     % easier way to write bracketed units in math mode
\newcommand{\p}[1]{\left(#1\right)}                 % easier way to write inside parentheses in math mode
\newcommand{\bracket}[1]{\left[#1\right]}                 % easier way to write inside brackets in math mode
\newcommand{\abs}[1]{\left|#1\right|}                 % easier way to write inside absolute value bars in math mode
\newcommand{\prob}[1]{P\p{#1}}                      % easier way to write the probability of smth --> P(smth)
\newcommand{\at}[2]{\bigg{|}_{#1}^{#2}}                  % easier way to write function evaluations at a point
\newcommand{\norm}[1]{\left|\left|#1\right|\right|}
\newcommand{\bMATtwo}[4]{\begin{bmatrix}
#1 & #2 \\
#3 & #4 \\
\end{bmatrix}}          % easier way to write a bracketed 2x2 matrix
\newcommand{\vMATtwo}[4]{\begin{vmatrix}
#1 & #2 \\
#3 & #4 \\
\end{vmatrix}}          % easier way to write a 2x2 matrix inside vertical lines
\newcommand{\bMATthree}[9]{\begin{bmatrix}
#1 & #2 & #3 \\
#4 & #5 & #6 \\
#7 & #8 & #9 \\
\end{bmatrix}}          % easier way to write a bracketed 3x3 matrix
\newcommand{\vMATthree}[9]{\begin{vmatrix}
#1 & #2 & #3 \\
#4 & #5 & #6 \\
#7 & #8 & #9 \\
\end{vmatrix}}          % easier way to write a 3x3 matrix inside vertical lines
\newcommand{\ro}[1]{\xrightarrow{#1}}   % easier way to write row operations
\newcommand{\var}[1]{\textrm{Var}\bracket{#1}} % easier way to write variance
\newcommand{\cov}[2]{\textrm{Cov}\bracket{#1,#2}} % easier way to write covariance
\newcommand{\expValue}[1]{\textrm{E}\bracket{#1}} % easier way to write expected value



\begin{document} % The writing for your homework should all come after this. 

\begin{center}
    \textbf{Take-Home Final Exam} % Full assignment name
\end{center}

% Enumerate starts a list of problems so you can put each homework problem after each item. 
\begin{enumerate}[wide, start=1, label={\bfseries Problem \arabic*:}, labelwidth=!, labelindent=0in] % You can change "Problem" to be whatever label you like. 

% Problem 1 begins here
    \item \hfill
    \begin{framed}
        \noindent Given system of nonlinear ODEs $\dot{x} = f(x,t)$:
        \begin{gather*}
            \begin{cases}
                \dot{x}_1(t) = -2x_1(t) - x_4^2(t) + 5x_2^2(t) + \cos(x_1(t)) \\
                \dot{x}_2(t) = -x_2(t) + x_3(t) \\
                \dot{x}_3(t) = -4x_3^2(t) + \sin(x_2(t)) \\
                \dot{x}_4(t) = -x_4(t)
            \end{cases}
        \end{gather*}
        \noindent for $x \in \mathcal{X} = \{0\le x_i\le2,i=1,2,3,4\}$
    \end{framed}
    \begin{enumerate}
% 1a
        \item \hfill
        \begin{framed}
            \noindent Find $K$ such that $\norm{f(x^{(1)}) - f(x^{(2)})} \le K\norm{x^{(1)} - x^{(2)}}$. \\ \\
            \noindent We have:
            \begin{gather*}
                x^{(1)} = \begin{bmatrix}
                    x_1^{(1)} \\ x_2^{(1)} \\ x_3^{(1)} \\ x_4^{(1)}
                \end{bmatrix} \textrm{ and } x^{(2)} = \begin{bmatrix}
                    x_1^{(2)} \\ x_2^{(2)} \\ x_3^{(2)} \\ x_4^{(2)}
                \end{bmatrix} \implies \norm{x^{(1)} - x^{(2)}} = \begin{Vmatrix}
                    x_1^{(1)} - x_1^{(2)} \\ x_2^{(1)} - x_2^{(2)} \\ x_3^{(1)} - x_3^{(2)} \\ x_4^{(1)} - x_4^{(2)}
                \end{Vmatrix} \\ \\
                f(x^{(1)}) = \begin{bmatrix}
                    -2x_1^{(1)} - \p{x_4^{(1)}}^2 + 5\p{x_2^{(1)}}^2 + \cos\p{x_1^{(1)}} \\ -x_2^{(1)} + x_3^{(1)} \\ -4\p{x_3^{(1)}}^2 + \sin\p{x_2^{(1)}} \\ -x_4^{(1)}
                \end{bmatrix} \\ \\ \textrm{and} \\ \\
                f(x^{(2)}) = \begin{bmatrix}
                    -2x_1^{(2)} - \p{x_4^{(2)}}^2 + 5\p{x_2^{(2)}}^2 + \cos\p{x_1^{(2)}} \\ -x_2^{(2)} + x_3^{(2)} \\ -4\p{x_3^{(2)}}^2 + \sin\p{x_2^{(2)}} \\ -x_4^{(2)}
                \end{bmatrix}
            \end{gather*}
            \begin{gather*}
                \implies \norm{f(x^{(1)}) - f(x^{(2)})} = \\ \\
                \begin{Vmatrix}
                    -2x_1^{(1)} - \p{x_4^{(1)}}^2 + 5\p{x_2^{(1)}}^2 + \cos\p{x_1^{(1)}} + 2x_1^{(2)} + \p{x_4^{(2)}}^2 - 5\p{x_2^{(2)}}^2 - \cos\p{x_1^{(2)}} \\ 
                    -x_2^{(1)} + x_3^{(1)} + x_2^{(2)} - x_3^{(2)} \\
                    -4\p{x_3^{(1)}}^2 + \sin\p{x_2^{(1)}} + 4\p{x_3^{(2)}}^2 - \sin\p{x_2^{(2)}} \\
                    -x_4^{(1)} + x_4^{(2)}
                \end{Vmatrix} \\
                = \begin{Vmatrix}
                    -2\p{x_1^{(1)}-x_1^{(2)}} - \p{x_4^{(1)}}^2 + \p{x_4^{(2)}}^2 + 5\p{\p{x_2^{(1)}}^2-\p{x_2^{(2)}}^2} + \cos x_1^{(1)}  - \cos x_1^{(2)} \\ 
                    -\p{x_2^{(1)}-x_2^{(2)}} + \p{x_3^{(1)}-x_3^{(2)}} \\
                    -4\p{\p{x_3^{(1)}}^2-\p{x_3^{(2)}}^2} + \sin\p{x_2^{(1)}}  - \sin\p{x_2^{(2)}} \\
                    -\p{x_4^{(1)} - x_4^{(2)}}
                \end{Vmatrix}
            \end{gather*}
            \noindent Using the sum to product trig identities and $\sin\p{\alpha x} \le \alpha x$:
            \begin{align*}
                \abs{\cos\p{x_1^{(1)}} - \cos\p{x_1^{(2)}}} &= \abs{-2\sin\p{\frac{1}{2}\p{x_1^{(1)}+x_1^{(2)}}}\sin\p{\frac{1}{2}\p{x_1^{(1)}-x_1^{(2)}}}} \\
                &\le \abs{2\sin\p{\frac{1}{2}\p{x_1^{(1)}-x_1^{(2)}}}} \\
                &\le \abs{2\cdot\frac{1}{2}\p{x_1^{(1)}-x_1^{(2)}}} \\ \\
                \abs{\sin\p{x_2^{(1)}}-\sin\p{x_2^{(2)}}} &\le \abs{x_2^{(1)}-x_2^{(2)}}
            \end{align*}
        \end{framed}
        \begin{framed}
            \noindent Upper bounding each row in $\norm{f(x^{(1)}) - f(x^{(2)})}$ yields:
            \begin{align*}
                \abs{f_1(x^{(1)}) - f_1(x^{(2)})} &\le \abs{-2\p{x_1^{(1)}-x_1^{(2)}}} + \abs{\p{x_4^{(2)}}^2-\p{x_4^{(1)}}^2} + \abs{5\p{\p{x_2^{(1)}}^2-\p{x_2^{(2)}}^2}} \\
                &\quad\quad + \abs{\p{x_1^{(1)}-x_1^{(2)}}} \\
                &\le 2\abs{x_1^{(1)}-x_1^{(2)}} + \abs{x_4^{(2)}-x_4^{(1)}}\cdot\abs{x_4^{(2)}+x_4^{(1)}} + 5\abs{x_2^{(1)}-x_2^{(2)}}\cdot\abs{x_2^{(1)}+x_2^{(2)}} \\
                &\quad\quad + \abs{x_1^{(1)}-x_1^{(2)}} \\
                &\le 2\abs{x_1^{(1)}-x_1^{(2)}} + \abs{x_4^{(2)}-x_4^{(1)}}\cdot4 + 5\abs{x_2^{(1)}-x_2^{(2)}}\cdot4 + \abs{x_1^{(1)}-x_1^{(2)}} \\
                &\le 3\abs{x_1^{(1)}-x_1^{(2)}} + 4\abs{x_4^{(2)}-x_4^{(1)}} + 20\abs{x_2^{(1)}-x_2^{(2)}} \\
                \abs{f_2(x^{(1)}) - f_2(x^{(2)})} &\le \abs{-\p{x_2^{(1)}-x_2^{(2)}}} + \abs{\p{x_3^{(1)}-x_3^{(2)}}} \\
                &\le \abs{x_2^{(1)}-x_2^{(2)}} + \abs{x_3^{(1)}-x_3^{(2)}} \\
                \abs{f_3(x^{(1)}) - f_3(x^{(2)})} &\le \abs{-4\p{\p{x_3^{(1)}}^2-\p{x_3^{(2)}}^2}} + \abs{x_2^{(1)}-x_2^{(2)}} \\ 
                &\le 4\abs{x_3^{(1)}-x_3^{(2)}}\cdot\abs{x_3^{(1)}+x_3^{(2)}} + \abs{x_2^{(1)}-x_2^{(2)}} \\ 
                &\le 16\abs{x_3^{(1)}-x_3^{(2)}} + \abs{x_2^{(1)}-x_2^{(2)}} \\ 
                \abs{f_4(x^{(1)}) - f_4(x^{(2)})} &\le \abs{x_4^{(1)} - x_4^{(2)}}
            \end{align*}
            \noindent Defining $\Delta x_i = x_i^{(1)} - x_i^{(2)}$ and rewriting:
            \begin{gather*}
                \norm{f(x^{(1)}) - f(x^{(2)})} \le \begin{Vmatrix}
                    3\abs{\Delta x_1} + 4\abs{\Delta x_4} + 20\abs{\Delta x_2} \\
                    \abs{\Delta x_2} + \abs{\Delta x_3} \\
                    16\abs{\Delta x_3} + \abs{\Delta x_2} \\
                    \abs{\Delta x_4}
                \end{Vmatrix}
            \end{gather*}
            \noindent Use the 2-norm:
            \begin{align*}
                ||f(&x^{(1)}) - f(x^{(2)})||_2 = \\
                &\quad \sqrt{\p{3\abs{\Delta x_1} + 4\abs{\Delta x_4} + 20\abs{\Delta x_2}}^2 + \p{\abs{\Delta x_2} + \abs{\Delta x_3}}^2 + \p{16\abs{\Delta x_3} + \abs{\Delta x_2}}^2 + \abs{\Delta x_4}^2}
            \end{align*}
            \noindent Square both sides:
            \begin{align*}
                ||f(&x^{(1)}) - f(x^{(2)})||_2^2 \le \\
                &\quad \sqrt{\p{3\abs{\Delta x_1} + 4\abs{\Delta x_4} + 20\abs{\Delta x_2}}^2 + \p{\abs{\Delta x_2} + \abs{\Delta x_3}}^2 + \p{16\abs{\Delta x_3} + \abs{\Delta x_2}}^2 + \abs{\Delta x_4}^2}
            \end{align*}
            \noindent Use dot product inequality to bound RHS ($\norm{a\cdot b}^2 \le \norm{a}^2\norm{b}^2$):
            \begin{align*}
                \norm{f(x^{(1)}) - f(x^{(2)})}_2^2 &\le \p{3^2+20^2+4^2}\p{\Delta x_1^2 + \Delta x_2^2 + \Delta x_4^2} \\
                &\quad\quad + \p{1^2 + 1^2}\p{\Delta x_2^2 + \Delta x_3^2} + \p{1^2 + 16^2}\p{\Delta x_2^2 + \Delta x_3^2} + \Delta x_4^2 \\
                &= 425\p{\Delta x_1^2 + \Delta x_2^2 + \Delta x_4^2} + 2\Delta x_2^2 + 2\Delta x_3^2 + 257\p{\Delta x_2^2 + \Delta x_3^2} + \Delta x_4^2 \\
                &= 425\Delta x_1^2 + 684\Delta x_2^2 + 259\Delta x_3^2 + 426\Delta x_4^2
            \end{align*}
            \noindent Take the square root:
            \begin{align*}
                \norm{f(x^{(1)}) - f(x^{(2)})}_2 &\le \sqrt{425\Delta x_1^2 + 684\Delta x_2^2 + 259\Delta x_3^2 + 426\Delta x_4^2} \\
                &\le \sqrt{684}\norm{\Delta x}_2 = K\norm{\Delta x}_2 = K\norm{x^{(1)} - x^{(2)}} \\ \\
                \implies K &= \sqrt{684} \approx 26.153
            \end{align*}
        \end{framed}
% 1b
        \item \hfill
        \begin{framed}
            \noindent To compute a Lipschitz constant using the Jacobian, we must have $\norm{J(x)}\le K$. For simplicity, I am choosing to use the Frobenius norm. \\ \\
            \noindent Construct $J(x)$:
            \begin{gather*}
                \frac{\partial f_1}{\partial x_1} = -2-\sin(x_1) \quad\quad \frac{\partial f_1}{\partial x_2} =10x_2 \quad\quad \frac{\partial f_1}{\partial x_3} = 0 \quad\quad \frac{\partial f_1}{\partial x_4} = -2x_4 \\
                \frac{\partial f_2}{\partial x_1} = \frac{\partial f_2}{\partial x_4} = 0 \quad\quad \frac{\partial f_2}{\partial x_2} = -1 \quad\quad \frac{\partial f_2}{\partial x_3} = 1 \\
                \frac{\partial f_3}{\partial x_1} = \frac{\partial f_3}{\partial x_4} = 0 \quad\quad \frac{\partial f_3}{\partial x_2} = \cos(x_2) \quad\quad \frac{\partial f_3}{\partial x_3} = -8x_3 \\
                \frac{\partial f_4}{\partial x_1} = \frac{\partial f_4}{\partial x_2} = \frac{\partial f_4}{\partial x_3} = 0 \quad\quad \frac{\partial f_4}{\partial x_4} = -1
            \end{gather*}
            \begin{gather*}
                \implies J(x) = \begin{bmatrix}
                    \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \frac{\partial f_1}{\partial x_3} & \frac{\partial f_1}{\partial x_4} \\
                    \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \frac{\partial f_2}{\partial x_3} & \frac{\partial f_2}{\partial x_4} \\
                    \frac{\partial f_3}{\partial x_1} & \frac{\partial f_3}{\partial x_2} & \frac{\partial f_3}{\partial x_3} & \frac{\partial f_3}{\partial x_4} \\
                    \frac{\partial f_4}{\partial x_1} & \frac{\partial f_4}{\partial x_2} & \frac{\partial f_4}{\partial x_3} & \frac{\partial f_4}{\partial x_4} \\
                \end{bmatrix} = \begin{bmatrix}
                    -2-\sin(x_1) & 10x_2 & 0 & -2x_4 \\
                    0 & -1 & 1 & 0 \\
                    0 & \cos(x_2) & -8x_3 & 0 \\
                    0 & 0 & 0 & -1
                \end{bmatrix}
            \end{gather*}
            $\norm{J}_F$ is upper bounded by the Frobenius norm of the matrix containing the maximum absolute values of each entry of $J$. Keeping in mind that $x \in \mathcal{X} = \{0\le x_i\le2,i=1,2,3,4\}$, we now have:
            \begin{gather*}
                \norm{J}_F = \begin{Vmatrix}
                    -2-\sin(x_1) & 10x_2 & 0 & -2x_4 \\
                    0 & -1 & 1 & 0 \\
                    0 & \cos(x_2) & -8x_3 & 0 \\
                    0 & 0 & 0 & -1
                \end{Vmatrix}_F \le \begin{Vmatrix}
                    2+\sin(2) & 20 & 0 & 4 \\
                    0 & 1 & 1 & 0 \\
                    0 & \cos(0) & 16 & 0 \\
                    0 & 0 & 0 & 1
                \end{Vmatrix}_F \\
            \end{gather*}
            \noindent Computing the norm yields:
            \begin{gather*}
                \norm{J}_F \le \sqrt{\p{2+\sin(2)}^2 + 20^2 + 4^2 + 1^2 + 1^2 + \cos^2(0) + 16^2 + 1^1} \approx 26.162 \le K \\ \\
                \implies K = 26.162
            \end{gather*} 
        \end{framed}
\newpage
% 1c
        \item \hfill
        \begin{framed}
            \noindent In order to numerically compute Lipschitz constants, I wrote a MATLAB script that takes different resolutions for evenly spaced sampling for $x(t)\in\mathcal{X}$, loops over a 4D grid of points (one for each $x_i$) with the given resolution, evaluates the Jacobian at the sampling point, computes $K$ by taking the 2-norm of the Jacobian, checks to see if that $K$-value is a new maximum, stores it if it is, repeats for all the sampling points and for all resolutions, and finally plots the results of $K$ for various resolutions $h$. \\ \\
            \noindent From Figure \ref{fig:K_vs_resolution}, which shows the results from the procedure above, we see that $K$ seems to be around 20.5 to 20.7, which matches our analytical results from parts (a) and (b). However, there is a dip for step sizes $h=0.8$ and $h=0.75$; this is an example of how sampling only a select number of points can miss other points that yield larger $K$-values. 
            \begin{center}
                \includegraphics[width=0.75\linewidth]{problem_1c_numerical_Lipschitz_constant_step_size_resolution.eps}
                \captionof{figure}{Plot of numerical $K$ versus resolution $h$.}
                \label{fig:K_vs_resolution}
            \end{center}
            \noindent For fun, I also added a section in my code that randomly samples points from $\mathcal{X}$ instead of sampling evenly-spaced points. Figure \ref{fig:K_vs_random_sampling} shows $K$ plotted against the number of randomly sampled points. From the plot, we see how important it is to have a large sampling size (showing convergence to a specific $K$), but it also shows how drastically increasing the sample size (i.e., from $10^4$ to $5\cdot 10^4$) may not lead to a much different $K$.
            \begin{center}
                \includegraphics[width=0.75\linewidth]{problem_1c_numerical_Lipschitz_constant_random_sampling.eps}
                \captionof{figure}{Plot of numerical $K$ versus number of randomly sampled points.}
                \label{fig:K_vs_random_sampling}
            \end{center}
        \end{framed}
        \begin{framed}
            \noindent My MATLAB script for problem 1 part (c) is shown below:
        \end{framed}
        \begin{center}
            \lstinputlisting[frame=single, numbers=left, style=Matlab-Pyglike]{problem1.m}
        \end{center}
    \end{enumerate}
\newpage
% Problem 2 begins here
    \item \hfill
    \begin{enumerate}
% 2a
        \item \hfill
        \begin{framed}
            \noindent The given system of linear ODEs can be written as:
            \begin{gather*}
                \dot{x} = Ax(t) + Bu(t) \\ \\
                \begin{bmatrix}
                    \dot{x}_1(t) \\ \dot{x}_2(t) \\ \dot{x}_3(t) \\ \dot{x}_4(t)
                \end{bmatrix} = \begin{bmatrix}
                    -2 & 4 & 3 & 0 \\
                    0 & -1 & 0 & 4 \\
                    0 & 0 & -3 & 0 \\
                    0 & 0 & 0 & -8
                \end{bmatrix} \begin{bmatrix}
                    x_1(t) \\ x_2(t) \\ x_3(t) \\ x_4(t)
                \end{bmatrix} + \begin{bmatrix}
                    1 & 0 \\
                    0 & 0 \\
                    0 & 1 \\
                    0 & 0 
                \end{bmatrix} \begin{bmatrix}
                    u_1(t) \\ u_2(t)
                \end{bmatrix}
            \end{gather*}
        \end{framed}
% 2b
        \item \hfill
        \begin{framed}
            \noindent Since $A$ is an upper triangular matrix, the eigenvalues of $A$ are its diagonal entries:
            \begin{gather*}
                \textrm{eig}\p{A} = \lambda = \begin{bmatrix}
                    -1 \\ -2 \\ -3 \\ -8
                \end{bmatrix}
            \end{gather*}
            \noindent Compute eigenvectors $\p{A - \lambda I_n}v = 0$: \\ \\
            \noindent For $\lambda_1 = -1$:
            \begin{gather*}
                \begin{bmatrix}
                    -2-(-1) & 4 & 3 & 0 \\
                    0 & -1-(-1) & 0 & 4 \\
                    0 & 0 & -3-(-1) & 0 \\
                    0 & 0 & 0 & -8-(-1)
                \end{bmatrix} \begin{bmatrix}
                    v_a \\ v_b \\v_c \\ v_d
                \end{bmatrix} = 0 \\ \\
                \implies \begin{cases}
                    -v_a + 4v_b + 3v_c + 0 = 0 \\
                    0 + 0 + 0 + 4v_d = 0 \\
                    0 + 0 - 2v_c + 0 = 0 \\
                    0 + 0 + 0 - 7v_d = 0 \\
                \end{cases} \\
                \implies v_c = v_d = 0 \implies \textrm{choose } v_b = 1 \implies v_a = 4 \\ \\
                \implies v_1 = \begin{bmatrix}
                    4 \\ 1 \\ 0 \\ 0
                \end{bmatrix}
            \end{gather*}
            \noindent For $\lambda_2 = -2$:
            \begin{gather*}
                \begin{bmatrix}
                    -2-(-2) & 4 & 3 & 0 \\
                    0 & -1-(-2) & 0 & 4 \\
                    0 & 0 & -3-(-2) & 0 \\
                    0 & 0 & 0 & -8-(-2)
                \end{bmatrix} \begin{bmatrix}
                    v_a \\ v_b \\v_c \\ v_d
                \end{bmatrix} = 0 \\ \\
                \implies \begin{cases}
                    0 + 4v_b + 3v_c + 0 = 0 \\
                    0 + v_b + 0 + 4v_d = 0 \\
                    0 + 0 - 1v_c + 0 = 0 \\
                    0 + 0 + 0 - 6v_d = 0 \\
                \end{cases} \\
                \implies v_c = v_d = 0 \implies v_b = 0 \implies \textrm{choose } v_a = 1 \\ \\
                \implies v_2 = \begin{bmatrix}
                    1 \\ 0 \\ 0 \\ 0
                \end{bmatrix}
            \end{gather*}
            \noindent For $\lambda_3 = -3$:
            \begin{gather*}
                \begin{bmatrix}
                    -2-(-3) & 4 & 3 & 0 \\
                    0 & -1-(-3) & 0 & 4 \\
                    0 & 0 & -3-(-3) & 0 \\
                    0 & 0 & 0 & -8-(-3)
                \end{bmatrix} \begin{bmatrix}
                    v_a \\ v_b \\v_c \\ v_d
                \end{bmatrix} = 0 \\ \\
                \implies \begin{cases}
                    v_a + 4v_b + 3v_c + 0 = 0 \\
                    0 + 2v_b + 0 + 4v_d = 0 \\
                    0 + 0 + 0 + 0 = 0 \\
                    0 + 0 + 0 - 5v_d = 0 \\
                \end{cases} \\
                \implies v_d = 0 \implies v_b = 0 \implies \textrm{choose } v_c = 1 \implies v_a = -3 \\ \\
                \implies v_3 = \begin{bmatrix}
                    -3 \\ 0 \\ 1 \\ 0
                \end{bmatrix}
            \end{gather*}
            \noindent For $\lambda_4 = -8$:
            \begin{gather*}
                \begin{bmatrix}
                    -2-(-8) & 4 & 3 & 0 \\
                    0 & -1-(-8) & 0 & 4 \\
                    0 & 0 & -3-(-8) & 0 \\
                    0 & 0 & 0 & -8-(-8)
                \end{bmatrix} \begin{bmatrix}
                    v_a \\ v_b \\v_c \\ v_d
                \end{bmatrix} = 0 \\ \\
                \implies \begin{cases}
                    6v_a + 4v_b + 3v_c + 0 = 0 \\
                    0 + 7v_b + 0 + 4v_d = 0 \\
                    0 + 0 + 5v_c + 0 = 0 \\
                    0 + 0 + 0 + 0 = 0 \\
                \end{cases} \\
                \implies v_c = 0 \implies \textrm{choose } v_d = 1 \implies v_b = -4/7 \implies v_a = 8/21 \\ \\
                \implies v_4 = \begin{bmatrix}
                    8/21 \\ -4/7 \\ 0 \\ 1
                \end{bmatrix}
            \end{gather*}
            \noindent Construct $T$ and compute $T^{-1}$ using Gauss-Jordan elimination:
            \begin{gather*}
                T = \begin{bmatrix}
                    v_1 & v_2 & v_3 & v_4
                \end{bmatrix} = \begin{bmatrix}
                    4 & 1 & -3 & 8/21 \\
                    1 & 0 & 0 & -4/7 \\
                    0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1
                \end{bmatrix}
            \end{gather*}
            \begin{gather*}
                \begin{bmatrix}
                    4 & 1 & -3 & 8/21 & | & 1 & 0 & 0 & 0 \\
                    1 & 0 & 0 & -4/7 & | & 0 & 1 & 0 & 0 \\
                    0 & 0 & 1 & 0 & | & 0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1 & | & 0 & 0 & 0 & 1
                \end{bmatrix} \ro{R_1 \leftarrow \frac{1}{4}R_1} \begin{bmatrix}
                    1 & 1/4 & -3/4 & 2/21 & | & 1/4 & 0 & 0 & 0 \\
                    1 & 0 & 0 & -4/7 & | & 0 & 1 & 0 & 0 \\
                    0 & 0 & 1 & 0 & | & 0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1 & | & 0 & 0 & 0 & 1
                \end{bmatrix} \\ \\
                \ro{R_2 \leftarrow R_2 - R_1} \begin{bmatrix}
                    1 & 1/4 & -3/4 & 2/21 & | & 1/4 & 0 & 0 & 0 \\
                    0 & -1/4 & 3/4 & -2/3 & | & -1/4 & 1 & 0 & 0 \\
                    0 & 0 & 1 & 0 & | & 0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1 & | & 0 & 0 & 0 & 1
                \end{bmatrix} \\ \\
                \ro{R_2 \leftarrow -4R_2}
                \begin{bmatrix}
                    1 & 1/4 & -3/4 & 2/21 & | & 1/4 & 0 & 0 & 0 \\
                    0 & 1 & -3 & 8/3 & | & 1 & -4 & 0 & 0 \\
                    0 & 0 & 1 & 0 & | & 0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1 & | & 0 & 0 & 0 & 1
                \end{bmatrix}
            \end{gather*}
            \begin{gather*}
                \ro{R_2 \leftarrow R_2 - \frac{8}{3}R_4} \begin{bmatrix}
                    1 & 1/4 & -3/4 & 2/21 & | & 1/4 & 0 & 0 & 0 \\
                    0 & 1 & -3 & 0 & | & 1 & -4 & 0 & -8/3 \\
                    0 & 0 & 1 & 0 & | & 0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1 & | & 0 & 0 & 0 & 1
                \end{bmatrix} \\ \\
                \ro{R_1 \leftarrow R_1 - \frac{2}{21}R_4} \begin{bmatrix}
                    1 & 1/4 & -3/4 & 0 & | & 1/4 & 0 & 0 & -2/21 \\
                    0 & 1 & -3 & 0 & | & 1 & -4 & 0 & -8/3 \\
                    0 & 0 & 1 & 0 & | & 0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1 & | & 0 & 0 & 0 & 1
                \end{bmatrix} \\ \\
                \ro{R_2 \leftarrow R_2 + 3R_3} \begin{bmatrix}
                    1 & 1/4 & -3/4 & 0 & | & 1/4 & 0 & 0 & -2/21 \\
                    0 & 1 & 0 & 0 & | & 1 & -4 & 3 & -8/3 \\
                    0 & 0 & 1 & 0 & | & 0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1 & | & 0 & 0 & 0 & 1
                \end{bmatrix} \\ \\
                \ro{R_1 \leftarrow R_1 + \frac{3}{4}R_3} \begin{bmatrix}
                    1 & 1/4 & 0 & 0 & | & 1/4 & 0 & 3/4 & -2/21 \\
                    0 & 1 & 0 & 0 & | & 1 & -4 & 3 & -8/3 \\
                    0 & 0 & 1 & 0 & | & 0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1 & | & 0 & 0 & 0 & 1
                \end{bmatrix} \\ \\
                \ro{R_1 \leftarrow R_1 - \frac{1}{4}R_2} \begin{bmatrix}
                    1 & 0 & 0 & 0 & | & 0 & 1 & 0 & 4/7 \\
                    0 & 1 & 0 & 0 & | & 1 & -4 & 3 & -8/3 \\
                    0 & 0 & 1 & 0 & | & 0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1 & | & 0 & 0 & 0 & 1
                \end{bmatrix} \\ \\
                \implies T^{-1} = \begin{bmatrix}
                    0 & 1 & 0 & 4/7 \\
                    1 & -4 & 3 & -8/3 \\
                    0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1
                \end{bmatrix}
            \end{gather*}
            \noindent Putting everything together:
            \begin{gather*}
                A = TDT^{-1} = \begin{bmatrix}
                    4 & 1 & -3 & 8/21 \\
                    1 & 0 & 0 & -4/7 \\
                    0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1
                \end{bmatrix} \begin{bmatrix}
                    -1 & 0 & 0 & 0\\
                    0 & -2 & 0 & 0 \\
                    0 & 0 & -3 & 0 \\
                    0 & 0 & 0 & -8
                \end{bmatrix} \begin{bmatrix}
                    0 & 1 & 0 & 4/7 \\
                    1 & -4 & 3 & -8/3 \\
                    0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1
                \end{bmatrix}
            \end{gather*}
        \end{framed}
\newpage
% 2c
        \item \hfill
        \begin{framed}
            \noindent Compute $e^{At}$:
            \begin{align*}
                e^{At} &= Te^{Dt}T^{-1} \\
                &= \begin{bmatrix}
                    4 & 1 & -3 & 8/21 \\
                    1 & 0 & 0 & -4/7 \\
                    0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1
                \end{bmatrix} \begin{bmatrix}
                    e^{-t} & 0 & 0 & 0\\
                    0 & e^{-2t} & 0 & 0 \\
                    0 & 0 & e^{-3t} & 0 \\
                    0 & 0 & 0 & e^{-8t}
                \end{bmatrix} \begin{bmatrix}
                    0 & 1 & 0 & 4/7 \\
                    1 & -4 & 3 & -8/3 \\
                    0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1
                \end{bmatrix}
            \end{align*}
        \end{framed}
% 2d
        \item \hfill
        \begin{framed}
            \noindent Solve the linear system of ODEs analytically, given:
            \begin{gather*}
                x(t_0) = x(0) = \begin{bmatrix}
                    1 & -1 & 2 & 0
                \end{bmatrix}^\top \quad\quad u(t) = \begin{bmatrix}
                    u_1(t) \\ u_2(t)
                \end{bmatrix} = \begin{bmatrix}
                    1-e^{2t} \\ 2\cos(t)
                \end{bmatrix}
            \end{gather*}
            \noindent To analytically solve, we must compute $x(t) = e^{At-t_0}x(t_0) + \int_{t_0}^{t}e^{A(t-\tau)}Bu(\tau)d\tau$. \\ \\
            \noindent First, compute $e^{At-t_0}x(t_0)$:
            \begin{align*}
                e^{At-t_0}x(t_0) &= e^{At-0}x(0) = e^{At}x(0) \\
                &= \begin{bmatrix}
                    4 & 1 & -3 & 8/21 \\
                    1 & 0 & 0 & -4/7 \\
                    0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1
                \end{bmatrix} \begin{bmatrix}
                    e^{-t} & 0 & 0 & 0\\
                    0 & e^{-2t} & 0 & 0 \\
                    0 & 0 & e^{-3t} & 0 \\
                    0 & 0 & 0 & e^{-8t}
                \end{bmatrix} \begin{bmatrix}
                    0 & 1 & 0 & 4/7 \\
                    1 & -4 & 3 & -8/3 \\
                    0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1
                \end{bmatrix} \begin{bmatrix}
                    1 \\ -1 \\ 2 \\ 0
                \end{bmatrix} \\
                &= \begin{bmatrix}
                    4 & 1 & -3 & 8/21 \\
                    1 & 0 & 0 & -4/7 \\
                    0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1
                \end{bmatrix} \begin{bmatrix}
                    e^{-t} & 0 & 0 & 0\\
                    0 & e^{-2t} & 0 & 0 \\
                    0 & 0 & e^{-3t} & 0 \\
                    0 & 0 & 0 & e^{-8t}
                \end{bmatrix} \begin{bmatrix}
                    -1 \\ 1+4+6=11 \\2 \\ 0
                \end{bmatrix} \\
                &= \begin{bmatrix}
                    4 & 1 & -3 & 8/21 \\
                    1 & 0 & 0 & -4/7 \\
                    0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1
                \end{bmatrix} \begin{bmatrix}
                    -e^{-t} \\ 11e^{-2t} \\ 2e^{-3t} \\ 0
                \end{bmatrix} \\
                &= \begin{bmatrix}
                    -4e^{-t} + 11e^{-2t} - 6e^{-3t} \\ -e^{-t} \\ 2e^{-3t} \\ 0
                \end{bmatrix}
            \end{align*}
            \noindent Now, compute the integrand:
            \begin{align*}
                e^{A(t-\tau)}Bu(\tau) &= T \begin{bmatrix}
                    e^{-(t-\tau)} & 0 & 0 & 0\\
                    0 & e^{-2(t-\tau)} & 0 & 0 \\
                    0 & 0 & e^{-3(t-\tau)} & 0 \\
                    0 & 0 & 0 & e^{-8(t-\tau)}
                \end{bmatrix} \begin{bmatrix}
                    0 & 1 & 0 & 4/7 \\
                    1 & -4 & 3 & -8/3 \\
                    0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1
                \end{bmatrix} \begin{bmatrix}
                    1-e^{2\tau} \\ 0 \\ 2\cos(\tau) \\ 0
                \end{bmatrix} \\
                &= T \begin{bmatrix}
                    e^{-(t-\tau)} & 0 & 0 & 0\\
                    0 & e^{-2(t-\tau)} & 0 & 0 \\
                    0 & 0 & e^{-3(t-\tau)} & 0 \\
                    0 & 0 & 0 & e^{-8(t-\tau)}
                \end{bmatrix} \begin{bmatrix}
                    0 \\ 1-e^{2\tau}+6\cos(\tau) \\ 2\cos(\tau) 0
                \end{bmatrix} \\
                &= T \begin{bmatrix}
                    0 \\ e^{-2(t-\tau)} - e^{2\tau}e^{-2(t-\tau)} + 6\cos(\tau)e^{-2(t-\tau)} \\ 2\cos(\tau)e^{-3(t-\tau)}
                \end{bmatrix} \\
                &= T \begin{bmatrix}
                    0 \\ e^{-2(t-\tau)}(1-e^{2\tau}+6\cos(\tau)) \\ 2\cos(\tau)e^{-3(t-\tau)} \\ 0
                \end{bmatrix}
            \end{align*}
            \noindent Since the leading matrix $(T)$ is a constant, I will factor it out of the integral and multiply it back in after evaluating the integral.
            \begin{align*}
                \int_{t_0}^{t}e^{A(t-\tau)}Bu(\tau)d\tau &= T \int_{0}^{t} \begin{bmatrix}
                    0 \\ e^{-2(t-\tau)}(1-e^{2\tau}+6\cos(\tau)) \\ 2\cos(\tau)e^{-3(t-\tau)} \\ 0
                \end{bmatrix} d\tau
            \end{align*}
            \noindent Evaluating each row separately yields:
            \begin{align*}
                I_2 &= \int_0^t \p{e^{-2(t-\tau)}(1-e^{2\tau}+6\cos(\tau))}d\tau \\
                &= e^{-2t}\int_0^t e^{2\tau}\p{1-e^{2\tau}+6\cos(\tau)}d\tau \\
                &= e^{-2t}\int_0^t \p{e^{2\tau}-e^{4\tau}+6e^{2\tau}\cos(\tau)}d\tau \\
                &= e^{-2t} \bracket{\frac{1}{2}e^{2\tau} - \frac{1}{4}e^{4\tau} + \frac{6e^{2\tau}}{2^2+1^2}\p{2\cos(\tau)+\sin(\tau)}}\at{0}{t} \\
                &= e^{-2\tau} \bracket{\p{\frac{1}{2}e^{2t}-\frac{1}{2}} - \frac{1}{4}\p{e^{4t}-1} + \frac{6}{5}\p{e^{2t}\p{2\cos(t)+\sin(t)}-\p{2+0}}} \\
                &= e^{-2t} \bracket{\frac{1}{2}e^{2t} -\frac{1}{2} -\frac{1}{4}e^{4t} +\frac{1}{4} + \frac{5}{12}e^{2t}\cos(t) + \frac{5}{6}e^{2t}\sin(t) -\frac{12}{5}} \\
                &= e^{-2t} \bracket{-\frac{53}{20} + \frac{1}{2}e^{2t} - \frac{1}{4}e^{4t} + \frac{12}{5}e^{2t}\cos(t) + \frac{6}{5}e^{2t}\sin(t)} \\
                &= -\frac{53}{20}e^{-2t} + \frac{1}{2} - \frac{1}{4}e^{2t} + \frac{12}{5}\cos(t) + \frac{6}{5}\sin(t)
            \end{align*}
            \begin{align*}
                I_3 &= \int_{0}^{t} 2\cos(\tau)e^{-3(t-\tau)}d\tau \\
                &= 2e^{-3t}\int_0^t e^{3\tau}\cos(\tau)d\tau \\
                &= 2e^{-3t} \bracket{\frac{e^{3\tau}}{3^2 + 1^2}\p{3\cos(\tau)+\sin(\tau)}}\at{0}{t} \\
                &= \frac{1}{5}e^{-3t}\bracket{e^{3t}\p{3\cos(t)+\sin(t)}-3-0} \\
                &= \frac{1}{5}\p{3\cos(t)+\sin(t)}-\frac{3}{5}e^{-3t} \\ \\
            I_1 &= I_4 = 0
            \end{align*}
            \noindent Multiply by $T$ and simplify:
            \begin{gather*}
                \int_{t_0}^{t}e^{A(t-\tau)}Bu(\tau)d\tau = \begin{bmatrix}
                    4 & 1 & -3 & 8/21 \\
                    1 & 0 & 0 & -4/7 \\
                    0 & 0 & 1 & 0 \\
                    0 & 0 & 0 & 1
                \end{bmatrix} \begin{bmatrix}
                    0 \\
                    -\frac{53}{20}e^{-2t} + \frac{1}{2} - \frac{1}{4}e^{2t} + \frac{12}{5}\cos(t) + \frac{6}{5}\sin(t) \\
                    \frac{1}{5}\p{3\cos(t)+\sin(t)}-\frac{3}{5}e^{-3t} \\
                    0
                \end{bmatrix} \\
                = \begin{bmatrix}
                    -\frac{53}{20}e^{-2t} + \frac{1}{2} - \frac{1}{4}e^{2t} + \frac{12}{5}\cos(t) + \frac{6}{5}\sin(t) - \frac{3}{5}\p{3\cos(t)+\sin(t)} + \frac{9}{5}e^{-3t} \\
                    0 \\
                    \frac{1}{5}\p{3\cos(t) + \sin(t)} - \frac{3}{5}e^{-3t} \\
                    0
                \end{bmatrix}
            \end{gather*}
            \noindent Putting everything together yields our solution:
            \begin{align*}
                x(t) &= \begin{bmatrix}
                    -4e^{-t} + 11e^{-2t} - 6e^{-3t} \\ -e^{-t} \\ 2e^{-3t} \\ 0
                \end{bmatrix} + \begin{bmatrix}
                    \frac{1}{2} + \frac{9}{5}e^{-3t} - \frac{53}{20}e^{-2t} - \frac{1}{4}e^{2t} + \frac{3}{5}\p{\cos(t)+\sin(t)} \\
                    0 \\
                    \frac{1}{5}\p{3\cos(t) + \sin(t)} - \frac{3}{5}e^{-3t} \\
                    0
                \end{bmatrix}
            \end{align*}
            \begin{gather*}
                \implies x(t) = \begin{bmatrix}
                    \frac{1}{2} - \frac{21}{5}e^{-3t} + \frac{167}{20}e^{-2t} - 4e^{-t} - \frac{1}{4}e^{2t} + \frac{3}{5}\p{\cos(t)+\sin(t)} \\
                    -e^{-t} \\
                    \frac{1}{5}\p{3\cos(t) + \sin(t)} + \frac{7}{5}e^{-3t} \\
                    0
                \end{bmatrix}
            \end{gather*}
        \end{framed}
\newpage
% 2e
        \item \hfill
        \begin{framed}
            \noindent Figure \ref{fig:2e_analytical_plots} contains my plots for each of the four states over a time span of $t\in\bracket{0,5}$ seconds with an increment of $h=0.01$ seconds.
            \begin{center}
                \includegraphics[width=1\linewidth]{problem_2e_analytical_solution_plots.eps}
                \captionof{figure}{Plots of the analytical solutions for each state over $t\in\bracket{0,2}$ seconds and $h=0.01$ seconds.}
                \label{fig:2e_analytical_plots}
            \end{center}
        \end{framed}
        \begin{framed}
            \noindent Below is my MATLAB code to produce these plots.
        \end{framed}
        \begin{center}
            \lstinputlisting[frame=single, numbers=left, style=Matlab-Pyglike]{problem_2e.m}
        \end{center}
% 2f
        \item \hfill
        \begin{framed}
            \noindent Figure \ref{fig:2f_numerical_plots} contains my plots for the numerical solutions of the four states over a time span of $t\in\bracket{0,5}$ seconds with an increment of $h=0.01$ seconds.
            \begin{center}
                \includegraphics[width=1\linewidth]{problem_2f_numerical_solution_plots.eps}
                \captionof{figure}{Plots of the numerical solutions for each state over $t\in\bracket{0,2}$ seconds and $h=0.01$ seconds.}
                \label{fig:2f_numerical_plots}
            \end{center}
        \end{framed}
        \begin{framed}
            \noindent Below is my MATLAB code to produce these plots.
        \end{framed}
        \begin{center}
            \lstinputlisting[frame=single, numbers=left, style=Matlab-Pyglike]{problem_2f.m}
        \end{center}
\newpage
% 2g
        \item \hfill
        \begin{framed}
            \noindent For ease of viewing, I have overlaid my plots for the analytical and numerical solutions for the system of ODEs in Figure \ref{fig:2g_overlaid_plots}. My plots from parts (e) and (f) look nearly identical.
            \begin{center}
                \includegraphics[width=1\linewidth]{problem_2g_overlaid_plots.eps}
                \captionof{figure}{Overlaid plots of the analytical (black, dashed) and numerical (colored, solid) solutions for each state over $t\in\bracket{0,2}$ seconds and $h=0.01$ seconds.}
                \label{fig:2g_overlaid_plots}
            \end{center}
        \end{framed}
        \begin{framed}
            \noindent My MATLAB code to overlay these plots is reproduced below.
        \end{framed}
        \begin{center}
            \lstinputlisting[frame=single, numbers=left, style=Matlab-Pyglike]{problem_2g.m}
        \end{center}
\newpage
% 2h
        \item \hfill
        \begin{framed}
            \noindent Using the \textit{ode45} function on MATLAB for the given system of ODEs and then plotting these states yields the following plots in Figure \ref{fig:2h_ode45_plots}
            \begin{center}
                \includegraphics[width=1\linewidth]{problem_2h_ODE45_solution_plots.eps}
                \captionof{figure}{Plots of the \textit{ode45} numerical solutions for each state over $t\in\bracket{0,2}$ seconds and $h=0.01$ seconds.}
                \label{fig:2h_ode45_plots}
            \end{center}
        \end{framed}
        \begin{framed}
            \noindent Below is my MATLAB code using the \textit{ode45} solver and to produce these plots.
        \end{framed}
        \begin{center}
            \lstinputlisting[frame=single, numbers=left, style=Matlab-Pyglike]{problem_2h.m}
        \end{center}
    \end{enumerate}
\newpage
% Problem 3 begins here
    \item \hfill
    \begin{enumerate}
% 3a
        \item \hfill
        \begin{framed}
            \noindent Derive the Adam-Bashforth four-step explicit method, given by:
            \begin{gather*}
                y_{i+1} = y_i + \frac{h}{24}\p{55f_i - 59f_{i-1} + 37f_{i-2} - 9f_{i-3}}.
            \end{gather*}
            \noindent We are given that $\dot{y}=f(t,y),t\in\bracket{a,b}$, and our initialization plus 3 more previous initializations. To derive this four-step explicit method, we must construct the third-degree Lagrange interpolation polynomial, integrate it for a single time step, and then add the result to our most recent initialization $y_i$. This produces a local truncation error of order $h^5$, which produces a global error of $h^4$. My MATLAB script to symbolically derive this method is shown at the end of this problem.
        \end{framed}
% 3b
        \item \hfill
        \begin{framed}
            \noindent Derive the Adam-Moulton two-step implicit method, given by:
            \begin{gather*}
                y_{i+1} = y_i + h\p{\frac{1}{12}f_{i-1} - \frac{2}{3}f_i + \frac{5}{12}f_{i+1}}.
            \end{gather*}
            \noindent Similarly to part (a), we are given an ODE, but only need 3 initializations with one of them being the solution ($y_{i+1}$) we are trying to predict. The same process applies, where we construct a second-degree Lagrange polynomial, integrate it for a single time step, and then add it to our most recent initialization $y_i$. This produces a a local error of order $h^4$, which produces a global error of order $h^3$.
        \end{framed}
        \begin{framed}
            \noindent Below is my MATLAB code for problem 3:
        \end{framed}
        \begin{center}
            \lstinputlisting[frame=single, numbers=left, style=Matlab-Pyglike]{problem3.m}
        \end{center}
        \begin{framed}
            \noindent Below is the output:
        \end{framed}
        \begin{center}
            \lstinputlisting[frame=single, numbers=left, style=Matlab-Pyglike]{problem3_output.m}
        \end{center}
    \end{enumerate}
\newpage
% Problem 4 begins here
    \item \hfill
    \begin{framed}
        \noindent Given data:
        \begin{gather*}
            \dot{h}(t) = -\frac{\pi d^2}{4A(h)}\sqrt{2g(h(t)+e)} \quad\quad h(t_0) = h(0) = 6 \quad\quad d = 0.25 \quad\quad e = 1
        \end{gather*}
        \begin{center}
           \begin{tabular}{|c|c|c|c|c|c|c|c|}
               \hline
               $h\units{m}$ & 6 & 5 & 4 & 3 & 2 & 1 & 0 \\
               \hline
               $A(h)$ in $10^4$ [m$^2$] & 1.17 & 0.97 & 0.67 & 0.45 & 0.32 & 0.18 & 0 \\
               \hline
           \end{tabular}
           \captionof{table}{Given data for the area $A(h(t))$ for different values of $h(t)$.}
           \label{tab:p4table}
       \end{center}
    \end{framed}
    \begin{enumerate}
% 4a
        \item \hfill
        \begin{framed}
            \noindent Using MATLAB's \textit{polyfit} and \textit{polyval} functions, I constructed a fifth degree polynomial for $A(h)$ using least squares regression: $A(h) = -10x^5 + 117.803x^4 - 388.636x^3 + 273.106x^2 + 1816.818x - 1.299$. 
        \end{framed}
% 4b
        \item \hfill
        \begin{framed}
            \noindent Figure \ref{fig:poly_regr} contains my plot of the given data points and my regression-based polynomial for $A(h)$.
            \begin{center}
                \includegraphics[width=0.7\linewidth]{4b_5th_degree_poly_regression.eps}
                \captionof{figure}{Data points (red) and fifth degree polynomial regression of $A(h)$.}
                \label{fig:poly_regr}
            \end{center}
        \end{framed}
% 4c
        \item \hfill
        \begin{framed}
            \noindent Using MATLAB, I plugged in my model for $A(h)$ into the given nonlinear ODE and solved it using the \textit{ode45} solver over a time span of $t\in\bracket{0,19\cdot3600}$ seconds. Plotting the solver's solution yields Figure \ref{fig:4_ode45_solution}.
            \begin{center}
                \includegraphics[width=0.75\linewidth]{4c_ode45_solution.eps}
                \captionof{figure}{Solution to the nonlinear ODE using MATLAB's \textit{ode45} solver.}
                \label{fig:4_ode45_solution}
            \end{center}
        \end{framed}
% 4d
        \item \hfill
        \begin{framed}
            \noindent Since there is a singularity at $h(t) = 0$, I approximated the root (the time when the pond is empty) by solving $h(t^*) \approx \varepsilon = 0.01$ using the given initial condition $h(0) = 6$. Using MATLAB, I set my tolerance equal to $\varepsilon$, found the index where the \textit{ode45} solution first reached my tolerance, and plugged that index back into the \textit{ode45} solution's time vector to find $t^* \approx 67650$ seconds.
        \end{framed}
% 4e
        \item \hfill
        \begin{framed}
            \noindent Trying MATLAB's \textit{dsolve} function to solve the nonlinear ODE yielded no analytical results. This means that the ODE is too complex to analytically solve, so we must resort to numerical methods!
        \end{framed}
        \begin{framed}
            \noindent My MATLAB script for problem 4 parts (a) through (e), followed by its output, is included at the end of this problem. In the output, I received a warning when plotting the \textit{ode45} solution due to the singularity and complex numbers; the second warning is the output of using the \textit{dsolve} command, which yielded no analytical solutions.
        \end{framed}
% 4f
        \item \hfill
        \begin{framed}
            \noindent My MATLAB script for problem 4 parts (f) and (g) is included at the end of this problem. This code implements seven numerical methods (first, second, and third order Taylor series approximations, modified Euler, midpoint Runge-Kutta (RK-2), RK-4, and RK-5) over various step sizes $dt = 100,10,1,0.1,0.01$ in order to solve the nonlinear ODE. Their absolute and relative errors when compared to the \textit{ode45} solution are also computed and plotted.
        \end{framed}
% 4g
        \item \hfill
        \begin{framed}
            \noindent For the seven numerical methods implemented in part (g), Figures \ref{fig:4g_errors_dt100} through \ref{fig:4g_errors_dt_point01} contain the plots of their absolute and relative errors when compared to the \textit{ode45} solution. \\ \\
            \noindent From these plots, we can immediately verify the existence of compounding errors (leading to global truncation errors) from the growth of the absolute error over time. Due to the singularity, we also have a huge spike in the relative error as it approaches that point. Within the Taylor series methods, I find that the absolute errors accumulate much more slowly as higher order TS methods are used. A similar pattern emerges within the RK methods as well, with RK-5 resulting in the least amount of error out of all the methods for every step size. More generally, it is also apparent that as the step size decreases, the numerical approximation gets closer to the actual solution, as shown by decreased errors for smaller step sizes.
            \begin{center}
                \includegraphics[width=1\linewidth]{4g_errors_dt100.eps}
                \captionof{figure}{Error plots for each numerical method to solve the nonlinear ODE with step size 100 seconds.}
                \label{fig:4g_errors_dt100}
            \end{center}
        \end{framed}
        \begin{framed}
            \begin{center}
                \includegraphics[width=1\linewidth]{4g_errors_dt10.eps}
                \captionof{figure}{Error plots for each numerical method to solve the nonlinear ODE with step size 10 seconds.}
                \label{fig:4g_errors_dt10}
            \end{center}
        \end{framed}
        \begin{framed}
            \begin{center}
                \includegraphics[width=1\linewidth]{4g_errors_dt1.eps}
                \captionof{figure}{Error plots for each numerical method to solve the nonlinear ODE with step size 1 second.}
                \label{fig:4g_errors_dt1}
            \end{center}
        \end{framed}
        \begin{framed}
            \begin{center}
                \includegraphics[width=1\linewidth]{4g_errors_dt_point1.eps}
                \captionof{figure}{Error plots for each numerical method to solve the nonlinear ODE with step size 0.1 seconds.}
                \label{fig:4g_errors_dt_point1}
            \end{center}
        \end{framed}
        \begin{framed}
            \begin{center}
                \includegraphics[width=1\linewidth]{4g_errors_dt_point01.eps}
                \captionof{figure}{Error plots for each numerical method to solve the nonlinear ODE with step size 0.01 seconds.}
                \label{fig:4g_errors_dt_point01}
            \end{center}
        \end{framed}
\newpage
% 4h
        \item \hfill
        \begin{framed}
            \noindent i) At the end of this problem is my MATLAB script that implements the explicit and implicit 4-step methods for step sizes $dt = 1000,100,10,1,0.1$, with the second, third, and fourth initial conditions approximated by the RK-4 method. For the implicit method, Newton's method was utilized to solve for the $i+1$ value of $h(t)$.
        \end{framed}
        \begin{framed}
            \noindent Figures \ref{fig:4hi_rk4_1000} through \ref{fig:4hi_rk4_point1} are the resulting error plots for the explicit and implicit methods that use RK-4, which is followed by a log-log plot of the computational times versus step size in Figure \ref{fig:4hi_rk4_comp_time}. From the errors plots, I notice that the implicit method performed about the same or slightly better than the explicit method. However, it doesn't seem like the step sizes had any significant effect on the error performance of the methods: the absolute errors all peak at around the same value. Looking at the computation time versus step size plot, it is clear that the implicit methods take much more time than the explicit methods, due to the addition of an iterative method to solve a nonlinear equation for the solution at the next time step (e.g. Newton's method).
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_i_rk4_dt1000.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using RK-4 with step size 1000 seconds.}
                \label{fig:4hi_rk4_1000}
            \end{center}
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_i_rk4_dt100.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using RK-4 with step size 100 seconds.}
                \label{fig:4hi_rk4_100}
            \end{center}
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_i_rk4_dt10.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using RK-4 with step size 10 seconds.}
                \label{fig:4hi_rk4_10}
            \end{center}
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_i_rk4_dt1.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using RK-4 with step size 1 seconds.}
                \label{fig:4hi_rk4_1}
            \end{center}
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_i_rk4_dt_point1.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using RK-4 with step size 0.1 seconds.}
                \label{fig:4hi_rk4_point1}
            \end{center}
            \begin{center}
                \includegraphics[width=0.75\linewidth]{4h_i_rk4_computional_times_vs_step_size.eps}
                \captionof{figure}{Log-log plot of the computational times vs step sizes for the explicit and implicit methods with RK-4 and Newton's method.}
                \label{fig:4hi_rk4_comp_time}
            \end{center}
        \end{framed}
        \begin{framed}
            \noindent Similarly, below are Figures \ref{fig:4hi_mod_euler_1000} through \ref{fig:4hi_mod_euler_point1}, which show the resulting error plots for the explicit and implicit methods but using the modified Euler methods instead of RK-4 to find the first four initial conditions. Just like before, the implicit method performed similarly to the explicit method, and step sizes had little effect on the error performance of the methods. \\ \\
            \noindent A plot of the computation time versus step size is shown in Figure \ref{fig:4hi_mod_euler_comp_time} when using modified Euler. It is crucial to note that the implicit method with modified Euler and step size $dt=0.1$ took over an hour to completely run, which is indicated in Figure \ref{fig:4hi_mod_euler_comp_time}. Overall, the computation times required for both methods when using modified Euler took significantly longer than those when using RK-4. However, the error performances for all four methods are very similar, so it would be advisable to use the explicit method with RK-4 (for any step size) solely due to computation times.
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_i_mod_euler_dt1000.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using modified Euler with step size 1000 seconds.}
                \label{fig:4hi_mod_euler_1000}
            \end{center}
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_i_mod_euler_dt100.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using modified Euler with step size 100 seconds.}
                \label{fig:4hi_mod_euler_100}
            \end{center}
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_i_mod_euler_dt10.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using modified Euler with step size 10 seconds.}
                \label{fig:4hi_mod_euler_10}
            \end{center}
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_i_mod_euler_dt1.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using modified Euler with step size 1 seconds.}
                \label{fig:4hi_mod_euler_1}
            \end{center}
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_i_mod_euler_dt_point1.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using modified Euler with step size 0.1 seconds.}
                \label{fig:4hi_mod_euler_point1}
            \end{center}
            \begin{center}
                \includegraphics[width=0.7\linewidth]{4h_i_mod_euler_computational_times_vs_step_size.eps}
                \captionof{figure}{Log-log plot of the computational times vs step sizes for the explicit and implicit methods with modified Euler and Newton's method.}
                \label{fig:4hi_mod_euler_comp_time}
            \end{center}
        \end{framed}
        \begin{framed}
            \noindent At the end of this problem is my MATLAB script for the explicit and implicit four-step methods that use modified Euler to compute the first four initial conditions.
        \end{framed}
        \begin{framed}
            \noindent ii) Replacing Newton's method in the previous codes for the implicit method with a fixed point iteration instead yields the MATLAB scripts labeled at the end of this problem.\\ \\
            \noindent Figures \ref{fig:4hii_rk4_1000} through \ref{fig:4hii_rk4_point1} are the error plots for the explicit and implicit methods when using RK-4 and a fixed point iteration (instead of Newton's method). Similarly, Figures \ref{fig:4hii_mod_euler_1000} through \ref{fig:4hii_mod_euler_point1} show those but when using modified Euler. Comparing the errors from these figures to those when Newton's method was used, I see no notable differences. However, when comparing the computational times from Figure \ref{fig:4hii_rk4_comp_time} and Figure \ref{fig:4hii_mod_euler_comp_time} to Figure \ref{fig:4hi_rk4_comp_time} and Figure \ref{fig:4hi_mod_euler_comp_time}, I notice that the run time is drastically lower for the implicit method with fixed point iteration. This is most likely due to Newton's method diverging most of the time when solving the nonlinear equation for the implicit method, while the fixed point iteration mostly converges.
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_ii_rk4_dt1000.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using RK-4 with step size 1000 seconds.}
                \label{fig:4hii_rk4_1000}
            \end{center}
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_ii_rk4_dt100.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using RK-4 with step size 100 seconds.}
                \label{fig:4hii_rk4_100}
            \end{center}
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_i_mod_euler_dt1.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using RK-4 with step size 10 seconds.}
                \label{fig:4hii_rk4_10}
            \end{center}
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_ii_rk4_dt1.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using RK-4 with step size 1 seconds.}
                \label{fig:4hii_rk4_1}
            \end{center}
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_ii_rk4_dt_point1.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using RK-4 with step size 0.1 seconds.}
                \label{fig:4hii_rk4_point1}
            \end{center}
            \begin{center}
                \includegraphics[width=0.75\linewidth]{4h_ii_rk4_comp_time_vs_step_size.eps}
                \captionof{figure}{Log-log plot of the computational times vs step sizes for the explicit and implicit methods with RK-4 and fixed-point iteration.}
                \label{fig:4hii_rk4_comp_time}
            \end{center}
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_ii_mod_euler_dt1000.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using modified Euler with step size 1000 seconds.}
                \label{fig:4hii_mod_euler_1000}
            \end{center}
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_ii_mod_euler_dt100.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using modified Euler with step size 100 seconds.}
                \label{fig:4hii_mod_euler_100}
            \end{center}
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_i_mod_euler_dt1.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using modified Euler with step size 10 seconds.}
                \label{fig:4hii_mod_euler_10}
            \end{center}
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_ii_mod_euler_dt1.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using modified Euler with step size 1 seconds.}
                \label{fig:4hii_mod_euler_1}
            \end{center}
            \begin{center}
                \includegraphics[width=0.8\linewidth]{4h_ii_mod_euler_dt_point1.eps}
                \captionof{figure}{Error plots for explicit and implicit methods using modified Euler with step size 0.1 seconds.}
                \label{fig:4hii_mod_euler_point1}
            \end{center}
            \begin{center}
                \includegraphics[width=0.75\linewidth]{4h_ii_mod_euler_comp_time_vs_step_size.eps}
                \captionof{figure}{Log-log plot of the computational times vs step sizes for the explicit and implicit methods with modified Euler and fixed point iteration.}
                \label{fig:4hii_mod_euler_comp_time}
            \end{center}
        \end{framed}
        \begin{framed}
            \noindent iii) Using the predictor-corrector method results in the error plots in Figure \ref{fig:4hiii_rk4_errors} and \ref{fig:4hiii_mod_euler_errors}, when using RK-4 and modified Euler, respectively. The plots of computational times versus step size are shown in Figures \ref{fig:4hiii_rk4_comp_times} and \ref{fig:4hiii_mod_euler_comp_times}. From these plots, I see that the error performances are still similar to the previous methods, but the computation times are similar to the explicit methods. Comparing to the implicit methods with Newton's method, the computation times for the predictor-corrector methods are about 10 times faster at worst, and about $10^3$ times faster at best; this is HUGE difference that makes this method much more efficient, given that it results in similar errors.
            \begin{center}
                \includegraphics[width=1\linewidth]{4h_iii_rk4_errors.eps}
                \captionof{figure}{Error plots for the predictor-corrector method using RK-4.}
                \label{fig:4hiii_rk4_errors}
            \end{center}
            \begin{center}
                \includegraphics[width=0.75\linewidth]{4h_iii_rk4_comp_time_vs_step_size.eps}
                \captionof{figure}{Log-log plot of the computational times vs step sizes for the explicit and implicit methods with RK-4.}
                \label{fig:4hiii_rk4_comp_times}
            \end{center}
            \begin{center}
                \includegraphics[width=1\linewidth]{4h_iii_mod_euler_errors.eps}
                \captionof{figure}{Error plots for the predictor-corrector method using modified Euler}
                \label{fig:4hiii_mod_euler_errors}
            \end{center}
            \begin{center}
                \includegraphics[width=0.75\linewidth]{4h_iii_mod_euler_comp_time_vs_step_size.eps}
                \captionof{figure}{Log-log plot of the computational times vs step sizes for the explicit and implicit methods with modified Euler.}
                \label{fig:4hiii_mod_euler_comp_times}
            \end{center}
        \end{framed}
    \end{enumerate}

% MATLAB CODES
\begin{framed}
    \noindent Below are my MATLAB scripts and their outputs (if any) for problem 4:
\end{framed}
        \begin{center}
            \lstinputlisting[frame=single, numbers=left, style=Matlab-Pyglike]{problem_4ae.m}
        \end{center}
        \begin{center}
            \lstinputlisting[frame=single, numbers=left, style=Matlab-Pyglike]{problem_4ae_output.m}
        \end{center}
        
        \begin{center}
            \lstinputlisting[frame=single, numbers=left, style=Matlab-Pyglike]{problem_4fg.m}
        \end{center}

        \begin{center}
            \lstinputlisting[frame=single, numbers=left, style=Matlab-Pyglike]{problem_4h_i_rk4.m}
        \end{center}

        \begin{center}
            \lstinputlisting[frame=single, numbers=left, style=Matlab-Pyglike]{problem_4h_i_modified_euler.m}
        \end{center}

        \begin{center}
            \lstinputlisting[frame=single, numbers=left, style=Matlab-Pyglike]{problem_4h_ii_rk4.m}
        \end{center}
        \begin{center}
            \lstinputlisting[frame=single, numbers=left, style=Matlab-Pyglike]{problem_4h_ii_modified_euler.m}
        \end{center}
    
% Image
%            \begin{center}
%                \includegraphics[width=0.75\linewidth]{}
%                \captionof{figure}{}
%                \label{fig:}
%            \end{center}
% Table
%           \begin{center}
%               \begin{tabular}{|c|c|c|c|}
%                   \hline
%                   Col1 & COl2 & Col3 & Col4 \\
%                   \hline
%                   7    & 6     & 5     & 8     \\
%                   7    & 6     & 6     & 7     \\
%                   7    & 6     & 7     & 6     \\
%                   7    & 7     & 5     & 9     \\
%                   \hline
%               \end{tabular}
%               \captionof{table}{}
%               \label{tab:}
%           \end{center}
% Colored MATLAB code:
%    \begin{center}
%        \lstinputlisting[frame=single, numbers=left, style=Matlab-Pyglike]{Q3.m}
%    \end{center}

\end{enumerate}

\end{document}