InitHEADmain

35 files changed, 3143 insertions, 0 deletions

diff --git a/Homework/math4610/homeworks/a.out b/Homework/math4610/homeworks/a.out
new file mode 100644
index 0000000..410a7d5
--- /dev/null
+++ b/Homework/math4610/homeworks/a.out
diff --git a/Homework/math4610/homeworks/hw-2.org b/Homework/math4610/homeworks/hw-2.org
new file mode 100644
index 0000000..f696ffb
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-2.org
@@ -0,0 +1,182 @@
+#+TITLE: HW 02
+#+AUTHOR: Elizabeth Hunt
+#+STARTUP: entitiespretty fold inlineimages
+#+LATEX_HEADER: \notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+#+LATEX: \setlength\parindent{0pt}
+#+OPTIONS: toc:nil
+
+* Question One
+Computing $\epsilon_{\text{mac}}$ for single precision numbers
+
+#+BEGIN_SRC lisp :session t :results table 
+  (load "../lizfcm.asd")
+  (ql:quickload :lizfcm)
+
+  (let ((domain-values (lizfcm.approx:compute-maceps (lambda (x) x)
+                                                     1.0
+                                                     1.0)))
+    (lizfcm.utils:table (:headers '("a" "h" "err")
+                         :domain-order (a h err)
+                         :domain-values domain-values)))
+#+END_SRC
+
+(with many rows truncated)
+
+|   a |             h |           err |
+| 1.0 |           1.0 |           1.0 |
+| 1.0 |           0.5 |           0.5 |
+| 1.0 |          0.25 |          0.25 |
+| 1.0 |         0.125 |         0.125 |
+| 1.0 |        0.0625 |        0.0625 |
+| 1.0 |       0.03125 |       0.03125 |
+| 1.0 | 1.9073486e-06 | 1.9073486e-06 |
+| 1.0 |  9.536743e-07 |  9.536743e-07 |
+| 1.0 | 4.7683716e-07 | 4.7683716e-07 |
+| 1.0 | 2.3841858e-07 | 2.3841858e-07 |
+| 1.0 | 1.1920929e-07 | 1.1920929e-07 |
+
+$\epsilon_{\text{mac single precision}}$ \approx 1.192(10^-7)
+
+* Question Two
+Computing $\epsilon_{\text{mac}}$ for double precision numbers:
+
+#+BEGIN_SRC lisp :session t :results table
+  (let ((domain-values (lizfcm.approx:compute-maceps (lambda (x) x)
+                                                     1.0d0
+                                                     1.0d0)))
+    (lizfcm.utils:table (:headers '("a" "h" "err")
+                         :domain-order (a h err)
+                         :domain-values domain-values)))
+#+END_SRC
+
+(with many rows truncated)
+|     a |                      h |                    err |
+| 1.0d0 |                  1.0d0 |                  1.0d0 |
+| 1.0d0 |                  0.5d0 |                  0.5d0 |
+| 1.0d0 |                 0.25d0 |                 0.25d0 |
+| 1.0d0 |                0.125d0 |                0.125d0 |
+| 1.0d0 |               0.0625d0 |               0.0625d0 |
+| 1.0d0 |              0.03125d0 |              0.03125d0 |
+| 1.0d0 |             0.015625d0 |             0.015625d0 |
+| 1.0d0 |            0.0078125d0 |            0.0078125d0 |
+| 1.0d0 |           0.00390625d0 |           0.00390625d0 |
+| 1.0d0 |          0.001953125d0 |          0.001953125d0 |
+| 1.0d0 |  7.105427357601002d-15 |  7.105427357601002d-15 |
+| 1.0d0 |  3.552713678800501d-15 |  3.552713678800501d-15 |
+| 1.0d0 | 1.7763568394002505d-15 | 1.7763568394002505d-15 |
+| 1.0d0 |  8.881784197001252d-16 |  8.881784197001252d-16 |
+| 1.0d0 |  4.440892098500626d-16 |  4.440892098500626d-16 |
+| 1.0d0 |  2.220446049250313d-16 |  2.220446049250313d-16 |
+
+Thus, $\epsilon_{\text{mac double precision}}$ \approx 2.220 \cdot 10^{-16}
+
+* Question Three - |v|_2
+#+BEGIN_SRC lisp :session t
+  (let ((vs '((1 1) (2 3) (4 5) (-1 2)))
+        (2-norm (lizfcm.vector:p-norm 2)))
+    (lizfcm.utils:table (:headers '("x" "y" "2norm")
+                         :domain-order (x y)
+                         :domain-values vs)
+      (funcall 2-norm (list x y))))
+#+END_SRC
+
+
+|  x | y |     2norm |
+|  1 | 1 | 1.4142135 |
+|  2 | 3 | 3.6055512 |
+|  4 | 5 | 6.4031243 |
+| -1 | 2 |  2.236068 |
+
+* Question Four - |v|_1
+#+BEGIN_SRC lisp :session t
+  (let ((vs '((1 1) (2 3) (4 5) (-1 2)))
+        (1-norm (lizfcm.vector:p-norm 1)))
+    (lizfcm.utils:table (:headers '("x" "y" "1norm")
+                         :domain-order (x y)
+                         :domain-values vs)
+      (funcall 1-norm (list x y))))
+#+END_SRC
+
+
+|  x | y | 1norm |
+|  1 | 1 |     2 |
+|  2 | 3 |     5 |
+|  4 | 5 |     9 |
+| -1 | 2 |     3 |
+
+* Question Five - |v|_{\infty}
+#+BEGIN_SRC lisp :session t
+  (let ((vs '((1 1) (2 3) (4 5) (-1 2))))
+    (lizfcm.utils:table (:headers '("x" "y" "max-norm")
+                         :domain-order (x y)
+                         :domain-values vs)
+      (lizfcm.vector:max-norm (list x y))))
+#+END_SRC
+
+
+|  x | y | infty-norm |
+|  1 | 1 |          1 |
+|  2 | 3 |          3 |
+|  4 | 5 |          5 |
+| -1 | 2 |          2 |
+
+* Question Six - ||v - u|| via |v|_{2}
+#+BEGIN_SRC lisp :session t
+  (let ((vs '((1 1) (2 3) (4 5) (-1 2)))
+        (vs2 '((7 9) (2 2) (8 -1) (4 4)))
+        (2-norm (lizfcm.vector:p-norm 2)))
+    (lizfcm.utils:table (:headers '("v1" "v2" "2-norm-d")
+                         :domain-order (v1 v2)
+                         :domain-values (mapcar (lambda (v1 v2)
+                                                  (list v1 v2))
+                                                vs
+                                                vs2))
+      (lizfcm.vector:distance v1 v2 2-norm)))
+#+END_SRC
+
+
+| v1     | v2     |    2-norm |
+| (1 1)  | (7 9)  |      10.0 |
+| (2 3)  | (2 2)  |       1.0 |
+| (4 5)  | (8 -1) | 7.2111025 |
+| (-1 2) | (4 4)  | 5.3851647 |
+
+* Question Seven - ||v - u|| via |v|_{1}
+#+BEGIN_SRC lisp :session t
+  (let ((vs '((1 1) (2 3) (4 5) (-1 2)))
+        (vs2 '((7 9) (2 2) (8 -1) (4 4)))
+        (1-norm (lizfcm.vector:p-norm 1)))
+    (lizfcm.utils:table (:headers '("v1" "v2" "1-norm-d")
+                         :domain-order (v1 v2)
+                         :domain-values (mapcar (lambda (v1 v2)
+                                                  (list v1 v2))
+                                                vs
+                                                vs2))
+      (lizfcm.vector:distance v1 v2 1-norm)))
+#+END_SRC
+
+
+| v1     | v2     | 1-norm-d |
+| (1 1)  | (7 9)  |       14 |
+| (2 3)  | (2 2)  |        1 |
+| (4 5)  | (8 -1) |       10 |
+| (-1 2) | (4 4)  |        7 |
+
+* Question Eight - ||v - u|| via |v|_{\infty}
+#+BEGIN_SRC lisp :session t
+  (let ((vs '((1 1) (2 3) (4 5) (-1 2)))
+        (vs2 '((7 9) (2 2) (8 -1) (4 4))))
+    (lizfcm.utils:table (:headers '("v1" "v2" "max-norm-d")
+                         :domain-order (v1 v2)
+                         :domain-values (mapcar (lambda (v1 v2)
+                                                  (list v1 v2))
+                                                vs
+                                                vs2))
+      (lizfcm.vector:distance v1 v2 'lizfcm.vector:max-norm)))
+#+END_SRC
+
+| v1     | v2     | max-norm-d |
+| (1 1)  | (7 9)  |         -6 |
+| (2 3)  | (2 2)  |          1 |
+| (4 5)  | (8 -1) |          6 |
+| (-1 2) | (4 4)  |         -2 |
diff --git a/Homework/math4610/homeworks/hw-2.pdf b/Homework/math4610/homeworks/hw-2.pdf
new file mode 100644
index 0000000..2dc4d28
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-2.pdf
diff --git a/Homework/math4610/homeworks/hw-2.tex b/Homework/math4610/homeworks/hw-2.tex
new file mode 100644
index 0000000..da8d3f5
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-2.tex
@@ -0,0 +1,246 @@
+% Created 2023-10-07 Sat 14:51
+% Intended LaTeX compiler: pdflatex
+\documentclass[11pt]{article}
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+\usepackage{graphicx}
+\usepackage{longtable}
+\usepackage{wrapfig}
+\usepackage{rotating}
+\usepackage[normalem]{ulem}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{capt-of}
+\usepackage{hyperref}
+\notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+\author{Elizabeth Hunt}
+\date{\today}
+\title{HW 02}
+\hypersetup{
+ pdfauthor={Elizabeth Hunt},
+ pdftitle={HW 02},
+ pdfkeywords={},
+ pdfsubject={},
+ pdfcreator={Emacs 28.2 (Org mode 9.7-pre)}, 
+ pdflang={English}}
+\begin{document}
+
+\maketitle
+\setlength\parindent{0pt}
+
+\section{Question One}
+\label{sec:org58b9af4}
+Computing \(\epsilon_{\text{mac}}\) for single precision numbers
+
+\begin{verbatim}
+(load "../lizfcm.asd")
+(ql:quickload :lizfcm)
+
+(let ((domain-values (lizfcm.approx:compute-maceps (lambda (x) x)
+                                                   1.0
+                                                   1.0)))
+  (lizfcm.utils:table (:headers '("a" "h" "err")
+                       :domain-order (a h err)
+                       :domain-values domain-values)))
+\end{verbatim}
+
+(with many rows truncated)
+
+\begin{center}
+\begin{tabular}{rrr}
+a & h & err\\[0pt]
+1.0 & 1.0 & 1.0\\[0pt]
+1.0 & 0.5 & 0.5\\[0pt]
+1.0 & 0.25 & 0.25\\[0pt]
+1.0 & 0.125 & 0.125\\[0pt]
+1.0 & 0.0625 & 0.0625\\[0pt]
+1.0 & 0.03125 & 0.03125\\[0pt]
+1.0 & 1.9073486e-06 & 1.9073486e-06\\[0pt]
+1.0 & 9.536743e-07 & 9.536743e-07\\[0pt]
+1.0 & 4.7683716e-07 & 4.7683716e-07\\[0pt]
+1.0 & 2.3841858e-07 & 2.3841858e-07\\[0pt]
+1.0 & 1.1920929e-07 & 1.1920929e-07\\[0pt]
+\end{tabular}
+\end{center}
+
+\(\epsilon_{\text{mac single precision}}\) \(\approx\) 1.192(10\textsuperscript{-7})
+
+\section{Question Two}
+\label{sec:org27557b4}
+Computing \(\epsilon_{\text{mac}}\) for double precision numbers:
+
+\begin{verbatim}
+(let ((domain-values (lizfcm.approx:compute-maceps (lambda (x) x)
+                                                   1.0d0
+                                                   1.0d0)))
+  (lizfcm.utils:table (:headers '("a" "h" "err")
+                       :domain-order (a h err)
+                       :domain-values domain-values)))
+\end{verbatim}
+
+(with many rows truncated)
+\begin{center}
+\begin{tabular}{rrr}
+a & h & err\\[0pt]
+1.0d0 & 1.0d0 & 1.0d0\\[0pt]
+1.0d0 & 0.5d0 & 0.5d0\\[0pt]
+1.0d0 & 0.25d0 & 0.25d0\\[0pt]
+1.0d0 & 0.125d0 & 0.125d0\\[0pt]
+1.0d0 & 0.0625d0 & 0.0625d0\\[0pt]
+1.0d0 & 0.03125d0 & 0.03125d0\\[0pt]
+1.0d0 & 0.015625d0 & 0.015625d0\\[0pt]
+1.0d0 & 0.0078125d0 & 0.0078125d0\\[0pt]
+1.0d0 & 0.00390625d0 & 0.00390625d0\\[0pt]
+1.0d0 & 0.001953125d0 & 0.001953125d0\\[0pt]
+1.0d0 & 7.105427357601002d-15 & 7.105427357601002d-15\\[0pt]
+1.0d0 & 3.552713678800501d-15 & 3.552713678800501d-15\\[0pt]
+1.0d0 & 1.7763568394002505d-15 & 1.7763568394002505d-15\\[0pt]
+1.0d0 & 8.881784197001252d-16 & 8.881784197001252d-16\\[0pt]
+1.0d0 & 4.440892098500626d-16 & 4.440892098500626d-16\\[0pt]
+1.0d0 & 2.220446049250313d-16 & 2.220446049250313d-16\\[0pt]
+\end{tabular}
+\end{center}
+
+Thus, \(\epsilon_{\text{mac double precision}}\) \(\approx\) 2.220 \(\cdot\) 10\textsuperscript{-16}
+
+\section{Question Three - |v|\textsubscript{2}}
+\label{sec:org59c6c10}
+\begin{verbatim}
+(let ((vs '((1 1) (2 3) (4 5) (-1 2)))
+      (2-norm (lizfcm.vector:p-norm 2)))
+  (lizfcm.utils:table (:headers '("x" "y" "2norm")
+                       :domain-order (x y)
+                       :domain-values vs)
+    (funcall 2-norm (list x y))))
+\end{verbatim}
+
+
+\begin{center}
+\begin{tabular}{rrr}
+x & y & 2norm\\[0pt]
+1 & 1 & 1.4142135\\[0pt]
+2 & 3 & 3.6055512\\[0pt]
+4 & 5 & 6.4031243\\[0pt]
+-1 & 2 & 2.236068\\[0pt]
+\end{tabular}
+\end{center}
+
+\section{Question Four - |v|\textsubscript{1}}
+\label{sec:org2b67b3e}
+\begin{verbatim}
+(let ((vs '((1 1) (2 3) (4 5) (-1 2)))
+      (1-norm (lizfcm.vector:p-norm 1)))
+  (lizfcm.utils:table (:headers '("x" "y" "1norm")
+                       :domain-order (x y)
+                       :domain-values vs)
+    (funcall 1-norm (list x y))))
+\end{verbatim}
+
+
+\begin{center}
+\begin{tabular}{rrr}
+x & y & 1norm\\[0pt]
+1 & 1 & 2\\[0pt]
+2 & 3 & 5\\[0pt]
+4 & 5 & 9\\[0pt]
+-1 & 2 & 3\\[0pt]
+\end{tabular}
+\end{center}
+
+\section{Question Five - |v|\textsubscript{\(\infty\)}}
+\label{sec:org922206e}
+\begin{verbatim}
+(let ((vs '((1 1) (2 3) (4 5) (-1 2))))
+  (lizfcm.utils:table (:headers '("x" "y" "max-norm")
+                       :domain-order (x y)
+                       :domain-values vs)
+    (lizfcm.vector:max-norm (list x y))))
+\end{verbatim}
+
+
+\begin{center}
+\begin{tabular}{rrr}
+x & y & infty-norm\\[0pt]
+1 & 1 & 1\\[0pt]
+2 & 3 & 3\\[0pt]
+4 & 5 & 5\\[0pt]
+-1 & 2 & 2\\[0pt]
+\end{tabular}
+\end{center}
+
+\section{Question Six - ||v - u|| via |v|\textsubscript{2}}
+\label{sec:org29ec18f}
+\begin{verbatim}
+(let ((vs '((1 1) (2 3) (4 5) (-1 2)))
+      (vs2 '((7 9) (2 2) (8 -1) (4 4)))
+      (2-norm (lizfcm.vector:p-norm 2)))
+  (lizfcm.utils:table (:headers '("v1" "v2" "2-norm-d")
+                       :domain-order (v1 v2)
+                       :domain-values (mapcar (lambda (v1 v2)
+                                                (list v1 v2))
+                                              vs
+                                              vs2))
+    (lizfcm.vector:distance v1 v2 2-norm)))
+\end{verbatim}
+
+
+\begin{center}
+\begin{tabular}{llr}
+v1 & v2 & 2-norm\\[0pt]
+(1 1) & (7 9) & 10.0\\[0pt]
+(2 3) & (2 2) & 1.0\\[0pt]
+(4 5) & (8 -1) & 7.2111025\\[0pt]
+(-1 2) & (4 4) & 5.3851647\\[0pt]
+\end{tabular}
+\end{center}
+
+\section{Question Seven - ||v - u|| via |v|\textsubscript{1}}
+\label{sec:org7a87810}
+\begin{verbatim}
+(let ((vs '((1 1) (2 3) (4 5) (-1 2)))
+      (vs2 '((7 9) (2 2) (8 -1) (4 4)))
+      (1-norm (lizfcm.vector:p-norm 1)))
+  (lizfcm.utils:table (:headers '("v1" "v2" "1-norm-d")
+                       :domain-order (v1 v2)
+                       :domain-values (mapcar (lambda (v1 v2)
+                                                (list v1 v2))
+                                              vs
+                                              vs2))
+    (lizfcm.vector:distance v1 v2 1-norm)))
+\end{verbatim}
+
+
+\begin{center}
+\begin{tabular}{llr}
+v1 & v2 & 1-norm-d\\[0pt]
+(1 1) & (7 9) & 14\\[0pt]
+(2 3) & (2 2) & 1\\[0pt]
+(4 5) & (8 -1) & 10\\[0pt]
+(-1 2) & (4 4) & 7\\[0pt]
+\end{tabular}
+\end{center}
+
+\section{Question Eight - ||v - u|| via |v|\textsubscript{\(\infty\)}}
+\label{sec:org0f3b64f}
+\begin{verbatim}
+(let ((vs '((1 1) (2 3) (4 5) (-1 2)))
+      (vs2 '((7 9) (2 2) (8 -1) (4 4))))
+  (lizfcm.utils:table (:headers '("v1" "v2" "max-norm-d")
+                       :domain-order (v1 v2)
+                       :domain-values (mapcar (lambda (v1 v2)
+                                                (list v1 v2))
+                                              vs
+                                              vs2))
+    (lizfcm.vector:distance v1 v2 'lizfcm.vector:max-norm)))
+\end{verbatim}
+
+\begin{center}
+\begin{tabular}{llr}
+v1 & v2 & max-norm-d\\[0pt]
+(1 1) & (7 9) & -6\\[0pt]
+(2 3) & (2 2) & 1\\[0pt]
+(4 5) & (8 -1) & 6\\[0pt]
+(-1 2) & (4 4) & -2\\[0pt]
+\end{tabular}
+\end{center}
+\end{document}
\ No newline at end of file
diff --git a/Homework/math4610/homeworks/hw-3.org b/Homework/math4610/homeworks/hw-3.org
new file mode 100644
index 0000000..8d7edcc
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-3.org
@@ -0,0 +1,242 @@
+#+TITLE: HW 03
+#+AUTHOR: Elizabeth Hunt
+#+STARTUP: entitiespretty fold inlineimages
+#+LATEX_HEADER: \notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+#+LATEX: \setlength\parindent{0pt}
+#+OPTIONS: toc:nil
+
+* Question One
+** Three Terms
+\begin{align*}
+Si_3(x) &= \int_0^x \frac{s - \frac{s^3}{3!} + \frac{s^5}{5!}}{s} dx \\
+&= x - \frac{x^3}{(3!)(3)} + \frac{x^5}{(5!)(5)}
+\end{align*}
+** Five Terms
+\begin{align*}
+Si_3(x) &= \int_0^x \frac{s - \frac{s^3}{3!} + \frac{s^5}{5!} - \frac{s^7}{7!} + \frac{s^9}{9!}}{s} dx \\
+&= x - \frac{x^3}{(3!)(3)} + \frac{x^5}{(5!)(5)} - \frac{x^7}{(7!)(7)} + \frac{s^9}{(9!)(9)}
+\end{align*}
+** Ten Terms
+\begin{align*}
+Si_{10}(x) &= \int_0^x \frac{s - \frac{s^3}{3!} + \frac{s^5}{5!} - \frac{s^7}{7!} + \frac{s^9}{9!} - \frac{s^{11}}{11!} + \frac{s^{13}}{13!} - \frac{s^{15}}{15!} + \frac{s^{17}}{17!} - \frac{s^{19}}{19!}}{s} ds \\
+&= x - \frac{x^3}{(3!)(3)} + \frac{x^5}{(5!)(5)} - \frac{x^7}{(7!)(7)} + \frac{s^9}{(9!)(9)} - \frac{s^{11}}{(11!)(11)} + \frac{s^{13}}{(13!)(13)} - \frac{s^{15}}{(15!)(15)} \\
+&+ \frac{s^{17}}{(17!)(17)} - \frac{s^{19}}{(19!)(19)}
+\end{align*}
+* Question Three
+For the second term in the difference quotient, we can expand the taylor series centered at x=a:
+
+\begin{equation*}
+f(x) = f(a) + f'(a)(x-a) + \frac{f''(a)}{2}(x-a)^2 + \cdots \\
+\end{equation*}
+
+Which we substitute into the difference quotient:
+
+\begin{equation*}
+\frac{f(a) - f(a - h)}{h} = \frac{f(a) - (f(a) + f'(a)(x-a) + \frac{f''(a)}{2}(x-a)^2 + \cdots)}{h}
+\end{equation*}
+
+And subs. $x=a-h$:
+
+\begin{align*}
+\frac{f(a) - (f(a) + f'(a)(x-a) + \frac{f''(a)}{2}(x-a)^2 + \cdots)}{h} &= -f'(a)(-1) + -\frac{1}{2}f''(a)h \\
+&= f'(a) - \frac{1}{2}f''(a)h + \cdots \\
+\end{align*}
+
+Which we now plug into the initial $e_{\text{abs}}$:
+
+\begin{align*}
+e_{\text{abs}} &= |f'(a) - \frac{f(a) - f(a - h)}{h}| \\
+&= |f'(a) - (f'(a) +  -\frac{f''(a)}{2}h + \cdots)| \\
+&= |- \frac{1}{2}f''(a)h + \cdots | \\
+\end{align*}
+
+With the Taylor Remainder theorem we can absorb the series following the second term:
+
+\begin{equation*}
+e_{\text{abs}} = |- \frac{1}{2}f''(a)h + \cdots | = |\frac{1}{2}f''(\xi)h| \leq Ch
+\end{equation*}
+
+Thus our error is bounded linearly with $h$.
+
+* Question Four
+For the first term in the difference quotient we know, from the given notes,
+
+\begin{equation*}
+f(a+h) = f(a) + f'(a)h + \frac{1}{2}f''(a)h^2 + \frac{1}{6}f'''(a)(h^3)
+\end{equation*}
+
+And from some of the work in Question Three,
+
+\begin{equation*}
+f(a - h) = f(a) + f'(a)(-h) + \frac{1}{2}f''(a)(-h)^2 + \frac{1}{6}f'''(a)(-h^3)
+\end{equation*}
+
+We can substitute immediately into $e_{\text{abs}} = |f'(a) - (\frac{f(a+h) - f(a-h)}{2h})|$:
+
+\begin{align*}
+e_{\text{abs}} &= |f'(a) - \frac{1}{2h}((f(a) + f'(a)h + \frac{1}{2}f''(a)h^2 + \cdots) - (f(a) - f'(a)h + \frac{1}{2}f''(a)h^2 + \cdots))| \\
+&= |f'(a) - \frac{1}{2h}(2f'(a)h + \frac{1}{6}f'''(a)h^3 + \cdots)| \\
+&= |f'(a) - f'(a) - \frac{1}{12}f'''(a)h^2 + \cdots| \\
+&= |-\frac{1}{12}f'''(a)h^2 + \cdots|
+\end{align*}
+
+Finally, with the Taylor Remainder theorem we can absorb the series following the third term:
+
+\begin{equation*}
+e_{\text{abs}} = |-\frac{1}{12}f'''(\xi)h^2| = |\frac{1}{12}f'''(\xi)h^2| \leq Ch^2
+\end{equation*}
+
+Meaning that as $h$ scales linearly, our error is bounded by $h^2$ as opposed to linearly as in Question Three.
+
+* Question Six
+** A
+#+BEGIN_SRC lisp
+  (load "../lizfcm.asd")
+  (ql:quickload :lizfcm)
+
+  (defun f (x)
+    (/ (- x 1) (+ x 1)))
+
+  (defun fprime (x)
+    (/ 2 (expt (+ x 1) 2)))
+
+  (let ((domain-values (loop for a from 0 to 2
+                             append 
+                             (loop for i from 0 to 9
+                                   for h = (/ 1.0 (expt 2 i))
+                                   collect (list a h)))))
+    (lizfcm.utils:table (:headers '("a" "h" "f'" "\\approx f'" "e_{\\text{abs}}")
+                         :domain-order (a h)
+                         :domain-values domain-values)
+      (fprime a)
+      (lizfcm.approx:fwd-derivative-at 'f a h)
+      (abs (- (fprime a)
+              (lizfcm.approx:fwd-derivative-at 'f a h)))))
+#+END_SRC
+
+#+RESULTS:
+| a |           h | f'  |       \approx f' |   e_{\text{abs}} |
+| 0 |         1.0 | 2   |        1.0 |           1.0 |
+| 0 |         0.5 | 2   |  1.3333333 |    0.66666675 |
+| 0 |        0.25 | 2   |  1.5999999 |     0.4000001 |
+| 0 |       0.125 | 2   |  1.7777777 |    0.22222233 |
+| 0 |      0.0625 | 2   |  1.8823528 |    0.11764717 |
+| 0 |     0.03125 | 2   |   1.939394 |   0.060606003 |
+| 0 |    0.015625 | 2   |  1.9692307 |   0.030769348 |
+| 0 |   0.0078125 | 2   |  1.9844971 |    0.01550293 |
+| 0 |  0.00390625 | 2   |   1.992218 |  0.0077819824 |
+| 0 | 0.001953125 | 2   |  1.9960938 |    0.00390625 |
+| 1 |         1.0 | 1/2 | 0.33333334 |    0.16666666 |
+| 1 |         0.5 | 1/2 |        0.4 |   0.099999994 |
+| 1 |        0.25 | 1/2 | 0.44444445 |   0.055555552 |
+| 1 |       0.125 | 1/2 | 0.47058824 |   0.029411763 |
+| 1 |      0.0625 | 1/2 |  0.4848485 |   0.015151501 |
+| 1 |     0.03125 | 1/2 |  0.4923077 |  0.0076923072 |
+| 1 |    0.015625 | 1/2 | 0.49612403 |  0.0038759708 |
+| 1 |   0.0078125 | 1/2 | 0.49805447 |  0.0019455254 |
+| 1 |  0.00390625 | 1/2 | 0.49902534 | 0.00097465515 |
+| 1 | 0.001953125 | 1/2 |  0.4995122 |  0.0004878044 |
+| 2 |         1.0 | 2/9 | 0.16666666 |   0.055555567 |
+| 2 |         0.5 | 2/9 | 0.19047618 |   0.031746045 |
+| 2 |        0.25 | 2/9 |  0.2051282 |   0.017094031 |
+| 2 |       0.125 | 2/9 | 0.21333337 |   0.008888856 |
+| 2 |      0.0625 | 2/9 | 0.21768713 |   0.004535094 |
+| 2 |     0.03125 | 2/9 | 0.21993065 |   0.002291575 |
+| 2 |    0.015625 | 2/9 | 0.22106934 |  0.0011528879 |
+| 2 |   0.0078125 | 2/9 | 0.22164536 | 0.00057686865 |
+| 2 |  0.00390625 | 2/9 | 0.22193146 | 0.00029076636 |
+| 2 | 0.001953125 | 2/9 | 0.22207642 | 0.00014580786 |
+
+* Question Nine
+** C
+
+#+BEGIN_SRC lisp
+  (load "../lizfcm.asd")
+  (ql:quickload :lizfcm)
+
+  (defun factorial (n)
+    (if (= n 0)
+        1
+        (* n (factorial (- n 1)))))
+
+  (defun taylor-term (n x)
+    (/ (* (expt (- 1) n)
+          (expt x (+ (* 2 n) 1)))
+       (* (factorial n)
+          (+ (* 2 n) 1))))
+
+  (defun f (x &optional (max-iterations 30))
+    (let ((sum 0.0))
+      (dotimes (n max-iterations)
+        (setq sum (+ sum (taylor-term n x))))
+      (* sum (/ 2 (sqrt pi)))))
+
+  (defun fprime (x)
+    (* (/ 2 (sqrt pi)) (exp (- 0 (* x x)))))
+
+  (let ((domain-values (loop for a from 0 to 1
+                             append 
+                             (loop for i from 0 to 9
+                                   for h = (/ 1.0 (expt 2 i))
+                                   collect (list a h)))))
+    (lizfcm.utils:table (:headers '("a" "h" "f'" "\\approx f'" "e_{\\text{abs}}")
+                         :domain-order (a h)
+                         :domain-values domain-values)
+      (fprime a)
+      (lizfcm.approx:central-derivative-at 'f a h)
+      (abs (- (fprime a)
+              (lizfcm.approx:central-derivative-at 'f a h)))))
+#+END_SRC
+
+
+| a |           h |                    f' |                  \approx f' |             e_{\text{abs}} |
+| 0 |         1.0 |  1.1283791670955126d0 |  0.8427006725464232d0 |   0.28567849454908933d0 |
+| 0 |         0.5 |  1.1283791670955126d0 |  1.0409997446922075d0 |    0.0873794224033051d0 |
+| 0 |        0.25 |  1.1283791670955126d0 |  1.1053055663206806d0 |  0.023073600774832004d0 |
+| 0 |       0.125 |  1.1283791670955126d0 |   1.122529655394656d0 |  0.005849511700856569d0 |
+| 0 |      0.0625 |  1.1283791670955126d0 |  1.1269116944798618d0 | 0.0014674726156507223d0 |
+| 0 |     0.03125 |  1.1283791670955126d0 |  1.1280120131008824d0 |   3.6715399463016496d-4 |
+| 0 |    0.015625 |  1.1283791670955126d0 |  1.1282873617826952d0 |    9.180531281738347d-5 |
+| 0 |   0.0078125 |  1.1283791670955126d0 |   1.128356232581468d0 |    2.293451404455915d-5 |
+| 0 |  0.00390625 |  1.1283791670955126d0 |  1.1283734502811613d0 |     5.71681435124205d-6 |
+| 0 | 0.001953125 |  1.1283791670955126d0 |  1.1283777547060847d0 |   1.4123894278572635d-6 |
+| 1 |         1.0 | 0.41510750774498784d0 |  0.4976611317561498d0 |   0.08255362401116195d0 |
+| 1 |         0.5 | 0.41510750774498784d0 | 0.44560523266293384d0 |     0.030497724917946d0 |
+| 1 |        0.25 | 0.41510750774498784d0 |  0.4234889628937013d0 |  0.008381455148713468d0 |
+| 1 |       0.125 | 0.41510750774498784d0 | 0.41725265825950153d0 |  0.002145150514513694d0 |
+| 1 |      0.0625 | 0.41510750774498784d0 | 0.41564710776310854d0 |    5.396000181207006d-4 |
+| 1 |     0.03125 | 0.41510750774498784d0 |  0.4152414157140871d0 |   1.3390796909928948d-4 |
+| 1 |    0.015625 | 0.41510750774498784d0 | 0.41514241394084905d0 |    3.490619586121735d-5 |
+| 1 |   0.0078125 | 0.41510750774498784d0 | 0.41510582632900395d0 |   1.6814159838896003d-6 |
+| 1 |  0.00390625 | 0.41510750774498784d0 |   0.415092913054238d0 |   1.4594690749825112d-5 |
+| 1 | 0.001953125 | 0.41510750774498784d0 |  0.4150670865046777d0 |   4.0421240310117845d-5 |
+
+* Question Twelve
+
+First we'll place a bound on $h$; looking at a graph of $f$ it's pretty obvious from the asymptotes that we don't want to go much further than $|h| = 2 - \frac{pi}{2}$.
+
+Following similar reasoning as Question Four, we can determine an optimal $h$ by computing $e_{\text{abs}}$ for the central difference, but now including a roundoff error for each time we run $f$
+such that $|f_{\text{machine}}(x) - f(x)| \le \epsilon_{\text{dblprec}}$ (we'll use double precision numbers, from HW 2 we know $\epsilon_{\text{dblprec}} \approx 2.22045 (10^{-16})$).
+
+We'll just assume $|f_{\text{machine}}(x) - f(x)| = \epsilon_{\text{dblprec}}$ so our new difference quotient becomes:
+
+\begin{align*}
+e_{\text{abs}} &= |f'(a) - (\frac{f(a+h) - f(a-h) + 2\epsilon_{\text{dblprec}}}{2h})| \\
+&= |\frac{1}{12}f'''(\xi)h^2 + \frac{\epsilon_{\text{dblprec}}}{h}|
+\end{align*}
+
+Because we bounded our $|h| = 2 - \frac{pi}{2}$ we'll find the maximum value of $f'''$ between $a - (2 - \frac{\pi}{2})$ and $a - (2 - \frac{\pi}{3})$. Using [[https://www.desmos.com/calculator/gen1zpohh2][desmos]] I found this to be -2.
+
+Thus, $e_{\text{abs}} \leq \frac{1}{6}h^2 + \frac{\epsilon_{\text{dblprec}}}{h}$. Finding the derivative:
+
+\begin{equation*}
+e' = \frac{1}{3}h - \frac{\epsilon_{\text{dblprec}}}{h^2}
+\end{equation*}
+
+And solving at $e' = 0$:
+
+\begin{equation*}
+\frac{1}{3}h = \frac{\epsilon_{\text{dblprec}}}{h^2} \Rightarrow h^3 = 3\epsilon_{\text{dblprec}} \Rightarrow h = (3\epsilon_{\text{dblprec}})^{1/3}
+\end{equation*}
+
+Which is $\approx (3(2.22045 (10^{-16}))^{\frac{1}{3}} \approx 8.7335 10^{-6}$.
diff --git a/Homework/math4610/homeworks/hw-3.pdf b/Homework/math4610/homeworks/hw-3.pdf
new file mode 100644
index 0000000..1f9bac6
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-3.pdf
diff --git a/Homework/math4610/homeworks/hw-3.tex b/Homework/math4610/homeworks/hw-3.tex
new file mode 100644
index 0000000..b3d029d
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-3.tex
@@ -0,0 +1,250 @@
+% Created 2023-10-07 Sat 14:49
+% Intended LaTeX compiler: pdflatex
+\documentclass[11pt]{article}
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+\usepackage{graphicx}
+\usepackage{longtable}
+\usepackage{wrapfig}
+\usepackage{rotating}
+\usepackage[normalem]{ulem}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{capt-of}
+\usepackage{hyperref}
+\notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+\author{Elizabeth Hunt}
+\date{\today}
+\title{HW 03}
+\hypersetup{
+ pdfauthor={Elizabeth Hunt},
+ pdftitle={HW 03},
+ pdfkeywords={},
+ pdfsubject={},
+ pdfcreator={Emacs 28.2 (Org mode 9.7-pre)}, 
+ pdflang={English}}
+\begin{document}
+
+\maketitle
+\setlength\parindent{0pt}
+
+\section{Question One}
+\label{sec:org6f2bd27}
+\subsection{Three Terms}
+\label{sec:orgeb827ff}
+\begin{align*}
+Si_3(x) &= \int_0^x \frac{s - \frac{s^3}{3!} + \frac{s^5}{5!}}{s} dx \\
+&= x - \frac{x^3}{(3!)(3)} + \frac{x^5}{(5!)(5)}
+\end{align*}
+\subsection{Five Terms}
+\label{sec:orge6a15e4}
+\begin{align*}
+Si_3(x) &= \int_0^x \frac{s - \frac{s^3}{3!} + \frac{s^5}{5!} - \frac{s^7}{7!} + \frac{s^9}{9!}}{s} dx \\
+&= x - \frac{x^3}{(3!)(3)} + \frac{x^5}{(5!)(5)} - \frac{x^7}{(7!)(7)} + \frac{s^9}{(9!)(9)}
+\end{align*}
+\subsection{Ten Terms}
+\label{sec:orge87e346}
+\begin{align*}
+Si_{10}(x) &= \int_0^x \frac{s - \frac{s^3}{3!} + \frac{s^5}{5!} - \frac{s^7}{7!} + \frac{s^9}{9!} - \frac{s^{11}}{11!} + \frac{s^{13}}{13!} - \frac{s^{15}}{15!} + \frac{s^{17}}{17!} - \frac{s^{19}}{19!}}{s} ds \\
+&= x - \frac{x^3}{(3!)(3)} + \frac{x^5}{(5!)(5)} - \frac{x^7}{(7!)(7)} + \frac{s^9}{(9!)(9)} - \frac{s^{11}}{(11!)(11)} + \frac{s^{13}}{(13!)(13)} - \frac{s^{15}}{(15!)(15)} \\
+&+ \frac{s^{17}}{(17!)(17)} - \frac{s^{19}}{(19!)(19)}
+\end{align*}
+\section{Question Three}
+\label{sec:org6e2f7fc}
+For the second term in the difference quotient, we can expand the taylor series centered at x=a:
+
+\begin{equation*}
+f(x) = f(a) + f'(a)(x-a) + \frac{f''(a)}{2}(x-a)^2 + \cdots \\
+\end{equation*}
+
+Which we substitute into the difference quotient:
+
+\begin{equation*}
+\frac{f(a) - f(a - h)}{h} = \frac{f(a) - (f(a) + f'(a)(x-a) + \frac{f''(a)}{2}(x-a)^2 + \cdots)}{h}
+\end{equation*}
+
+And subs. \(x=a-h\):
+
+\begin{align*}
+\frac{f(a) - (f(a) + f'(a)(x-a) + \frac{f''(a)}{2}(x-a)^2 + \cdots)}{h} &= -f'(a)(-1) + -\frac{1}{2}f''(a)h \\
+&= f'(a) - \frac{1}{2}f''(a)h + \cdots \\
+\end{align*}
+
+Which we now plug into the initial \(e_{\text{abs}}\):
+
+\begin{align*}
+e_{\text{abs}} &= |f'(a) - \frac{f(a) - f(a - h)}{h}| \\
+&= |f'(a) - (f'(a) +  -\frac{f''(a)}{2}h + \cdots)| \\
+&= |- \frac{1}{2}f''(a)h + \cdots | \\
+\end{align*}
+
+With the Taylor Remainder theorem we can absorb the series following the second term:
+
+\begin{equation*}
+e_{\text{abs}} = |- \frac{1}{2}f''(a)h + \cdots | = |\frac{1}{2}f''(\xi)h| \leq Ch
+\end{equation*}
+
+Thus our error is bounded linearly with \(h\).
+
+\section{Question Four}
+\label{sec:orga7d02a2}
+For the first term in the difference quotient we know, from the given notes,
+
+\begin{equation*}
+f(a+h) = f(a) + f'(a)h + \frac{1}{2}f''(a)h^2 + \frac{1}{6}f'''(a)(h^3)
+\end{equation*}
+
+And from some of the work in Question Three,
+
+\begin{equation*}
+f(a - h) = f(a) + f'(a)(-h) + \frac{1}{2}f''(a)(-h)^2 + \frac{1}{6}f'''(a)(-h^3)
+\end{equation*}
+
+We can substitute immediately into \(e_{\text{abs}} = |f'(a) - (\frac{f(a+h) - f(a-h)}{2h})|\):
+
+\begin{align*}
+e_{\text{abs}} &= |f'(a) - \frac{1}{2h}((f(a) + f'(a)h + \frac{1}{2}f''(a)h^2 + \cdots) - (f(a) - f'(a)h + \frac{1}{2}f''(a)h^2 + \cdots))| \\
+&= |f'(a) - \frac{1}{2h}(2f'(a)h + \frac{1}{6}f'''(a)h^3 + \cdots)| \\
+&= |f'(a) - f'(a) - \frac{1}{12}f'''(a)h^2 + \cdots| \\
+&= |-\frac{1}{12}f'''(a)h^2 + \cdots|
+\end{align*}
+
+Finally, with the Taylor Remainder theorem we can absorb the series following the third term:
+
+\begin{equation*}
+e_{\text{abs}} = |-\frac{1}{12}f'''(\xi)h^2| = |\frac{1}{12}f'''(\xi)h^2| \leq Ch^2
+\end{equation*}
+
+Meaning that as \(h\) scales linearly, our error is bounded by \(h^2\) as opposed to linearly as in Question Three.
+
+\section{Question Six}
+\label{sec:org7b05811}
+\subsection{A}
+\label{sec:org8341a77}
+\begin{verbatim}
+(load "../lizfcm.asd")
+(ql:quickload :lizfcm)
+
+(defun f (x)
+  (/ (- x 1) (+ x 1)))
+
+(defun fprime (x)
+  (/ 2 (expt (+ x 1) 2)))
+
+(let ((domain-values (loop for a from 0 to 2
+                           append 
+                           (loop for i from 0 to 9
+                                 for h = (/ 1.0 (expt 2 i))
+                                 collect (list a h)))))
+  (lizfcm.utils:table (:headers '("a" "h" "f'" "\\approx f'" "e_{\\text{abs}}")
+                       :domain-order (a h)
+                       :domain-values domain-values)
+    (fprime a)
+    (lizfcm.approx:fwd-derivative-at 'f a h)
+    (abs (- (fprime a)
+            (lizfcm.approx:fwd-derivative-at 'f a h)))))
+\end{verbatim}
+
+
+\section{Question Nine}
+\label{sec:orgeb1839f}
+\subsection{C}
+\label{sec:org5691277}
+
+\begin{verbatim}
+(load "../lizfcm.asd")
+(ql:quickload :lizfcm)
+
+(defun factorial (n)
+  (if (= n 0)
+      1
+      (* n (factorial (- n 1)))))
+
+(defun taylor-term (n x)
+  (/ (* (expt (- 1) n)
+        (expt x (+ (* 2 n) 1)))
+     (* (factorial n)
+        (+ (* 2 n) 1))))
+
+(defun f (x &optional (max-iterations 30))
+  (let ((sum 0.0))
+    (dotimes (n max-iterations)
+      (setq sum (+ sum (taylor-term n x))))
+    (* sum (/ 2 (sqrt pi)))))
+
+(defun fprime (x)
+  (* (/ 2 (sqrt pi)) (exp (- 0 (* x x)))))
+
+(let ((domain-values (loop for a from 0 to 1
+                           append 
+                           (loop for i from 0 to 9
+                                 for h = (/ 1.0 (expt 2 i))
+                                 collect (list a h)))))
+  (lizfcm.utils:table (:headers '("a" "h" "f'" "\\approx f'" "e_{\\text{abs}}")
+                       :domain-order (a h)
+                       :domain-values domain-values)
+    (fprime a)
+    (lizfcm.approx:central-derivative-at 'f a h)
+    (abs (- (fprime a)
+            (lizfcm.approx:central-derivative-at 'f a h)))))
+\end{verbatim}
+
+
+\begin{center}
+\begin{tabular}{rrrrr}
+a & h & f' & \(\approx\) f' & e\textsubscript{\text{abs}}\\[0pt]
+0 & 1.0 & 1.1283791670955126d0 & 0.8427006725464232d0 & 0.28567849454908933d0\\[0pt]
+0 & 0.5 & 1.1283791670955126d0 & 1.0409997446922075d0 & 0.0873794224033051d0\\[0pt]
+0 & 0.25 & 1.1283791670955126d0 & 1.1053055663206806d0 & 0.023073600774832004d0\\[0pt]
+0 & 0.125 & 1.1283791670955126d0 & 1.122529655394656d0 & 0.005849511700856569d0\\[0pt]
+0 & 0.0625 & 1.1283791670955126d0 & 1.1269116944798618d0 & 0.0014674726156507223d0\\[0pt]
+0 & 0.03125 & 1.1283791670955126d0 & 1.1280120131008824d0 & 3.6715399463016496d-4\\[0pt]
+0 & 0.015625 & 1.1283791670955126d0 & 1.1282873617826952d0 & 9.180531281738347d-5\\[0pt]
+0 & 0.0078125 & 1.1283791670955126d0 & 1.128356232581468d0 & 2.293451404455915d-5\\[0pt]
+0 & 0.00390625 & 1.1283791670955126d0 & 1.1283734502811613d0 & 5.71681435124205d-6\\[0pt]
+0 & 0.001953125 & 1.1283791670955126d0 & 1.1283777547060847d0 & 1.4123894278572635d-6\\[0pt]
+1 & 1.0 & 0.41510750774498784d0 & 0.4976611317561498d0 & 0.08255362401116195d0\\[0pt]
+1 & 0.5 & 0.41510750774498784d0 & 0.44560523266293384d0 & 0.030497724917946d0\\[0pt]
+1 & 0.25 & 0.41510750774498784d0 & 0.4234889628937013d0 & 0.008381455148713468d0\\[0pt]
+1 & 0.125 & 0.41510750774498784d0 & 0.41725265825950153d0 & 0.002145150514513694d0\\[0pt]
+1 & 0.0625 & 0.41510750774498784d0 & 0.41564710776310854d0 & 5.396000181207006d-4\\[0pt]
+1 & 0.03125 & 0.41510750774498784d0 & 0.4152414157140871d0 & 1.3390796909928948d-4\\[0pt]
+1 & 0.015625 & 0.41510750774498784d0 & 0.41514241394084905d0 & 3.490619586121735d-5\\[0pt]
+1 & 0.0078125 & 0.41510750774498784d0 & 0.41510582632900395d0 & 1.6814159838896003d-6\\[0pt]
+1 & 0.00390625 & 0.41510750774498784d0 & 0.415092913054238d0 & 1.4594690749825112d-5\\[0pt]
+1 & 0.001953125 & 0.41510750774498784d0 & 0.4150670865046777d0 & 4.0421240310117845d-5\\[0pt]
+\end{tabular}
+\end{center}
+
+
+\section{Question Twelve}
+\label{sec:orgc55bfd1}
+
+First we'll place a bound on \(h\); looking at a graph of \(f\) it's pretty obvious from the asymptotes that we don't want to go much further than \(|h| = 2 - \frac{pi}{2}\).
+
+Following similar reasoning as Question Four, we can determine an optimal \(h\) by computing \(e_{\text{abs}}\) for the central difference, but now including a roundoff error for each time we run \(f\)
+such that \(|f_{\text{machine}}(x) - f(x)| \le \epsilon_{\text{dblprec}}\) (we'll use double precision numbers, from HW 2 we know \(\epsilon_{\text{dblprec}} \approx 2.22045 (10^{-16})\)).
+
+We'll just assume \(|f_{\text{machine}}(x) - f(x)| = \epsilon_{\text{dblprec}}\) so our new difference quotient becomes:
+
+\begin{align*}
+e_{\text{abs}} &= |f'(a) - (\frac{f(a+h) - f(a-h) + 2\epsilon_{\text{dblprec}}}{2h})| \\
+&= |\frac{1}{12}f'''(\xi)h^2 + \frac{\epsilon_{\text{dblprec}}}{h}|
+\end{align*}
+
+Because we bounded our \(|h| = 2 - \frac{pi}{2}\) we'll find the maximum value of \(f'''\) between \(a - (2 - \frac{\pi}{2})\) and \(a - (2 - \frac{\pi}{3})\). Using \href{https://www.desmos.com/calculator/gen1zpohh2}{desmos} I found this to be -2.
+
+Thus, \(e_{\text{abs}} \leq \frac{1}{6}h^2 + \frac{\epsilon_{\text{dblprec}}}{h}\). Finding the derivative:
+
+\begin{equation*}
+e' = \frac{1}{3}h - \frac{\epsilon_{\text{dblprec}}}{h^2}
+\end{equation*}
+
+And solving at \(e' = 0\):
+
+\begin{equation*}
+\frac{1}{3}h = \frac{\epsilon_{\text{dblprec}}}{h^2} \Rightarrow h^3 = 3\epsilon_{\text{dblprec}} \Rightarrow h = (3\epsilon_{\text{dblprec}})^{1/3}
+\end{equation*}
+
+Which is \(\approx (3(2.22045 (10^{-16}))^{\frac{1}{3}} \approx 8.7335 10^{-6}\).
+\end{document}
\ No newline at end of file
diff --git a/Homework/math4610/homeworks/hw-4.org b/Homework/math4610/homeworks/hw-4.org
new file mode 100644
index 0000000..8884f63
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-4.org
@@ -0,0 +1,34 @@
+#+TITLE: Homework 4
+#+AUTHOR: Elizabeth Hunt
+#+LATEX_HEADER: \notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+#+LATEX: \setlength\parindent{0pt}
+#+OPTIONS: toc:nil
+
+* Question 1
+See the attached LIZFCM Software Manual.
+
+* Question 2, 3, 4
+#+attr_latex: :width 350px
+[[./img/make_run.png]]
+
+* Question 5
+#+attr_latex: :width 350px
+[[./img/test_routine_1.png]]
+
+#+attr_latex: :width 350px
+[[./img/test_routine_2.png]]
+
+* Question 6
+See the LIZFCM Software Manual.
+
+* Question 7
+See ~src/matrix.c -> lu_decomp, fsubst, bsubst, solve_matrix~
+
+* Question 8
+See ~test/main.c -> lines 109 - 113~ in correspondence to the run in Question 5
+
+* Question 9
+See ~test/main.c -> lines 118 - 121~ in correspondence to the run in Question 5
+
+* Question 10
+See the TOC on the first page of the LIZFCM Software Manual.
diff --git a/Homework/math4610/homeworks/hw-4.pdf b/Homework/math4610/homeworks/hw-4.pdf
new file mode 100644
index 0000000..35176a1
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-4.pdf
diff --git a/Homework/math4610/homeworks/hw-4.tex b/Homework/math4610/homeworks/hw-4.tex
new file mode 100644
index 0000000..7a83378
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-4.tex
@@ -0,0 +1,70 @@
+% Created 2023-10-13 Fri 21:11
+% Intended LaTeX compiler: pdflatex
+\documentclass[11pt]{article}
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+\usepackage{graphicx}
+\usepackage{longtable}
+\usepackage{wrapfig}
+\usepackage{rotating}
+\usepackage[normalem]{ulem}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{capt-of}
+\usepackage{hyperref}
+\notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+\author{Elizabeth Hunt}
+\date{\today}
+\title{Homework 4}
+\hypersetup{
+ pdfauthor={Elizabeth Hunt},
+ pdftitle={Homework 4},
+ pdfkeywords={},
+ pdfsubject={},
+ pdfcreator={Emacs 28.2 (Org mode 9.7-pre)}, 
+ pdflang={English}}
+\begin{document}
+
+\maketitle
+\setlength\parindent{0pt}
+
+\section{Question 1}
+\label{sec:orgf7348d4}
+See the attached LIZFCM Software Manual.
+
+\section{Question 2, 3, 4}
+\label{sec:orgaf52510}
+\begin{center}
+\includegraphics[width=350px]{./img/make_run.png}
+\end{center}
+
+\section{Question 5}
+\label{sec:orgd0fe6e8}
+\begin{center}
+\includegraphics[width=350px]{./img/test_routine_1.png}
+\end{center}
+
+\begin{center}
+\includegraphics[width=350px]{./img/test_routine_2.png}
+\end{center}
+
+\section{Question 6}
+\label{sec:org9e2023c}
+See the LIZFCM Software Manual.
+
+\section{Question 7}
+\label{sec:org6c11571}
+See \texttt{src/matrix.c -> lu\_decomp, fsubst, bsubst, solve\_matrix}
+
+\section{Question 8}
+\label{sec:org9ba7792}
+See \texttt{test/main.c -> lines 109 - 113} in correspondence to the run in Question 5
+
+\section{Question 9}
+\label{sec:org3cff888}
+See \texttt{test/main.c -> lines 118 - 121} in correspondence to the run in Question 5
+
+\section{Question 10}
+\label{sec:org522eabc}
+See the TOC on the first page of the LIZFCM Software Manual.
+\end{document}
\ No newline at end of file
diff --git a/Homework/math4610/homeworks/hw-5.org b/Homework/math4610/homeworks/hw-5.org
new file mode 100644
index 0000000..a2339f9
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-5.org
@@ -0,0 +1,59 @@
+#+TITLE: Homework 5
+#+AUTHOR: Elizabeth Hunt
+#+LATEX_HEADER: \notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+#+LATEX: \setlength\parindent{0pt}
+#+OPTIONS: toc:nil
+
+* Question One
+See LIZFCM \rightarrow Matrix Routines \rightarrow ~lu decomp~ & ~bsubst~.
+
+The test ~UTEST(matrix, lu_decomp)~ is a unit test for the ~lu_decomp~ routine,
+and ~UTEST(matrix, bsubst)~ verifies back substitution on an upper triangular
+3 \times 3 matrix with a known solution that can be verified manually.
+
+Both can be found in ~tests/matrix.t.c~.
+
+* Question Two
+Unless the following are met, the resulting solution will be garbage.
+
+1. The matrix $U$ must be not be singular.
+2. $U$ must be square (or it will fail the ~assert~).
+3. The system created by $Ux = b$ must be consistent.
+4. $U$ is (quite obviously) in upper-triangular form.
+
+Thus, the actual calculation performing the $LU$ decomposition
+(in ~lu_decomp~) does a sanity
+check for 1-3 will fail an assert, should a point along the diagonal (pivot) be
+zero, or the matrix be non-factorable.
+
+* Question Three
+See LIZFCM \rightarrow Matrix Routines \rightarrow ~fsubst~.
+
+~UTEST(matrix, fsubst)~ verifies forward substitution on a lower triangular 3 \times 3
+matrix with a known solution that can be verified manually.
+
+* Question Four
+
+See LIZFCM \rightarrow Matrix Routines \rightarrow ~gaussian_elimination~ and ~solve_gaussian_elimination~.
+
+* Question Five
+See LIZFCM \rightarrow Matrix Routines \rightarrow ~m_dot_v~, and the ~UTEST(matrix, m_dot_v)~ in
+~tests/matrix.t.c~.
+
+* Question Six
+See ~UTEST(matrix, solve_gaussian_elimination)~ in ~tests/matrix.t.c~, which generates a diagonally dominant 10 \times 10 matrix
+and shows that the solution is consistent with the initial matrix, according to the steps given. Then,
+we do a dot product between each row of the diagonally dominant matrix and the solution vector to ensure
+it is near equivalent to the input vector.
+
+* Question Seven
+See ~UTEST(matrix, solve_matrix_lu_bsubst)~ which does the same test in Question Six with the solution according to
+~solve_matrix_lu_bsubst~ as shown in the Software Manual.
+
+* Question Eight
+No, since the time complexity for Gaussian Elimination is always less than that of the LU factorization solution by $O(n^2)$ operations
+(in LU factorization we perform both backwards and forwards substitutions proceeding the LU decomp, in Gaussian Elimination we only need
+back substitution).
+
+* Question Nine, Ten
+See LIZFCM Software manual and shared library in ~dist~ after compiling.
diff --git a/Homework/math4610/homeworks/hw-5.pdf b/Homework/math4610/homeworks/hw-5.pdf
new file mode 100644
index 0000000..a7773bc
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-5.pdf
diff --git a/Homework/math4610/homeworks/hw-5.tex b/Homework/math4610/homeworks/hw-5.tex
new file mode 100644
index 0000000..98cca2e
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-5.tex
@@ -0,0 +1,95 @@
+% Created 2023-11-01 Wed 20:49
+% Intended LaTeX compiler: pdflatex
+\documentclass[11pt]{article}
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+\usepackage{graphicx}
+\usepackage{longtable}
+\usepackage{wrapfig}
+\usepackage{rotating}
+\usepackage[normalem]{ulem}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{capt-of}
+\usepackage{hyperref}
+\notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+\author{Elizabeth Hunt}
+\date{\today}
+\title{Homework 5}
+\hypersetup{
+ pdfauthor={Elizabeth Hunt},
+ pdftitle={Homework 5},
+ pdfkeywords={},
+ pdfsubject={},
+ pdfcreator={Emacs 28.2 (Org mode 9.7-pre)}, 
+ pdflang={English}}
+\begin{document}
+
+\maketitle
+\setlength\parindent{0pt}
+
+\section{Question One}
+\label{sec:org4e80298}
+See LIZFCM \(\rightarrow\) Matrix Routines \(\rightarrow\) \texttt{lu decomp} \& \texttt{bsubst}.
+
+The test \texttt{UTEST(matrix, lu\_decomp)} is a unit test for the \texttt{lu\_decomp} routine,
+and \texttt{UTEST(matrix, bsubst)} verifies back substitution on an upper triangular
+3 \texttimes{} 3 matrix with a known solution that can be verified manually.
+
+Both can be found in \texttt{tests/matrix.t.c}.
+
+\section{Question Two}
+\label{sec:orga73d05c}
+Unless the following are met, the resulting solution will be garbage.
+
+\begin{enumerate}
+\item The matrix \(U\) must be not be singular.
+\item \(U\) must be square (or it will fail the \texttt{assert}).
+\item The system created by \(Ux = b\) must be consistent.
+\item \(U\) is (quite obviously) in upper-triangular form.
+\end{enumerate}
+
+Thus, the actual calculation performing the \(LU\) decomposition
+(in \texttt{lu\_decomp}) does a sanity
+check for 1-3 will fail an assert, should a point along the diagonal (pivot) be
+zero, or the matrix be non-factorable.
+
+\section{Question Three}
+\label{sec:org35163c5}
+See LIZFCM \(\rightarrow\) Matrix Routines \(\rightarrow\) \texttt{fsubst}.
+
+\texttt{UTEST(matrix, fsubst)} verifies forward substitution on a lower triangular 3 \texttimes{} 3
+matrix with a known solution that can be verified manually.
+
+\section{Question Four}
+\label{sec:org79d9061}
+
+See LIZFCM \(\rightarrow\) Matrix Routines \(\rightarrow\) \texttt{gaussian\_elimination} and \texttt{solve\_gaussian\_elimination}.
+
+\section{Question Five}
+\label{sec:orgc6ac464}
+See LIZFCM \(\rightarrow\) Matrix Routines \(\rightarrow\) \texttt{m\_dot\_v}, and the \texttt{UTEST(matrix, m\_dot\_v)} in
+\texttt{tests/matrix.t.c}.
+
+\section{Question Six}
+\label{sec:org66fedab}
+See \texttt{UTEST(matrix, solve\_gaussian\_elimination)} in \texttt{tests/matrix.t.c}, which generates a diagonally dominant 10 \texttimes{} 10 matrix
+and shows that the solution is consistent with the initial matrix, according to the steps given. Then,
+we do a dot product between each row of the diagonally dominant matrix and the solution vector to ensure
+it is near equivalent to the input vector.
+
+\section{Question Seven}
+\label{sec:org6897ff2}
+See \texttt{UTEST(matrix, solve\_matrix\_lu\_bsubst)} which does the same test in Question Six with the solution according to
+\texttt{solve\_matrix\_lu\_bsubst} as shown in the Software Manual.
+
+\section{Question Eight}
+\label{sec:org5d529dd}
+No, since the time complexity for Gaussian Elimination is always less than that of the LU factorization solution by \(O(n^2)\) operations
+(in LU factorization we perform both backwards and forwards substitutions proceeding the LU decomp, in Gaussian Elimination we only need
+back substitution).
+
+\section{Question Nine, Ten}
+\label{sec:org0fb8e09}
+See LIZFCM Software manual and shared library in \texttt{dist} after compiling.
+\end{document}
\ No newline at end of file
diff --git a/Homework/math4610/homeworks/hw-6.org b/Homework/math4610/homeworks/hw-6.org
new file mode 100644
index 0000000..eebc0c2
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-6.org
@@ -0,0 +1,199 @@
+#+TITLE: Homework 6
+#+AUTHOR: Elizabeth Hunt
+#+LATEX_HEADER: \notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+#+LATEX: \setlength\parindent{0pt}
+#+OPTIONS: toc:nil
+
+* Question One
+
+For $g(x) = x + f(x)$ then we know $g'(x) = 1 + 2x - 5$ and thus $|g'(x)| \lt 1$ is only true
+on the interval $(1.5, 2.5)$, and for $g(x) = x - f(x)$ then we know $g'(x) = 1 - (2x - 5)$
+and thus $|g'(x)| < 1$ is only true on the interval $(2.5, 3.5)$.
+
+Because we know the roots of $f$ are $2, 3$ ($f(x) = (x-2)(x-3)$) then we can only be
+certain that $g(x) = x + f(x)$ will converge to the root $2$ if we pick an initial
+guess between $(1.5, 2.5)$, and likewise for $g(x) = x - f(x)$, $3$:
+
+#+BEGIN_SRC c
+  // tests/roots.t.c
+  UTEST(root, fixed_point_iteration_method) {
+    // x^2 - 5x + 6 = (x - 3)(x - 2)
+    double expect_x1 = 3.0;
+    double expect_x2 = 2.0;
+
+    double tolerance = 0.001;
+    uint64_t max_iterations = 10;
+
+    double x_0 = 1.55; // 1.5 < 1.55 < 2.5
+    // g1(x) = x + f(x)
+    double root1 =
+        fixed_point_iteration_method(&f2, &g1, x_0, tolerance, max_iterations);
+    EXPECT_NEAR(root1, expect_x2, tolerance);
+
+    // g2(x) = x - f(x)
+    x_0 = 3.4; // 2.5 < 3.4 < 3.5
+    double root2 =
+        fixed_point_iteration_method(&f2, &g2, x_0, tolerance, max_iterations);
+    EXPECT_NEAR(root2, expect_x1, tolerance);
+  }
+#+END_SRC
+
+And by this method passing in ~tests/roots.t.c~ we know they converged within ~tolerance~ before
+10 iterations.
+
+* Question Two
+
+Yes, we showed that for $\epsilon = 1$ in Question One, we can converge upon a root in the range $(2.5, 3.5)$, and
+when $\epsilon = -1$ we can converge upon a root in the range $(1.5, 2.5)$.
+
+See the above unit tests in Question One for each $\epsilon$.
+
+* Question Three
+
+See ~test/roots.t.c -> UTEST(root, bisection_with_error_assumption)~
+and the software manual entry ~bisect_find_root_with_error_assumption~.
+
+* Question Four
+
+See ~test/roots.t.c -> UTEST(root, fixed_point_newton_method)~
+and the software manual entry ~fixed_point_newton_method~.
+
+* Question Five
+
+See ~test/roots.t.c -> UTEST(root, fixed_point_secant_method)~
+and the software manual entry ~fixed_point_secant_method~.
+
+* Question Six
+
+See ~test/roots.t.c -> UTEST(root, fixed_point_bisection_secant_method)~
+and the software manual entry ~fixed_point_bisection_secant_method~.
+
+* Question Seven
+
+The existance of ~test/roots.t.c~'s compilation into ~dist/lizfcm.test~ via ~make~
+shows that the compiled ~lizfcm.a~ contains the root methods mentioned; a user
+could link the library and use them, as we do in Question Eight.
+
+* Question Eight
+
+The given ODE $\frac{dP}{dt} = \alpha P - \beta P$ has a trivial solution by separation:
+
+\begin{equation*}
+P(t) = C e^{t(\alpha - \beta)}
+\end{equation*}
+
+And
+
+\begin{equation*}
+P_0 = P(0) = C e^0 = C
+\end{equation*}
+
+So $P(t) = P_0 e^{t(\alpha - \beta)}$.
+
+We're trying to find $t$ such that $P(t) = P_\infty$, thus we're finding roots of $P(t) - P_\infty$.
+
+The following code (in ~homeworks/hw_6_p_8.c~) produces this output:
+
+\begin{verbatim}
+$ gcc -I../inc/ -Wall hw_6_p_8.c ../lib/lizfcm.a -lm -o hw_6_p_8 && ./hw_6_p_8
+a ~ 27.269515; P(27.269515) - P_infty = -0.000000
+b ~ 40.957816; P(40.957816) - P_infty = -0.000000
+c ~ 40.588827; P(40.588827) - P_infty = -0.000000
+d ~ 483.611967; P(483.611967) - P_infty = -0.000000
+e ~ 4.894274; P(4.894274) - P_infty = -0.000000
+
+\end{verbatim}
+
+#+BEGIN_SRC c
+// compile & test w/
+//  \--> gcc -I../inc/ -Wall hw_6_p_8.c ../lib/lizfcm.a -lm -o hw_6_p_8
+//  \--> ./hw_6_p_8
+
+#include "lizfcm.h"
+#include <math.h>
+#include <stdio.h>
+
+double a(double t) {
+  double alpha = 0.1;
+  double beta = 0.001;
+  double p_0 = 2;
+  double p_infty = 29.75;
+
+  return p_0 * exp(t * (alpha - beta)) - p_infty;
+}
+
+double b(double t) {
+  double alpha = 0.1;
+  double beta = 0.001;
+  double p_0 = 2;
+  double p_infty = 115.35;
+
+  return p_0 * exp(t * (alpha - beta)) - p_infty;
+}
+
+double c(double t) {
+  double alpha = 0.1;
+  double beta = 0.0001;
+  double p_0 = 2;
+  double p_infty = 115.35;
+
+  return p_0 * exp(t * (alpha - beta)) - p_infty;
+}
+
+double d(double t) {
+  double alpha = 0.01;
+  double beta = 0.001;
+  double p_0 = 2;
+  double p_infty = 155.346;
+
+  return p_0 * exp(t * (alpha - beta)) - p_infty;
+}
+
+double e(double t) {
+  double alpha = 0.1;
+  double beta = 0.01;
+  double p_0 = 100;
+  double p_infty = 155.346;
+
+  return p_0 * exp(t * (alpha - beta)) - p_infty;
+}
+
+int main() {
+  uint64_t max_iterations = 1000;
+  double tolerance = 0.0000001;
+
+  Array_double *ivt_range = find_ivt_range(&a, -5.0, 3.0, 1000);
+  double approx_a = fixed_point_secant_bisection_method(
+      &a, ivt_range->data[0], ivt_range->data[1], tolerance, max_iterations);
+
+  free_vector(ivt_range);
+  ivt_range = find_ivt_range(&b, -5.0, 3.0, 1000);
+  double approx_b = fixed_point_secant_bisection_method(
+      &b, ivt_range->data[0], ivt_range->data[1], tolerance, max_iterations);
+
+  free_vector(ivt_range);
+  ivt_range = find_ivt_range(&c, -5.0, 3.0, 1000);
+  double approx_c = fixed_point_secant_bisection_method(
+      &c, ivt_range->data[0], ivt_range->data[1], tolerance, max_iterations);
+
+  free_vector(ivt_range);
+  ivt_range = find_ivt_range(&d, -5.0, 3.0, 1000);
+  double approx_d = fixed_point_secant_bisection_method(
+      &d, ivt_range->data[0], ivt_range->data[1], tolerance, max_iterations);
+
+  free_vector(ivt_range);
+  ivt_range = find_ivt_range(&e, -5.0, 3.0, 1000);
+  double approx_e = fixed_point_secant_bisection_method(
+      &e, ivt_range->data[0], ivt_range->data[1], tolerance, max_iterations);
+
+  printf("a ~ %f; P(%f) = %f\n", approx_a, approx_a, a(approx_a));
+  printf("b ~ %f; P(%f) = %f\n", approx_b, approx_b, b(approx_b));
+  printf("c ~ %f; P(%f) = %f\n", approx_c, approx_c, c(approx_c));
+  printf("d ~ %f; P(%f) = %f\n", approx_d, approx_d, d(approx_d));
+  printf("e ~ %f; P(%f) = %f\n", approx_e, approx_e, e(approx_e));
+
+  return 0;
+}
+#+END_SRC
+
+
diff --git a/Homework/math4610/homeworks/hw-6.pdf b/Homework/math4610/homeworks/hw-6.pdf
new file mode 100644
index 0000000..c056102
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-6.pdf
diff --git a/Homework/math4610/homeworks/hw-6.tex b/Homework/math4610/homeworks/hw-6.tex
new file mode 100644
index 0000000..1a0ddc4
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-6.tex
@@ -0,0 +1,223 @@
+% Created 2023-11-11 Sat 13:13
+% Intended LaTeX compiler: pdflatex
+\documentclass[11pt]{article}
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+\usepackage{graphicx}
+\usepackage{longtable}
+\usepackage{wrapfig}
+\usepackage{rotating}
+\usepackage[normalem]{ulem}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{capt-of}
+\usepackage{hyperref}
+\notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+\author{Elizabeth Hunt}
+\date{\today}
+\title{Homework 6}
+\hypersetup{
+ pdfauthor={Elizabeth Hunt},
+ pdftitle={Homework 6},
+ pdfkeywords={},
+ pdfsubject={},
+ pdfcreator={Emacs 29.1 (Org mode 9.7-pre)}, 
+ pdflang={English}}
+\begin{document}
+
+\maketitle
+\setlength\parindent{0pt}
+\section{Question One}
+\label{sec:org206b859}
+
+For \(g(x) = x + f(x)\) then we know \(g'(x) = 1 + 2x - 5\) and thus \(|g'(x)| \lt 1\) is only true
+on the interval \((1.5, 2.5)\), and for \(g(x) = x - f(x)\) then we know \(g'(x) = 1 - (2x - 5)\)
+and thus \(|g'(x)| < 1\) is only true on the interval \((2.5, 3.5)\).
+
+Because we know the roots of \(f\) are \(2, 3\) (\(f(x) = (x-2)(x-3)\)) then we can only be
+certain that \(g(x) = x + f(x)\) will converge to the root \(2\) if we pick an initial
+guess between \((1.5, 2.5)\), and likewise for \(g(x) = x - f(x)\), \(3\):
+
+\begin{verbatim}
+// tests/roots.t.c
+UTEST(root, fixed_point_iteration_method) {
+  // x^2 - 5x + 6 = (x - 3)(x - 2)
+  double expect_x1 = 3.0;
+  double expect_x2 = 2.0;
+
+  double tolerance = 0.001;
+  uint64_t max_iterations = 10;
+
+  double x_0 = 1.55; // 1.5 < 1.55 < 2.5
+  // g1(x) = x + f(x)
+  double root1 =
+      fixed_point_iteration_method(&f2, &g1, x_0, tolerance, max_iterations);
+  EXPECT_NEAR(root1, expect_x2, tolerance);
+
+  // g2(x) = x - f(x)
+  x_0 = 3.4; // 2.5 < 3.4 < 3.5
+  double root2 =
+      fixed_point_iteration_method(&f2, &g2, x_0, tolerance, max_iterations);
+  EXPECT_NEAR(root2, expect_x1, tolerance);
+}
+\end{verbatim}
+
+And by this method passing in \texttt{tests/roots.t.c} we know they converged within \texttt{tolerance} before
+10 iterations.
+\section{Question Two}
+\label{sec:orga0f5b42}
+
+Yes, we showed that for \(\epsilon = 1\) in Question One, we can converge upon a root in the range \((2.5, 3.5)\), and
+when \(\epsilon = -1\) we can converge upon a root in the range \((1.5, 2.5)\).
+
+See the above unit tests in Question One for each \(\epsilon\).
+\section{Question Three}
+\label{sec:org19aa326}
+
+See \texttt{test/roots.t.c -> UTEST(root, bisection\_with\_error\_assumption)}
+and the software manual entry \texttt{bisect\_find\_root\_with\_error\_assumption}.
+\section{Question Four}
+\label{sec:org722aa6a}
+
+See \texttt{test/roots.t.c -> UTEST(root, fixed\_point\_newton\_method)}
+and the software manual entry \texttt{fixed\_point\_newton\_method}.
+\section{Question Five}
+\label{sec:org587ee52}
+
+See \texttt{test/roots.t.c -> UTEST(root, fixed\_point\_secant\_method)}
+and the software manual entry \texttt{fixed\_point\_secant\_method}.
+\section{Question Six}
+\label{sec:org79bf754}
+
+See \texttt{test/roots.t.c -> UTEST(root, fixed\_point\_bisection\_secant\_method)}
+and the software manual entry \texttt{fixed\_point\_bisection\_secant\_method}.
+\section{Question Seven}
+\label{sec:org4cb47e5}
+
+The existance of \texttt{test/roots.t.c}'s compilation into \texttt{dist/lizfcm.test} via \texttt{make}
+shows that the compiled \texttt{lizfcm.a} contains the root methods mentioned; a user
+could link the library and use them, as we do in Question Eight.
+\section{Question Eight}
+\label{sec:org4a8160d}
+
+The given ODE \(\frac{dP}{dt} = \alpha P - \beta P\) has a trivial solution by separation:
+
+\begin{equation*}
+P(t) = C e^{t(\alpha - \beta)}
+\end{equation*}
+
+And
+
+\begin{equation*}
+P_0 = P(0) = C e^0 = C
+\end{equation*}
+
+So \(P(t) = P_0 e^{t(\alpha - \beta)}\).
+
+We're trying to find \(t\) such that \(P(t) = P_\infty\), thus we're finding roots of \(P(t) - P_\infty\).
+
+The following code (in \texttt{homeworks/hw\_6\_p\_8.c}) produces this output:
+
+\begin{verbatim}
+$ gcc -I../inc/ -Wall hw_6_p_8.c ../lib/lizfcm.a -lm -o hw_6_p_8 && ./hw_6_p_8
+
+a ~ 27.303411; P(27.303411) - P_infty = -0.000000
+b ~ 40.957816; P(40.957816) - P_infty = -0.000000
+c ~ 40.588827; P(40.588827) - P_infty = -0.000000
+d ~ 483.611967; P(483.611967) - P_infty = -0.000000
+e ~ 4.894274; P(4.894274) - P_infty = -0.000000
+
+\end{verbatim}
+
+\begin{verbatim}
+// compile & test w/
+//  \--> gcc -I../inc/ -Wall hw_6_p_8.c ../lib/lizfcm.a -lm -o hw_6_p_8
+//  \--> ./hw_6_p_8
+
+#include "lizfcm.h"
+#include <math.h>
+#include <stdio.h>
+
+double a(double t) {
+  double alpha = 0.1;
+  double beta = 0.001;
+  double p_0 = 2;
+  double p_infty = 29.85;
+
+  return p_0 * exp(t * (alpha - beta)) - p_infty;
+}
+
+double b(double t) {
+  double alpha = 0.1;
+  double beta = 0.001;
+  double p_0 = 2;
+  double p_infty = 115.35;
+
+  return p_0 * exp(t * (alpha - beta)) - p_infty;
+}
+
+double c(double t) {
+  double alpha = 0.1;
+  double beta = 0.0001;
+  double p_0 = 2;
+  double p_infty = 115.35;
+
+  return p_0 * exp(t * (alpha - beta)) - p_infty;
+}
+
+double d(double t) {
+  double alpha = 0.01;
+  double beta = 0.001;
+  double p_0 = 2;
+  double p_infty = 155.346;
+
+  return p_0 * exp(t * (alpha - beta)) - p_infty;
+}
+
+double e(double t) {
+  double alpha = 0.1;
+  double beta = 0.01;
+  double p_0 = 100;
+  double p_infty = 155.346;
+
+  return p_0 * exp(t * (alpha - beta)) - p_infty;
+}
+
+int main() {
+  uint64_t max_iterations = 1000;
+  double tolerance = 0.0000001;
+
+  Array_double *ivt_range = find_ivt_range(&a, -5.0, 3.0, 1000);
+  double approx_a = fixed_point_secant_bisection_method(
+      &a, ivt_range->data[0], ivt_range->data[1], tolerance, max_iterations);
+
+  free_vector(ivt_range);
+  ivt_range = find_ivt_range(&b, -5.0, 3.0, 1000);
+  double approx_b = fixed_point_secant_bisection_method(
+      &b, ivt_range->data[0], ivt_range->data[1], tolerance, max_iterations);
+
+  free_vector(ivt_range);
+  ivt_range = find_ivt_range(&c, -5.0, 3.0, 1000);
+  double approx_c = fixed_point_secant_bisection_method(
+      &c, ivt_range->data[0], ivt_range->data[1], tolerance, max_iterations);
+
+  free_vector(ivt_range);
+  ivt_range = find_ivt_range(&d, -5.0, 3.0, 1000);
+  double approx_d = fixed_point_secant_bisection_method(
+      &d, ivt_range->data[0], ivt_range->data[1], tolerance, max_iterations);
+
+  free_vector(ivt_range);
+  ivt_range = find_ivt_range(&e, -5.0, 3.0, 1000);
+  double approx_e = fixed_point_secant_bisection_method(
+      &e, ivt_range->data[0], ivt_range->data[1], tolerance, max_iterations);
+
+  printf("a ~ %f; P(%f) = %f\n", approx_a, approx_a, a(approx_a));
+  printf("b ~ %f; P(%f) = %f\n", approx_b, approx_b, b(approx_b));
+  printf("c ~ %f; P(%f) = %f\n", approx_c, approx_c, c(approx_c));
+  printf("d ~ %f; P(%f) = %f\n", approx_d, approx_d, d(approx_d));
+  printf("e ~ %f; P(%f) = %f\n", approx_e, approx_e, e(approx_e));
+
+  return 0;
+}
+\end{verbatim}
+\end{document}
diff --git a/Homework/math4610/homeworks/hw-7.org b/Homework/math4610/homeworks/hw-7.org
new file mode 100644
index 0000000..2c28af2
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-7.org
@@ -0,0 +1,76 @@
+#+TITLE: Homework 7
+#+AUTHOR: Elizabeth Hunt
+#+LATEX_HEADER: \notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+#+LATEX: \setlength\parindent{0pt}
+#+OPTIONS: toc:nil
+
+* Question One
+See ~UTEST(eigen, dominant_eigenvalue)~ in ~test/eigen.t.c~ and the entry
+~Eigen-Adjacent -> dominant_eigenvalue~ in the LIZFCM API documentation.
+* Question Two
+See ~UTEST(eigen, leslie_matrix_dominant_eigenvalue)~ in ~test/eigen.t.c~
+and the entry ~Eigen-Adjacent -> leslie_matrix~ in the LIZFCM API
+documentation.
+* Question Three
+See ~UTEST(eigen, least_dominant_eigenvalue)~ in ~test/eigen.t.c~ which
+finds the least dominant eigenvalue on the matrix:
+
+\begin{bmatrix}
+2 & 2 & 4 \\
+1 & 4 & 7 \\
+0 & 2 & 6 
+\end{bmatrix}
+
+which has eigenvalues: $5 + \sqrt{17}, 2, 5 - \sqrt{17}$ and should thus produce $5 - \sqrt{17}$.
+
+See also the entry ~Eigen-Adjacent -> least_dominant_eigenvalue~ in the LIZFCM API
+documentation.
+* Question Four
+See ~UTEST(eigen, shifted_eigenvalue)~ in ~test/eigen.t.c~ which
+finds the least dominant eigenvalue on the matrix:
+
+\begin{bmatrix}
+2 & 2 & 4 \\
+1 & 4 & 7 \\
+0 & 2 & 6 
+\end{bmatrix}
+
+which has eigenvalues: $5 + \sqrt{17}, 2, 5 - \sqrt{17}$ and should thus produce $2.0$.
+
+With the initial guess: $[0.5, 1.0, 0.75]$.
+
+See also the entry ~Eigen-Adjacent -> shift_inverse_power_eigenvalue~ in the LIZFCM API
+documentation.
+* Question Five
+See ~UTEST(eigen, partition_find_eigenvalues)~ in ~test/eigen.t.c~ which
+finds the eigenvalues in a partition of 10 on the matrix:
+
+\begin{bmatrix}
+2 & 2 & 4 \\
+1 & 4 & 7 \\
+0 & 2 & 6 
+\end{bmatrix}
+
+which has eigenvalues: $5 + \sqrt{17}, 2, 5 - \sqrt{17}$, and should produce all three from
+the partitions when given the guesses $[0.5, 1.0, 0.75]$ from the questions above.
+
+See also the entry ~Eigen-Adjacent -> partition_find_eigenvalues~ in the LIZFCM API
+documentation.
+
+* Question Six
+Consider we have the results of two methods developed in this homework: ~least_dominant_eigenvalue~, and ~dominant_eigenvalue~
+into ~lambda_0~, ~lambda_n~, respectively. Also assume that we have the method implemented as we've introduced,
+~shift_inverse_power_eigenvalue~.
+
+Then, we begin at the midpoint of ~lambda_0~ and ~lambda_n~, and compute the
+~new_lambda = shift_inverse_power_eigenvalue~
+with a shift at the midpoint, and some given initial guess.
+
+1. If the result is equal (or within some tolerance) to ~lambda_n~ then the closest eigenvalue to the midpoint
+   is still the dominant eigenvalue, and thus the next most dominant will be on the left. Set ~lambda_n~
+   to the midpoint and reiterate.
+2. If the result is greater or equal to ~lambda_0~ we know an eigenvalue of greater or equal magnitude
+   exists on the right. So, we set ~lambda_0~ to this eigenvalue associated with the midpoint, and
+   re-iterate.
+3. Continue re-iterating until we hit some given maximum number of iterations. Finally we will return
+   ~new_lambda~.
diff --git a/Homework/math4610/homeworks/hw-7.pdf b/Homework/math4610/homeworks/hw-7.pdf
new file mode 100644
index 0000000..4003f59
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-7.pdf
diff --git a/Homework/math4610/homeworks/hw-7.tex b/Homework/math4610/homeworks/hw-7.tex
new file mode 100644
index 0000000..be3fde4
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-7.tex
@@ -0,0 +1,107 @@
+% Created 2023-11-27 Mon 15:13
+% Intended LaTeX compiler: pdflatex
+\documentclass[11pt]{article}
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+\usepackage{graphicx}
+\usepackage{longtable}
+\usepackage{wrapfig}
+\usepackage{rotating}
+\usepackage[normalem]{ulem}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{capt-of}
+\usepackage{hyperref}
+\notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+\author{Elizabeth Hunt}
+\date{\today}
+\title{Homework 7}
+\hypersetup{
+ pdfauthor={Elizabeth Hunt},
+ pdftitle={Homework 7},
+ pdfkeywords={},
+ pdfsubject={},
+ pdfcreator={Emacs 29.1 (Org mode 9.7-pre)}, 
+ pdflang={English}}
+\begin{document}
+
+\maketitle
+\setlength\parindent{0pt}
+\section{Question One}
+\label{sec:org8ef0ee6}
+See \texttt{UTEST(eigen, dominant\_eigenvalue)} in \texttt{test/eigen.t.c} and the entry
+\texttt{Eigen-Adjacent -> dominant\_eigenvalue} in the LIZFCM API documentation.
+\section{Question Two}
+\label{sec:orgbdba5c1}
+See \texttt{UTEST(eigen, leslie\_matrix\_dominant\_eigenvalue)} in \texttt{test/eigen.t.c}
+and the entry \texttt{Eigen-Adjacent -> leslie\_matrix} in the LIZFCM API
+documentation.
+\section{Question Three}
+\label{sec:org19b04f4}
+See \texttt{UTEST(eigen, least\_dominant\_eigenvalue)} in \texttt{test/eigen.t.c} which
+finds the least dominant eigenvalue on the matrix:
+
+\begin{bmatrix}
+2 & 2 & 4 \\
+1 & 4 & 7 \\
+0 & 2 & 6 
+\end{bmatrix}
+
+which has eigenvalues: \(5 + \sqrt{17}, 2, 5 - \sqrt{17}\) and should thus produce \(5 - \sqrt{17}\).
+
+See also the entry \texttt{Eigen-Adjacent -> least\_dominant\_eigenvalue} in the LIZFCM API
+documentation.
+\section{Question Four}
+\label{sec:orgc58d42d}
+See \texttt{UTEST(eigen, shifted\_eigenvalue)} in \texttt{test/eigen.t.c} which
+finds the least dominant eigenvalue on the matrix:
+
+\begin{bmatrix}
+2 & 2 & 4 \\
+1 & 4 & 7 \\
+0 & 2 & 6 
+\end{bmatrix}
+
+which has eigenvalues: \(5 + \sqrt{17}, 2, 5 - \sqrt{17}\) and should thus produce \(2.0\).
+
+With the initial guess: \([0.5, 1.0, 0.75]\).
+
+See also the entry \texttt{Eigen-Adjacent -> shift\_inverse\_power\_eigenvalue} in the LIZFCM API
+documentation.
+\section{Question Five}
+\label{sec:orga369221}
+See \texttt{UTEST(eigen, partition\_find\_eigenvalues)} in \texttt{test/eigen.t.c} which
+finds the eigenvalues in a partition of 10 on the matrix:
+
+\begin{bmatrix}
+2 & 2 & 4 \\
+1 & 4 & 7 \\
+0 & 2 & 6 
+\end{bmatrix}
+
+which has eigenvalues: \(5 + \sqrt{17}, 2, 5 - \sqrt{17}\), and should produce all three from
+the partitions when given the guesses \([0.5, 1.0, 0.75]\) from the questions above.
+
+See also the entry \texttt{Eigen-Adjacent -> partition\_find\_eigenvalues} in the LIZFCM API
+documentation.
+\section{Question Six}
+\label{sec:orgadc3078}
+Consider we have the results of two methods developed in this homework: \texttt{least\_dominant\_eigenvalue}, and \texttt{dominant\_eigenvalue}
+into \texttt{lambda\_0}, \texttt{lambda\_n}, respectively. Also assume that we have the method implemented as we've introduced,
+\texttt{shift\_inverse\_power\_eigenvalue}.
+
+Then, we begin at the midpoint of \texttt{lambda\_0} and \texttt{lambda\_n}, and compute the
+\texttt{new\_lambda = shift\_inverse\_power\_eigenvalue}
+with a shift at the midpoint, and some given initial guess.
+
+\begin{enumerate}
+\item If the result is equal (or within some tolerance) to \texttt{lambda\_n} then the closest eigenvalue to the midpoint
+is still the dominant eigenvalue, and thus the next most dominant will be on the left. Set \texttt{lambda\_n}
+to the midpoint and reiterate.
+\item If the result is greater or equal to \texttt{lambda\_0} we know an eigenvalue of greater or equal magnitude
+exists on the right. So, we set \texttt{lambda\_0} to this eigenvalue associated with the midpoint, and
+re-iterate.
+\item Continue re-iterating until we hit some given maximum number of iterations. Finally we will return
+\texttt{new\_lambda}.
+\end{enumerate}
+\end{document}
diff --git a/Homework/math4610/homeworks/hw-8.org b/Homework/math4610/homeworks/hw-8.org
new file mode 100644
index 0000000..10c7dd8
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-8.org
@@ -0,0 +1,311 @@
+#+TITLE: Homework 8
+#+AUTHOR: Elizabeth Hunt
+#+LATEX_HEADER: \notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+#+LATEX: \setlength\parindent{0pt}
+#+OPTIONS: toc:nil
+
+* Question One
+See ~UTEST(jacobi, solve_jacobi)~ in ~test/jacobi.t.c~ and the entry
+~Jacobi / Gauss-Siedel -> solve_jacobi~ in the LIZFCM API documentation.
+* Question Two
+We cannot just perform the Jacobi algorithm on a Leslie matrix since
+it is obviously not diagonally dominant - which is a requirement. It is
+certainly not always the case, but, if a Leslie matrix $L$ is invertible, we can
+first perform gaussian elimination on $L$ augmented with $n_{k+1}$
+to obtain $n_k$ with the Jacobi method. See ~UTEST(jacobi, leslie_solve)~
+in ~test/jacobi.t.c~ for an example wherein this method is tested on a Leslie
+matrix to recompute a given initial population distribution.
+
+In terms of accuracy, an LU factorization and back substitution approach will
+always be as correct as possible within the limits of computation; it's a
+direct solution method. It's simply the nature of the Jacobi algorithm being
+a convergent solution that determines its accuracy.
+
+LU factorization also performs in order $O(n^3)$ runtime for an $n \times n$
+matrix, whereas the Jacobi algorithm runs in order $O(k n^2) = O(n^2)$ on average
+but with the con that $k$ is given by some function on both the convergence criteria and the number of
+nonzero entries in the matrix - which might end up worse in some cases than the LU decomp approach.
+
+* Question Three
+See ~UTEST(jacobi, gauss_siedel_solve)~ in ~test/jacobi.t.c~ which runs the same
+unit test as ~UTEST(jacobi, solve_jacobi)~ but using the
+~Jacobi / Gauss-Siedel -> gauss_siedel_solve~ method as documented in the LIZFCM API reference.
+
+* Question Four, Five
+We produce the following operation counts (by hackily adding the operation count as the last element
+to the solution vector) and errors - the sum of each vector elements' absolute value away from 1.0
+using the proceeding patch and unit test.
+
+|  N | JAC opr |  JAC err | GS opr |   GS err | LU opr |   LU err |
+|  5 |    1622 | 0.001244 |    577 | 0.000098 |    430 | 0.000000 |
+|  6 |    2812 | 0.001205 |    775 | 0.000080 |    681 | 0.000000 |
+|  7 |    5396 | 0.001187 |    860 | 0.000178 |   1015 | 0.000000 |
+|  8 |    5618 | 0.001468 |   1255 | 0.000121 |   1444 | 0.000000 |
+|  9 |    7534 | 0.001638 |   1754 | 0.000091 |   1980 | 0.000000 |
+| 10 |   10342 | 0.001425 |   1847 | 0.000435 |   2635 | 0.000000 |
+| 11 |   12870 | 0.001595 |   2185 | 0.000368 |   3421 | 0.000000 |
+| 12 |   17511 | 0.001860 |   2912 | 0.000322 |   4350 | 0.000000 |
+| 13 |   16226 | 0.001631 |   3362 | 0.000270 |   5434 | 0.000000 |
+| 14 |   34333 | 0.001976 |   3844 | 0.000121 |   6685 | 0.000000 |
+| 15 |   38474 | 0.001922 |   4358 | 0.000311 |   8115 | 0.000000 |
+| 16 |   40405 | 0.002061 |   4904 | 0.000204 |   9736 | 0.000000 |
+| 17 |   58518 | 0.002125 |   5482 | 0.000311 |  11560 | 0.000000 |
+| 18 |   68079 | 0.002114 |   6092 | 0.000279 |  13599 | 0.000000 |
+| 19 |   95802 | 0.002159 |   6734 | 0.000335 |  15865 | 0.000000 |
+| 20 |   85696 | 0.002141 |   7408 | 0.000289 |  18370 | 0.000000 |
+| 21 |   89026 | 0.002316 |   8114 | 0.000393 |  21126 | 0.000000 |
+| 22 |  101537 | 0.002344 |   8852 | 0.000414 |  24145 | 0.000000 |
+| 23 |  148040 | 0.002323 |   9622 | 0.000230 |  27439 | 0.000000 |
+| 24 |  137605 | 0.002348 |  10424 | 0.000213 |  31020 | 0.000000 |
+| 25 |  169374 | 0.002409 |  11258 | 0.000894 |  34900 | 0.000000 |
+| 26 |  215166 | 0.002502 |  12124 | 0.000564 |  39091 | 0.000000 |
+| 27 |  175476 | 0.002616 |  13022 | 0.000535 |  43605 | 0.000000 |
+| 28 |  268454 | 0.002651 |  13952 | 0.000690 |  48454 | 0.000000 |
+| 29 |  267034 | 0.002697 |  14914 | 0.000675 |  53650 | 0.000000 |
+| 30 |  277193 | 0.002686 |  15908 | 0.000542 |  59205 | 0.000000 |
+| 31 |  336792 | 0.002736 |  16934 | 0.000390 |  65131 | 0.000000 |
+| 32 |  293958 | 0.002741 |  17992 | 0.000660 |  71440 | 0.000000 |
+| 33 |  323638 | 0.002893 |  19082 | 0.001072 |  78144 | 0.000000 |
+| 34 |  375104 | 0.003001 |  20204 | 0.001018 |  85255 | 0.000000 |
+| 35 |  436092 | 0.003004 |  21358 | 0.000912 |  92785 | 0.000000 |
+| 36 |  538143 | 0.003005 |  22544 | 0.000954 | 100746 | 0.000000 |
+| 37 |  511886 | 0.003029 |  23762 | 0.000462 | 109150 | 0.000000 |
+| 38 |  551332 | 0.003070 |  25012 | 0.000996 | 118009 | 0.000000 |
+| 39 |  592750 | 0.003110 |  26294 | 0.000989 | 127335 | 0.000000 |
+| 40 |  704208 | 0.003165 |  27608 | 0.000583 | 137140 | 0.000000 |
+
+#+BEGIN_SRC
+diff --git a/src/matrix.c b/src/matrix.c
+index 901a426..af5529f 100644
+--- a/src/matrix.c
++++ b/src/matrix.c
+@@ -144,20 +144,54 @@ Array_double *solve_matrix_lu_bsubst(Matrix_double *m, Array_double *b) {
+   assert(b->size == m->rows);
+   assert(m->rows == m->cols);
+
++  double opr = 0;
++
++  opr += b->size;
+   Array_double *x = copy_vector(b);
++
++  size_t n = m->rows;
++  opr += n * n;     // (u copy)
++  opr += n * n;     // l_empty
++  opr += n * n + n; // copy + put_identity_diagonal
++  opr += n;         // pivot check
++  opr += m->cols;
++  for (size_t x = 0; x < m->cols; x++) {
++    opr += (m->rows - (x + 1));
++    for (size_t y = x + 1; y < m->rows; y++) {
++      opr += 1;
++      opr += 2;     // -factor
++      opr += 4 * n; // scale, add_v, free_vector
++      opr += 1;     // -factor
++    }
++  }
++  opr += n;
+   Matrix_double **u_l = lu_decomp(m);
++
+   Matrix_double *u = u_l[0];
+   Matrix_double *l = u_l[1];
+
++  opr += n;
++  for (int64_t row = n - 1; row >= 0; row--) {
++    opr += 2 * (n - row);
++    opr += 1;
++  }
+   Array_double *b_fsub = fsubst(l, b);
++
++  opr += n;
++  for (size_t x = 0; x < n; x++) {
++    opr += 2 * (x + 1);
++    opr += 1; // /= l->data[row]->data[row]
++  }
+   x = bsubst(u, b_fsub);
+-  free_vector(b_fsub);
+
++  free_vector(b_fsub);
+   free_matrix(u);
+   free_matrix(l);
+   free(u_l);
+
+-  return x;
++  Array_double *copy = add_element(x, opr);
++  free_vector(x);
++  return copy;
+ }
+
+ Matrix_double *gaussian_elimination(Matrix_double *m) {
+@@ -231,18 +265,36 @@ Array_double *jacobi_solve(Matrix_double *m, Array_double *b,
+   assert(b->size == m->cols);
+   size_t iter = max_iterations;
+
++  double opr = 0;
++
++  opr += 2 * b->size; // to initialize two vectors with the same dim of b twice
+   Array_double *x_k = InitArrayWithSize(double, b->size, 0.0);
+   Array_double *x_k_1 =
+       InitArrayWithSize(double, b->size, rand_from(0.1, 10.0));
+
++  // add since these wouldn't be accounter for after the loop
++  opr += 1; // iter decrement
++  opr +=
++      3 * x_k_1->size; // 1 to perform x_k_1, x_k and 2 to perform ||x_k_1||_2
+   while ((--iter) > 0 && l2_distance(x_k_1, x_k) > l2_convergence_tolerance) {
++    opr += 1; // iter decrement
++    opr +=
++        3 * x_k_1->size; // 1 to perform x_k_1, x_k and 2 to perform ||x_k_1||_2
++
++    opr += m->rows; // row for add oprs
+     for (size_t i = 0; i < m->rows; i++) {
+       double delta = 0.0;
++
++      opr += m->cols;
+       for (size_t j = 0; j < m->cols; j++) {
+         if (i == j)
+           continue;
++
++        opr += 1;
+         delta += m->data[i]->data[j] * x_k->data[j];
+       }
++
++      opr += 2;
+       x_k_1->data[i] = (b->data[i] - delta) / m->data[i]->data[i];
+     }
+
+@@ -251,8 +303,9 @@ Array_double *jacobi_solve(Matrix_double *m, Array_double *b,
+     x_k_1 = tmp;
+   }
+
+-  free_vector(x_k);
+-  return x_k_1;
++  Array_double *copy = add_element(x_k_1, opr);
++  free_vector(x_k_1);
++  return copy;
+ }
+
+ Array_double *gauss_siedel_solve(Matrix_double *m, Array_double *b,
+@@ -262,30 +315,48 @@ Array_double *gauss_siedel_solve(Matrix_double *m, Array_double *b,
+   assert(b->size == m->cols);
+   size_t iter = max_iterations;
+
++  double opr = 0;
++
++  opr += 2 * b->size; // to initialize two vectors with the same dim of b twice
+   Array_double *x_k = InitArrayWithSize(double, b->size, 0.0);
+   Array_double *x_k_1 =
+       InitArrayWithSize(double, b->size, rand_from(0.1, 10.0));
+
+   while ((--iter) > 0) {
++    opr += 1; // iter decrement
++
++    opr += x_k->size; // copy oprs
+     for (size_t i = 0; i < x_k->size; i++)
+       x_k->data[i] = x_k_1->data[i];
+
++    opr += m->rows; // row for add oprs
+     for (size_t i = 0; i < m->rows; i++) {
+       double delta = 0.0;
++
++      opr += m->cols;
+       for (size_t j = 0; j < m->cols; j++) {
+         if (i == j)
+           continue;
++
++        opr += 1;
+         delta += m->data[i]->data[j] * x_k_1->data[j];
+       }
++
++      opr += 2;
+       x_k_1->data[i] = (b->data[i] - delta) / m->data[i]->data[i];
+     }
+
++    opr +=
++        3 * x_k_1->size; // 1 to perform x_k_1, x_k and 2 to perform ||x_k_1||_2
+     if (l2_distance(x_k_1, x_k) <= l2_convergence_tolerance)
+       break;
+   }
+
+   free_vector(x_k);
+-  return x_k_1;
++
++  Array_double *copy = add_element(x_k_1, opr);
++  free_vector(x_k_1);
++  return copy;
+ }
+#+END_SRC
+
+
+And this unit test:
+#+BEGIN_SRC c
+UTEST(hw_8, p4_5) {
+  printf("| N | JAC opr | JAC err | GS opr | GS err | LU opr | LU err | \n");
+
+  for (size_t i = 5; i < 100; i++) {
+    Matrix_double *m = generate_ddm(i);
+    double oprs[3] = {0.0, 0.0, 0.0};
+    double errs[3] = {0.0, 0.0, 0.0};
+
+    Array_double *b_1 = InitArrayWithSize(double, m->rows, 1.0);
+    Array_double *b = m_dot_v(m, b_1);
+    double tolerance = 0.001;
+    size_t max_iter = 400;
+
+    // JACOBI
+    {
+      Array_double *solution_with_opr_count =
+          jacobi_solve(m, b, tolerance, max_iter);
+      Array_double *solution = slice_element(solution_with_opr_count,
+                                             solution_with_opr_count->size - 1);
+
+      for (size_t i = 0; i < solution->size; i++)
+        errs[0] += fabs(solution->data[i] - 1.0);
+
+      oprs[0] =
+          solution_with_opr_count->data[solution_with_opr_count->size - 1];
+
+      free_vector(solution);
+      free_vector(solution_with_opr_count);
+    }
+
+    // GAUSS-SIEDEL
+    {
+      Array_double *solution_with_opr_count =
+          gauss_siedel_solve(m, b, tolerance, max_iter);
+      Array_double *solution = slice_element(solution_with_opr_count,
+                                             solution_with_opr_count->size - 1);
+
+      for (size_t i = 0; i < solution->size; i++)
+        errs[1] += fabs(solution->data[i] - 1.0);
+
+      oprs[1] =
+          solution_with_opr_count->data[solution_with_opr_count->size - 1];
+
+      free_vector(solution);
+      free_vector(solution_with_opr_count);
+    }
+
+    // LU-BSUBST
+    {
+      Array_double *solution_with_opr_count = solve_matrix_lu_bsubst(m, b);
+      Array_double *solution = slice_element(solution_with_opr_count,
+                                             solution_with_opr_count->size - 1);
+
+      for (size_t i = 0; i < solution->size; i++)
+        errs[2] += fabs(solution->data[i] - 1.0);
+
+      oprs[2] =
+          solution_with_opr_count->data[solution_with_opr_count->size - 1];
+
+      free_vector(solution);
+      free_vector(solution_with_opr_count);
+    }
+    free_matrix(m);
+    free_vector(b_1);
+    free_vector(b);
+
+    printf("| %zu | %f | %f | %f | %f | %f | %f | \n", i, oprs[0], errs[0],
+           oprs[1], errs[1], oprs[2], errs[2]);
+  }
+}
+#+END_SRC
diff --git a/Homework/math4610/homeworks/hw-8.pdf b/Homework/math4610/homeworks/hw-8.pdf
new file mode 100644
index 0000000..c14bb2e
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-8.pdf
diff --git a/Homework/math4610/homeworks/hw-8.tex b/Homework/math4610/homeworks/hw-8.tex
new file mode 100644
index 0000000..9071f5b
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-8.tex
@@ -0,0 +1,344 @@
+% Created 2023-12-09 Sat 22:06
+% Intended LaTeX compiler: pdflatex
+\documentclass[11pt]{article}
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+\usepackage{graphicx}
+\usepackage{longtable}
+\usepackage{wrapfig}
+\usepackage{rotating}
+\usepackage[normalem]{ulem}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{capt-of}
+\usepackage{hyperref}
+\notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+\author{Elizabeth Hunt}
+\date{\today}
+\title{Homework 8}
+\hypersetup{
+ pdfauthor={Elizabeth Hunt},
+ pdftitle={Homework 8},
+ pdfkeywords={},
+ pdfsubject={},
+ pdfcreator={Emacs 28.2 (Org mode 9.7-pre)}, 
+ pdflang={English}}
+\begin{document}
+
+\maketitle
+\setlength\parindent{0pt}
+
+\section{Question One}
+\label{sec:org800c743}
+See \texttt{UTEST(jacobi, solve\_jacobi)} in \texttt{test/jacobi.t.c} and the entry
+\texttt{Jacobi / Gauss-Siedel -> solve\_jacobi} in the LIZFCM API documentation.
+\section{Question Two}
+\label{sec:org6121bef}
+We cannot just perform the Jacobi algorithm on a Leslie matrix since
+it is obviously not diagonally dominant - which is a requirement. It is
+certainly not always the case, but, if a Leslie matrix \(L\) is invertible, we can
+first perform gaussian elimination on \(L\) augmented with \(n_{k+1}\)
+to obtain \(n_k\) with the Jacobi method. See \texttt{UTEST(jacobi, leslie\_solve)}
+in \texttt{test/jacobi.t.c} for an example wherein this method is tested on a Leslie
+matrix to recompute a given initial population distribution.
+
+In terms of accuracy, an LU factorization and back substitution approach will
+always be as correct as possible within the limits of computation; it's a
+direct solution method. It's simply the nature of the Jacobi algorithm being
+a convergent solution that determines its accuracy.
+
+LU factorization also performs in order \(O(n^3)\) runtime for an \(n \times n\)
+matrix, whereas the Jacobi algorithm runs in order \(O(k n^2) = O(n^2)\) on average
+but with the con that \(k\) is given by some function on both the convergence criteria and the number of
+nonzero entries in the matrix - which might end up worse in some cases than the LU decomp approach.
+
+\section{Question Three}
+\label{sec:org11282e6}
+See \texttt{UTEST(jacobi, gauss\_siedel\_solve)} in \texttt{test/jacobi.t.c} which runs the same
+unit test as \texttt{UTEST(jacobi, solve\_jacobi)} but using the
+\texttt{Jacobi / Gauss-Siedel -> gauss\_siedel\_solve} method as documented in the LIZFCM API reference.
+
+\section{Question Four, Five}
+\label{sec:org22b52a9}
+We produce the following operation counts (by hackily adding the operation count as the last element
+to the solution vector) and errors - the sum of each vector elements' absolute value away from 1.0
+using the proceeding patch and unit test.
+
+\begin{center}
+\begin{tabular}{rrrrrrr}
+N & JAC opr & JAC err & GS opr & GS err & LU opr & LU err\\[0pt]
+5 & 1622 & 0.001244 & 577 & 0.000098 & 430 & 0.000000\\[0pt]
+6 & 2812 & 0.001205 & 775 & 0.000080 & 681 & 0.000000\\[0pt]
+7 & 5396 & 0.001187 & 860 & 0.000178 & 1015 & 0.000000\\[0pt]
+8 & 5618 & 0.001468 & 1255 & 0.000121 & 1444 & 0.000000\\[0pt]
+9 & 7534 & 0.001638 & 1754 & 0.000091 & 1980 & 0.000000\\[0pt]
+10 & 10342 & 0.001425 & 1847 & 0.000435 & 2635 & 0.000000\\[0pt]
+11 & 12870 & 0.001595 & 2185 & 0.000368 & 3421 & 0.000000\\[0pt]
+12 & 17511 & 0.001860 & 2912 & 0.000322 & 4350 & 0.000000\\[0pt]
+13 & 16226 & 0.001631 & 3362 & 0.000270 & 5434 & 0.000000\\[0pt]
+14 & 34333 & 0.001976 & 3844 & 0.000121 & 6685 & 0.000000\\[0pt]
+15 & 38474 & 0.001922 & 4358 & 0.000311 & 8115 & 0.000000\\[0pt]
+16 & 40405 & 0.002061 & 4904 & 0.000204 & 9736 & 0.000000\\[0pt]
+17 & 58518 & 0.002125 & 5482 & 0.000311 & 11560 & 0.000000\\[0pt]
+18 & 68079 & 0.002114 & 6092 & 0.000279 & 13599 & 0.000000\\[0pt]
+19 & 95802 & 0.002159 & 6734 & 0.000335 & 15865 & 0.000000\\[0pt]
+20 & 85696 & 0.002141 & 7408 & 0.000289 & 18370 & 0.000000\\[0pt]
+21 & 89026 & 0.002316 & 8114 & 0.000393 & 21126 & 0.000000\\[0pt]
+22 & 101537 & 0.002344 & 8852 & 0.000414 & 24145 & 0.000000\\[0pt]
+23 & 148040 & 0.002323 & 9622 & 0.000230 & 27439 & 0.000000\\[0pt]
+24 & 137605 & 0.002348 & 10424 & 0.000213 & 31020 & 0.000000\\[0pt]
+25 & 169374 & 0.002409 & 11258 & 0.000894 & 34900 & 0.000000\\[0pt]
+26 & 215166 & 0.002502 & 12124 & 0.000564 & 39091 & 0.000000\\[0pt]
+27 & 175476 & 0.002616 & 13022 & 0.000535 & 43605 & 0.000000\\[0pt]
+28 & 268454 & 0.002651 & 13952 & 0.000690 & 48454 & 0.000000\\[0pt]
+29 & 267034 & 0.002697 & 14914 & 0.000675 & 53650 & 0.000000\\[0pt]
+30 & 277193 & 0.002686 & 15908 & 0.000542 & 59205 & 0.000000\\[0pt]
+31 & 336792 & 0.002736 & 16934 & 0.000390 & 65131 & 0.000000\\[0pt]
+32 & 293958 & 0.002741 & 17992 & 0.000660 & 71440 & 0.000000\\[0pt]
+33 & 323638 & 0.002893 & 19082 & 0.001072 & 78144 & 0.000000\\[0pt]
+34 & 375104 & 0.003001 & 20204 & 0.001018 & 85255 & 0.000000\\[0pt]
+35 & 436092 & 0.003004 & 21358 & 0.000912 & 92785 & 0.000000\\[0pt]
+36 & 538143 & 0.003005 & 22544 & 0.000954 & 100746 & 0.000000\\[0pt]
+37 & 511886 & 0.003029 & 23762 & 0.000462 & 109150 & 0.000000\\[0pt]
+38 & 551332 & 0.003070 & 25012 & 0.000996 & 118009 & 0.000000\\[0pt]
+39 & 592750 & 0.003110 & 26294 & 0.000989 & 127335 & 0.000000\\[0pt]
+40 & 704208 & 0.003165 & 27608 & 0.000583 & 137140 & 0.000000\\[0pt]
+\end{tabular}
+\end{center}
+
+\begin{verbatim}
+diff --git a/src/matrix.c b/src/matrix.c
+index 901a426..af5529f 100644
+--- a/src/matrix.c
++++ b/src/matrix.c
+@@ -144,20 +144,54 @@ Array_double *solve_matrix_lu_bsubst(Matrix_double *m, Array_double *b) {
+   assert(b->size == m->rows);
+   assert(m->rows == m->cols);
+
++  double opr = 0;
++
++  opr += b->size;
+   Array_double *x = copy_vector(b);
++
++  size_t n = m->rows;
++  opr += n * n;     // (u copy)
++  opr += n * n;     // l_empty
++  opr += n * n + n; // copy + put_identity_diagonal
++  opr += n;         // pivot check
++  opr += m->cols;
++  for (size_t x = 0; x < m->cols; x++) {
++    opr += (m->rows - (x + 1));
++    for (size_t y = x + 1; y < m->rows; y++) {
++      opr += 1;
++      opr += 2;     // -factor
++      opr += 4 * n; // scale, add_v, free_vector
++      opr += 1;     // -factor
++    }
++  }
++  opr += n;
+   Matrix_double **u_l = lu_decomp(m);
++
+   Matrix_double *u = u_l[0];
+   Matrix_double *l = u_l[1];
+
++  opr += n;
++  for (int64_t row = n - 1; row >= 0; row--) {
++    opr += 2 * (n - row);
++    opr += 1;
++  }
+   Array_double *b_fsub = fsubst(l, b);
++
++  opr += n;
++  for (size_t x = 0; x < n; x++) {
++    opr += 2 * (x + 1);
++    opr += 1; // /= l->data[row]->data[row]
++  }
+   x = bsubst(u, b_fsub);
+-  free_vector(b_fsub);
+
++  free_vector(b_fsub);
+   free_matrix(u);
+   free_matrix(l);
+   free(u_l);
+
+-  return x;
++  Array_double *copy = add_element(x, opr);
++  free_vector(x);
++  return copy;
+ }
+
+ Matrix_double *gaussian_elimination(Matrix_double *m) {
+@@ -231,18 +265,36 @@ Array_double *jacobi_solve(Matrix_double *m, Array_double *b,
+   assert(b->size == m->cols);
+   size_t iter = max_iterations;
+
++  double opr = 0;
++
++  opr += 2 * b->size; // to initialize two vectors with the same dim of b twice
+   Array_double *x_k = InitArrayWithSize(double, b->size, 0.0);
+   Array_double *x_k_1 =
+       InitArrayWithSize(double, b->size, rand_from(0.1, 10.0));
+
++  // add since these wouldn't be accounter for after the loop
++  opr += 1; // iter decrement
++  opr +=
++      3 * x_k_1->size; // 1 to perform x_k_1, x_k and 2 to perform ||x_k_1||_2
+   while ((--iter) > 0 && l2_distance(x_k_1, x_k) > l2_convergence_tolerance) {
++    opr += 1; // iter decrement
++    opr +=
++        3 * x_k_1->size; // 1 to perform x_k_1, x_k and 2 to perform ||x_k_1||_2
++
++    opr += m->rows; // row for add oprs
+     for (size_t i = 0; i < m->rows; i++) {
+       double delta = 0.0;
++
++      opr += m->cols;
+       for (size_t j = 0; j < m->cols; j++) {
+         if (i == j)
+           continue;
++
++        opr += 1;
+         delta += m->data[i]->data[j] * x_k->data[j];
+       }
++
++      opr += 2;
+       x_k_1->data[i] = (b->data[i] - delta) / m->data[i]->data[i];
+     }
+
+@@ -251,8 +303,9 @@ Array_double *jacobi_solve(Matrix_double *m, Array_double *b,
+     x_k_1 = tmp;
+   }
+
+-  free_vector(x_k);
+-  return x_k_1;
++  Array_double *copy = add_element(x_k_1, opr);
++  free_vector(x_k_1);
++  return copy;
+ }
+
+ Array_double *gauss_siedel_solve(Matrix_double *m, Array_double *b,
+@@ -262,30 +315,48 @@ Array_double *gauss_siedel_solve(Matrix_double *m, Array_double *b,
+   assert(b->size == m->cols);
+   size_t iter = max_iterations;
+
++  double opr = 0;
++
++  opr += 2 * b->size; // to initialize two vectors with the same dim of b twice
+   Array_double *x_k = InitArrayWithSize(double, b->size, 0.0);
+   Array_double *x_k_1 =
+       InitArrayWithSize(double, b->size, rand_from(0.1, 10.0));
+
+   while ((--iter) > 0) {
++    opr += 1; // iter decrement
++
++    opr += x_k->size; // copy oprs
+     for (size_t i = 0; i < x_k->size; i++)
+       x_k->data[i] = x_k_1->data[i];
+
++    opr += m->rows; // row for add oprs
+     for (size_t i = 0; i < m->rows; i++) {
+       double delta = 0.0;
++
++      opr += m->cols;
+       for (size_t j = 0; j < m->cols; j++) {
+         if (i == j)
+           continue;
++
++        opr += 1;
+         delta += m->data[i]->data[j] * x_k_1->data[j];
+       }
++
++      opr += 2;
+       x_k_1->data[i] = (b->data[i] - delta) / m->data[i]->data[i];
+     }
+
++    opr +=
++        3 * x_k_1->size; // 1 to perform x_k_1, x_k and 2 to perform ||x_k_1||_2
+     if (l2_distance(x_k_1, x_k) <= l2_convergence_tolerance)
+       break;
+   }
+
+   free_vector(x_k);
+-  return x_k_1;
++
++  Array_double *copy = add_element(x_k_1, opr);
++  free_vector(x_k_1);
++  return copy;
+ }
+\end{verbatim}
+
+
+And this unit test:
+\begin{verbatim}
+UTEST(hw_8, p4_5) {
+  printf("| N | JAC opr | JAC err | GS opr | GS err | LU opr | LU err | \n");
+
+  for (size_t i = 5; i < 100; i++) {
+    Matrix_double *m = generate_ddm(i);
+    double oprs[3] = {0.0, 0.0, 0.0};
+    double errs[3] = {0.0, 0.0, 0.0};
+
+    Array_double *b_1 = InitArrayWithSize(double, m->rows, 1.0);
+    Array_double *b = m_dot_v(m, b_1);
+    double tolerance = 0.001;
+    size_t max_iter = 400;
+
+    // JACOBI
+    {
+      Array_double *solution_with_opr_count =
+          jacobi_solve(m, b, tolerance, max_iter);
+      Array_double *solution = slice_element(solution_with_opr_count,
+                                             solution_with_opr_count->size - 1);
+
+      for (size_t i = 0; i < solution->size; i++)
+        errs[0] += fabs(solution->data[i] - 1.0);
+
+      oprs[0] =
+          solution_with_opr_count->data[solution_with_opr_count->size - 1];
+
+      free_vector(solution);
+      free_vector(solution_with_opr_count);
+    }
+
+    // GAUSS-SIEDEL
+    {
+      Array_double *solution_with_opr_count =
+          gauss_siedel_solve(m, b, tolerance, max_iter);
+      Array_double *solution = slice_element(solution_with_opr_count,
+                                             solution_with_opr_count->size - 1);
+
+      for (size_t i = 0; i < solution->size; i++)
+        errs[1] += fabs(solution->data[i] - 1.0);
+
+      oprs[1] =
+          solution_with_opr_count->data[solution_with_opr_count->size - 1];
+
+      free_vector(solution);
+      free_vector(solution_with_opr_count);
+    }
+
+    // LU-BSUBST
+    {
+      Array_double *solution_with_opr_count = solve_matrix_lu_bsubst(m, b);
+      Array_double *solution = slice_element(solution_with_opr_count,
+                                             solution_with_opr_count->size - 1);
+
+      for (size_t i = 0; i < solution->size; i++)
+        errs[2] += fabs(solution->data[i] - 1.0);
+
+      oprs[2] =
+          solution_with_opr_count->data[solution_with_opr_count->size - 1];
+
+      free_vector(solution);
+      free_vector(solution_with_opr_count);
+    }
+    free_matrix(m);
+    free_vector(b_1);
+    free_vector(b);
+
+    printf("| %zu | %f | %f | %f | %f | %f | %f | \n", i, oprs[0], errs[0],
+           oprs[1], errs[1], oprs[2], errs[2]);
+  }
+}
+\end{verbatim}
+\end{document}
\ No newline at end of file
diff --git a/Homework/math4610/homeworks/hw-9.org b/Homework/math4610/homeworks/hw-9.org
new file mode 100644
index 0000000..de58d2a
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-9.org
@@ -0,0 +1,222 @@
+#+TITLE: Homework 9
+#+AUTHOR: Elizabeth Hunt
+#+LATEX_HEADER: \notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+#+LATEX: \setlength\parindent{0pt}
+#+OPTIONS: toc:nil
+
+* Question One
+
+With a ~matrix_dimension~ set to 700, I consistently see about a 3x improvement in performance on my
+10-thread machine. The serial implementation gives an average ~0.189s~ total runtime, while the below
+parallel implementation runs in about ~0.066s~ after the cpu cache has filled on the first run.
+
+#+BEGIN_SRC c
+#include <math.h>
+#include <omp.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+
+#define matrix_dimension 700
+
+int n = matrix_dimension;
+float sum;
+
+int main() {
+  float A[n][n];
+  float x0[n];
+  float b[n];
+  float x1[n];
+  float res[n];
+
+  srand((unsigned int)(time(NULL)));
+
+  // not worth parallellization - rand() is not thread-safe
+  for (int i = 0; i < n; i++) {
+    for (int j = 0; j < n; j++) {
+      A[i][j] = ((float)rand() / (float)(RAND_MAX) * 5.0);
+    }
+    x0[i] = ((float)rand() / (float)(RAND_MAX) * 5.0);
+  }
+
+#pragma omp parallel for private(sum)
+  for (int i = 0; i < n; i++) {
+    sum = 0.0;
+    for (int j = 0; j < n; j++) {
+      sum += fabs(A[i][j]);
+    }
+    A[i][i] += sum;
+  }
+
+#pragma omp parallel for private(sum)
+  for (int i = 0; i < n; i++) {
+    sum = 0.0;
+    for (int j = 0; j < n; j++) {
+      sum += A[i][j];
+    }
+    b[i] = sum;
+  }
+
+  float tol = 0.0001;
+  float error = 10.0 * tol;
+  int maxiter = 100;
+  int iter = 0;
+
+  while (error > tol && iter < maxiter) {
+#pragma omp parallel for
+    for (int i = 0; i < n; i++) {
+      float temp_sum = b[i];
+      for (int j = 0; j < n; j++) {
+        temp_sum -= A[i][j] * x0[j];
+      }
+      res[i] = temp_sum;
+      x1[i] = x0[i] + res[i] / A[i][i];
+    }
+
+    sum = 0.0;
+#pragma omp parallel for reduction(+ : sum)
+    for (int i = 0; i < n; i++) {
+      float val = x1[i] - x0[i];
+      sum += val * val;
+    }
+    error = sqrt(sum);
+
+#pragma omp parallel for
+    for (int i = 0; i < n; i++) {
+      x0[i] = x1[i];
+    }
+
+    iter++;
+  }
+
+  for (int i = 0; i < n; i++)
+    printf("x[%d] = %6f \t res[%d] = %6f\n", i, x1[i], i, res[i]);
+
+  return 0;
+}
+
+#+END_SRC
+
+* Question Two
+
+I only see lowerings in performance (likely due to overhead) on my machine using OpenMP until
+~matrix_dimension~ becomes quite large, about ~300~ in testing. At ~matrix_dimension=1000~, I see another
+about 3x improvement in total runtime (including initialization & I/O which was untouched, so, even further
+improvements could be made) on my 10-thread machine; from around ~0.174~ seconds to ~.052~.
+
+#+BEGIN_SRC c
+  #include <math.h>
+  #include <stdio.h>
+  #include <stdlib.h>
+  #include <time.h>
+
+  #ifdef _OPENMP
+  #include <omp.h>
+  #else
+  #define omp_get_num_threads() 0
+  #define omp_set_num_threads(int) 0
+  #define omp_get_thread_num() 0
+  #endif
+
+  #define matrix_dimension 1000
+
+  int n = matrix_dimension;
+  float ynrm;
+
+  int main() {
+    float A[n][n];
+    float v0[n];
+    float v1[n];
+    float y[n];
+    //
+    // create a matrix
+    //
+    // not worth parallellization - rand() is not thread-safe
+    srand((unsigned int)(time(NULL)));
+    float a = 5.0;
+    for (int i = 0; i < n; i++) {
+      for (int j = 0; j < n; j++) {
+        A[i][j] = ((float)rand() / (float)(RAND_MAX)*a);
+      }
+      v0[i] = ((float)rand() / (float)(RAND_MAX)*a);
+    }
+    //
+    // modify the diagonal entries for diagonal dominance
+    // --------------------------------------------------
+    //
+    for (int i = 0; i < n; i++) {
+      float sum = 0.0;
+      for (int j = 0; j < n; j++) {
+        sum = sum + fabs(A[i][j]);
+      }
+      A[i][i] = A[i][i] + sum;
+    }
+    //
+    // generate a vector of ones
+    // -------------------------
+    //
+    for (int j = 0; j < n; j++) {
+      v0[j] = 1.0;
+    }
+    //
+    // power iteration test
+    // --------------------
+    //
+    float tol = 0.0000001;
+    float error = 10.0 * tol;
+    float lam1, lam0;
+    int maxiter = 100;
+    int iter = 0;
+
+    while (error > tol && iter < maxiter) {
+  #pragma omp parallel for
+      for (int i = 0; i < n; i++) {
+        y[i] = 0;
+        for (int j = 0; j < n; j++) {
+          y[i] = y[i] + A[i][j] * v0[j];
+        }
+      }
+
+      ynrm = 0.0;
+  #pragma omp parallel for reduction(+ : ynrm)
+      for (int i = 0; i < n; i++) {
+        ynrm += y[i] * y[i];
+      }
+      ynrm = sqrt(ynrm);
+
+  #pragma omp parallel for
+      for (int i = 0; i < n; i++) {
+        v1[i] = y[i] / ynrm;
+      }
+
+  #pragma omp parallel for
+      for (int i = 0; i < n; i++) {
+        y[i] = 0.0;
+        for (int j = 0; j < n; j++) {
+          y[i] += A[i][j] * v1[j];
+        }
+      }
+
+      lam1 = 0.0;
+  #pragma omp parallel for reduction(+ : lam1)
+      for (int i = 0; i < n; i++) {
+        lam1 += v1[i] * y[i];
+      }
+
+      error = fabs(lam1 - lam0);
+      lam0 = lam1;
+
+  #pragma omp parallel for
+      for (int i = 0; i < n; i++) {
+        v0[i] = v1[i];
+      }
+
+      iter++;
+    }
+
+    printf("in %d iterations, eigenvalue = %f\n", iter, lam1);
+  }
+#+END_SRC
+
+* Question Three
+[[https://static.simponic.xyz/lizfcm.pdf]]
diff --git a/Homework/math4610/homeworks/hw-9.pdf b/Homework/math4610/homeworks/hw-9.pdf
new file mode 100644
index 0000000..4753a3b
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-9.pdf
diff --git a/Homework/math4610/homeworks/hw-9.tex b/Homework/math4610/homeworks/hw-9.tex
new file mode 100644
index 0000000..9d5693a
--- /dev/null
+++ b/Homework/math4610/homeworks/hw-9.tex
@@ -0,0 +1,250 @@
+% Created 2023-12-11 Mon 19:24
+% Intended LaTeX compiler: pdflatex
+\documentclass[11pt]{article}
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+\usepackage{graphicx}
+\usepackage{longtable}
+\usepackage{wrapfig}
+\usepackage{rotating}
+\usepackage[normalem]{ulem}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{capt-of}
+\usepackage{hyperref}
+\notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+\author{Elizabeth Hunt}
+\date{\today}
+\title{Homework 9}
+\hypersetup{
+ pdfauthor={Elizabeth Hunt},
+ pdftitle={Homework 9},
+ pdfkeywords={},
+ pdfsubject={},
+ pdfcreator={Emacs 28.2 (Org mode 9.7-pre)}, 
+ pdflang={English}}
+\begin{document}
+
+\maketitle
+\setlength\parindent{0pt}
+
+\section{Question One}
+\label{sec:org69bed2d}
+
+With a \texttt{matrix\_dimension} set to 700, I consistently see about a 3x improvement in performance on my
+10-thread machine. The serial implementation gives an average \texttt{0.189s} total runtime, while the below
+parallel implementation runs in about \texttt{0.066s} after the cpu cache has filled on the first run.
+
+\begin{verbatim}
+#include <math.h>
+#include <omp.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+
+#define matrix_dimension 700
+
+int n = matrix_dimension;
+float sum;
+
+int main() {
+  float A[n][n];
+  float x0[n];
+  float b[n];
+  float x1[n];
+  float res[n];
+
+  srand((unsigned int)(time(NULL)));
+
+  // not worth parallellization - rand() is not thread-safe
+  for (int i = 0; i < n; i++) {
+    for (int j = 0; j < n; j++) {
+      A[i][j] = ((float)rand() / (float)(RAND_MAX) * 5.0);
+    }
+    x0[i] = ((float)rand() / (float)(RAND_MAX) * 5.0);
+  }
+
+#pragma omp parallel for private(sum)
+  for (int i = 0; i < n; i++) {
+    sum = 0.0;
+    for (int j = 0; j < n; j++) {
+      sum += fabs(A[i][j]);
+    }
+    A[i][i] += sum;
+  }
+
+#pragma omp parallel for private(sum)
+  for (int i = 0; i < n; i++) {
+    sum = 0.0;
+    for (int j = 0; j < n; j++) {
+      sum += A[i][j];
+    }
+    b[i] = sum;
+  }
+
+  float tol = 0.0001;
+  float error = 10.0 * tol;
+  int maxiter = 100;
+  int iter = 0;
+
+  while (error > tol && iter < maxiter) {
+#pragma omp parallel for
+    for (int i = 0; i < n; i++) {
+      float temp_sum = b[i];
+      for (int j = 0; j < n; j++) {
+        temp_sum -= A[i][j] * x0[j];
+      }
+      res[i] = temp_sum;
+      x1[i] = x0[i] + res[i] / A[i][i];
+    }
+
+    sum = 0.0;
+#pragma omp parallel for reduction(+ : sum)
+    for (int i = 0; i < n; i++) {
+      float val = x1[i] - x0[i];
+      sum += val * val;
+    }
+    error = sqrt(sum);
+
+#pragma omp parallel for
+    for (int i = 0; i < n; i++) {
+      x0[i] = x1[i];
+    }
+
+    iter++;
+  }
+
+  for (int i = 0; i < n; i++)
+    printf("x[%d] = %6f \t res[%d] = %6f\n", i, x1[i], i, res[i]);
+
+  return 0;
+}
+
+\end{verbatim}
+
+\section{Question Two}
+\label{sec:orgbeace25}
+
+I only see lowerings in performance (likely due to overhead) on my machine using OpenMP until
+\texttt{matrix\_dimension} becomes quite large, about \texttt{300} in testing. At \texttt{matrix\_dimension=1000}, I see another
+about 3x improvement in total runtime (including initialization \& I/O which was untouched, so, even further
+improvements could be made) on my 10-thread machine; from around \texttt{0.174} seconds to \texttt{.052}.
+
+\begin{verbatim}
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+
+#ifdef _OPENMP
+#include <omp.h>
+#else
+#define omp_get_num_threads() 0
+#define omp_set_num_threads(int) 0
+#define omp_get_thread_num() 0
+#endif
+
+#define matrix_dimension 1000
+
+int n = matrix_dimension;
+float ynrm;
+
+int main() {
+  float A[n][n];
+  float v0[n];
+  float v1[n];
+  float y[n];
+  //
+  // create a matrix
+  //
+  // not worth parallellization - rand() is not thread-safe
+  srand((unsigned int)(time(NULL)));
+  float a = 5.0;
+  for (int i = 0; i < n; i++) {
+    for (int j = 0; j < n; j++) {
+      A[i][j] = ((float)rand() / (float)(RAND_MAX)*a);
+    }
+    v0[i] = ((float)rand() / (float)(RAND_MAX)*a);
+  }
+  //
+  // modify the diagonal entries for diagonal dominance
+  // --------------------------------------------------
+  //
+  for (int i = 0; i < n; i++) {
+    float sum = 0.0;
+    for (int j = 0; j < n; j++) {
+      sum = sum + fabs(A[i][j]);
+    }
+    A[i][i] = A[i][i] + sum;
+  }
+  //
+  // generate a vector of ones
+  // -------------------------
+  //
+  for (int j = 0; j < n; j++) {
+    v0[j] = 1.0;
+  }
+  //
+  // power iteration test
+  // --------------------
+  //
+  float tol = 0.0000001;
+  float error = 10.0 * tol;
+  float lam1, lam0;
+  int maxiter = 100;
+  int iter = 0;
+
+  while (error > tol && iter < maxiter) {
+#pragma omp parallel for
+    for (int i = 0; i < n; i++) {
+      y[i] = 0;
+      for (int j = 0; j < n; j++) {
+        y[i] = y[i] + A[i][j] * v0[j];
+      }
+    }
+
+    ynrm = 0.0;
+#pragma omp parallel for reduction(+ : ynrm)
+    for (int i = 0; i < n; i++) {
+      ynrm += y[i] * y[i];
+    }
+    ynrm = sqrt(ynrm);
+
+#pragma omp parallel for
+    for (int i = 0; i < n; i++) {
+      v1[i] = y[i] / ynrm;
+    }
+
+#pragma omp parallel for
+    for (int i = 0; i < n; i++) {
+      y[i] = 0.0;
+      for (int j = 0; j < n; j++) {
+        y[i] += A[i][j] * v1[j];
+      }
+    }
+
+    lam1 = 0.0;
+#pragma omp parallel for reduction(+ : lam1)
+    for (int i = 0; i < n; i++) {
+      lam1 += v1[i] * y[i];
+    }
+
+    error = fabs(lam1 - lam0);
+    lam0 = lam1;
+
+#pragma omp parallel for
+    for (int i = 0; i < n; i++) {
+      v0[i] = v1[i];
+    }
+
+    iter++;
+  }
+
+  printf("in %d iterations, eigenvalue = %f\n", iter, lam1);
+}
+\end{verbatim}
+
+\section{Question Three}
+\label{sec:org33439e0}
+\url{https://static.simponic.xyz/lizfcm.pdf}
+\end{document}
\ No newline at end of file
diff --git a/Homework/math4610/homeworks/hw_6_p_8 b/Homework/math4610/homeworks/hw_6_p_8
new file mode 100644
index 0000000..46b58a2
--- /dev/null
+++ b/Homework/math4610/homeworks/hw_6_p_8
diff --git a/Homework/math4610/homeworks/hw_6_p_8.c b/Homework/math4610/homeworks/hw_6_p_8.c
new file mode 100644
index 0000000..56f199f
--- /dev/null
+++ b/Homework/math4610/homeworks/hw_6_p_8.c
@@ -0,0 +1,89 @@
+// compile & test w/
+//  \--> gcc -I../inc/ -Wall hw_6_p_8.c ../lib/lizfcm.a -lm -o hw_6_p_8
+//  \--> ./hw_6_p_8
+
+#include "lizfcm.h"
+#include <math.h>
+#include <stdio.h>
+
+double a(double t) {
+  double alpha = 0.1;
+  double beta = 0.001;
+  double p_0 = 2;
+  double p_infty = 29.75;
+
+  return p_0 * exp(t * (alpha - beta)) - p_infty;
+}
+
+double b(double t) {
+  double alpha = 0.1;
+  double beta = 0.001;
+  double p_0 = 2;
+  double p_infty = 115.35;
+
+  return p_0 * exp(t * (alpha - beta)) - p_infty;
+}
+
+double c(double t) {
+  double alpha = 0.1;
+  double beta = 0.0001;
+  double p_0 = 2;
+  double p_infty = 115.35;
+
+  return p_0 * exp(t * (alpha - beta)) - p_infty;
+}
+
+double d(double t) {
+  double alpha = 0.01;
+  double beta = 0.001;
+  double p_0 = 2;
+  double p_infty = 155.346;
+
+  return p_0 * exp(t * (alpha - beta)) - p_infty;
+}
+
+double e(double t) {
+  double alpha = 0.1;
+  double beta = 0.01;
+  double p_0 = 100;
+  double p_infty = 155.346;
+
+  return p_0 * exp(t * (alpha - beta)) - p_infty;
+}
+
+int main() {
+  uint64_t max_iterations = 1000;
+  double tolerance = 0.0000001;
+
+  Array_double *ivt_range = find_ivt_range(&a, -5.0, 3.0, 1000);
+  double approx_a = fixed_point_secant_bisection_method(
+      &a, ivt_range->data[0], ivt_range->data[1], tolerance, max_iterations);
+
+  free_vector(ivt_range);
+  ivt_range = find_ivt_range(&b, -5.0, 3.0, 1000);
+  double approx_b = fixed_point_secant_bisection_method(
+      &b, ivt_range->data[0], ivt_range->data[1], tolerance, max_iterations);
+
+  free_vector(ivt_range);
+  ivt_range = find_ivt_range(&c, -5.0, 3.0, 1000);
+  double approx_c = fixed_point_secant_bisection_method(
+      &c, ivt_range->data[0], ivt_range->data[1], tolerance, max_iterations);
+
+  free_vector(ivt_range);
+  ivt_range = find_ivt_range(&d, -5.0, 3.0, 1000);
+  double approx_d = fixed_point_secant_bisection_method(
+      &d, ivt_range->data[0], ivt_range->data[1], tolerance, max_iterations);
+
+  free_vector(ivt_range);
+  ivt_range = find_ivt_range(&e, -5.0, 3.0, 1000);
+  double approx_e = fixed_point_secant_bisection_method(
+      &e, ivt_range->data[0], ivt_range->data[1], tolerance, max_iterations);
+
+  printf("a ~ %f; P(%f) - P_infty = %f\n", approx_a, approx_a, a(approx_a));
+  printf("b ~ %f; P(%f) - P_infty = %f\n", approx_b, approx_b, b(approx_b));
+  printf("c ~ %f; P(%f) - P_infty = %f\n", approx_c, approx_c, c(approx_c));
+  printf("d ~ %f; P(%f) - P_infty = %f\n", approx_d, approx_d, d(approx_d));
+  printf("e ~ %f; P(%f) - P_infty = %f\n", approx_e, approx_e, e(approx_e));
+
+  return 0;
+}
diff --git a/Homework/math4610/homeworks/img/make_run.png b/Homework/math4610/homeworks/img/make_run.png
new file mode 100644
index 0000000..f20015e
--- /dev/null
+++ b/Homework/math4610/homeworks/img/make_run.png
diff --git a/Homework/math4610/homeworks/img/test_routine_1.png b/Homework/math4610/homeworks/img/test_routine_1.png
new file mode 100644
index 0000000..57ce59f
--- /dev/null
+++ b/Homework/math4610/homeworks/img/test_routine_1.png
diff --git a/Homework/math4610/homeworks/img/test_routine_2.png b/Homework/math4610/homeworks/img/test_routine_2.png
new file mode 100644
index 0000000..b116620
--- /dev/null
+++ b/Homework/math4610/homeworks/img/test_routine_2.png
diff --git a/Homework/math4610/homeworks/virtualization/hw1.pdf b/Homework/math4610/homeworks/virtualization/hw1.pdf
new file mode 100644
index 0000000..00d84a4
--- /dev/null
+++ b/Homework/math4610/homeworks/virtualization/hw1.pdf
diff --git a/Homework/math4610/homeworks/virtualization/img/htop.png b/Homework/math4610/homeworks/virtualization/img/htop.png
new file mode 100644
index 0000000..880befa
--- /dev/null
+++ b/Homework/math4610/homeworks/virtualization/img/htop.png
diff --git a/Homework/math4610/homeworks/virtualization/img/no_virtualization.png b/Homework/math4610/homeworks/virtualization/img/no_virtualization.png
new file mode 100644
index 0000000..2322456
--- /dev/null
+++ b/Homework/math4610/homeworks/virtualization/img/no_virtualization.png
diff --git a/Homework/math4610/homeworks/virtualization/virtual_machines.md b/Homework/math4610/homeworks/virtualization/virtual_machines.md
new file mode 100644
index 0000000..c6d3e12
--- /dev/null
+++ b/Homework/math4610/homeworks/virtualization/virtual_machines.md
@@ -0,0 +1,41 @@
+* Elizabeth Hunt (A02364151), MATH 4610
+
+## Virtual Machines
+
+**Question 1**
+
+Run the Linux OS as a virtual machine, or run the application in a containerized Linux environment (which
+is the same abstraction).
+
+**Question 2**
+
+A native system virtual machine has dedicated hardware to run the hypervisor, while a hosted system 
+virtual machine runs a hypervisor as a process in the operating system.
+
+**Question 3**
+
+A virtual machine hosts an entire operating system and requires users to perform configuration if they
+want to run an application, whereas a Virtual Appliance is built to provide an easy plug-and-play virtual 
+machine image built to run some specific software stack.
+
+**Question 4**
+
+In a large application sense, containerizing services into their own virtual machines allows for easier
+replication, scaling, and networking. Instead of running several smaller servers, one large server can
+host several applications in parallel. This provides a good seperation of concern. And, if one service 
+goes down, the whole system does not go down with it.
+
+Locally, it can help in development when targeting another operating system. Virtual machines can be
+used to verify builds without installing a whole other operating system.
+
+**Question 5**
+
+A virtual machine monitor is just another term for a hypervisor, so, see question 2.
+
+**Question 6**
+
+The three components of a virtual machine are:
+
+1. The host
+2. The virtualization layer
+3. The guest
diff --git a/Homework/math4610/homeworks/virtualization/virtualization.md b/Homework/math4610/homeworks/virtualization/virtualization.md
new file mode 100644
index 0000000..4d03637
--- /dev/null
+++ b/Homework/math4610/homeworks/virtualization/virtualization.md
@@ -0,0 +1,103 @@
+## Virtualization
+
+**Question 1** 
+
+I use an Apple Silicon Mac which is based on the ARM architecture - so it's necessary to use 
+[Multipass](https://multipass.run/), as native virtualization is _not available to us_.
+
+![No Virtualization Strings](./img/no_virtualization.png)
+
+**Question 2**
+
+One of the downsides of running a virtual machine, as opposed to a hosted virtual instance, is that local
+resources are used. On a laptop especially, this increases power draw, draining the battery. Additionally,
+the security of mind provided by "faster disaster recovery", as discussed in the article, is not as 
+necessary for consumer applications on personal machines as servers. Finally, virtual machines are 
+inherently slower in compute due to general overhead.
+
+**Question 3** 
+
+![htop resources](./img/htop.png)
+
+**Question 4** 
+
+In a large application sense, containerizing services into their own virtual machines allows for easier
+replication, scaling, and networking. Instead of running several smaller servers, one large server can
+host several applications in parallel.
+
+Locally, it can help in development when targeting another operating system. Virtual machines can be
+used to verify builds without installing a whole other operating system.
+
+**Question 5**
+
+A native system virtual machine has dedicated hardware to run the hypervisor, while a hosted system 
+virtual machine runs a hypervisor as a process in the operating system.
+
+**Question 6**
+
+1. Easier networking between "servers"
+2. Efficient resource use
+
+**Question 7**
+
+A Virtual Appliance is built to provide an easy plug-and-play virtual machine image built to run some 
+specific software stack.
+
+**Question 8** 
+
+A Virtual Appliance would be desirable to eliminate maintenance and configuration overhead when running an
+application. In my own experience, I've used a form of virtual appliances - "Docker Containers", to easily
+spin up multiple versions of small services at work.
+
+**Question 9** What are 2 benefits of Virtualization?
+
+See question 6.
+
+**Question 10** 
+
+See question 4.
+
+**Question 11** 
+
+See question 8.
+
+**Question 12**  What are the three main types of virtualization?
+
+1. Full virtualization
+2. Para virtualization
+3. OS-level virtualization
+
+**Question 13** What you should know about virtualization?
+
+How to create a virtual machine, and maintain it.
+
+**Question 14** What is the weakness of virtualization?
+
+Inherent overhead in all system operations.
+
+**Question 15** What are the six areas of virtualization?
+
+Source: [HiTechNectar](https://www.hitechnectar.com/blogs/virtualization-types)
+
+1. Application - run individual applications in a seperate environment than a host OS
+2. Data - abstract exact location and formatting information away from retrieval of data
+3. Desktop - hosts a desktop environment virtually on another machine (reminds me of mainframes).
+4. Network - physical networking tools are abstracted into software resources
+5. Server - division of a server into multiple guest operating systems
+6. Storage - abstraction over multiple storage arrays into a single pool
+
+**Question 16** What is the biggest challenge in virtualization?
+
+Resource distribution is a big one; it's difficult to keep track of several resources on a host machine 
+and ensure a Virtual Machine accesses them correctly.
+
+**Question 17** What is the risk of using virtualization?
+
+The biggest risk of using virtualization is sandbox escape vulnerabilities. Although mostly research and
+proof-of-concept, highly skilled engineers can theoretically craft exploits to escape the sandbox of the
+VM and directly mess with the host operating system.
+
+**Question 18**
+
+When (question 17) is trusted; sandboxing. Virtualization should supply no access to resources within the
+host operating system.