simponic/lizfcm
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

hw6

Browse Source
This commit is contained in:
Elizabeth Hunt 2023-11-11 13:15:59 -07:00
parent 586d8056c1
commit 3f1f18b149
Signed by: simponic
GPG Key ID: 52B3774857EB24B1
14 changed files with 1158 additions and 202 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

188
doc/software_manual.org
View File

@ -1,4 +1,4 @@
#+TITLE: LIZFCM Software Manual (v0.2)
#+TITLE: LIZFCM Software Manual (v0.3)
#+AUTHOR: Elizabeth Hunt
#+LATEX_HEADER: \notindent \notag \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
#+LATEX: \setlength\parindent{0pt}
@ -110,8 +110,8 @@ double central_derivative_at(double (*f)(double), double a, double h) {
double x2 = a + h;
double x1 = a - h;
double y2 = (*f)(x2);
double y1 = (*f)(x1);
double y2 = f(x2);
double y1 = f(x1);
return (y2 - y1) / (x2 - x1);
}
@ -136,8 +136,8 @@ double forward_derivative_at(double (*f)(double), double a, double h) {
double x2 = a + h;
double x1 = a;
double y2 = (*f)(x2);
double y1 = (*f)(x1);
double y2 = f(x2);
double y1 = f(x1);
return (y2 - y1) / (x2 - x1);
}
@ -162,8 +162,8 @@ double backward_derivative_at(double (*f)(double), double a, double h) {
double x2 = a;
double x1 = a - h;
double y2 = (*f)(x2);
double y1 = (*f)(x1);
double y2 = f(x2);
double y1 = f(x1);
return (y2 - y1) / (x2 - x1);
}
@ -761,38 +761,42 @@ void format_matrix_into(Matrix_double *m, char *s) {
+ Input: a pointer to a oneary function taking a double and producing a double, the beginning point
in $R$ to search for a range, a ~delta~ step that is taken, and a ~max_steps~ number of maximum
iterations to perform.
+ Output: a pair of ~double~'s representing a closed closed interval ~[beginning, end]~
+ Output: a pair of ~double~'s in an ~Array_double~ representing a closed closed interval ~[beginning, end]~
#+BEGIN_SRC c
double *find_ivt_range(double (*f)(double), double start_x, double delta,
size_t max_steps) {
double *range = malloc(sizeof(double) * 2);
// f is well defined at start_x + delta*n for all n on the integer range [0,
// max_iterations]
Array_double *find_ivt_range(double (*f)(double), double start_x, double delta,
size_t max_iterations) {
double a = start_x;
while (f(a) * f(start_x) >= 0 && max_steps-- > 0)
while (f(a) * f(a + delta) >= 0 && max_iterations > 0) {
max_iterations--;
a += delta;
}
if (max_steps == 0 && f(a) * f(start_x) > 0)
double end = a + delta;
double begin = a - delta;
if (max_iterations == 0 && f(begin) * f(end) >= 0)
return NULL;
range[0] = start_x;
range[1] = a + delta;
return range;
return InitArray(double, {begin, end});
}
#+END_SRC
*** ~bisect_find_root~
+ Author: Elizabeth Hunt
+ Name(s): ~bisect_find_root~
+ Input: a one-ary function taking a double and producing a double, a closed interval represented
by ~a~ and ~b~: ~[a, b]~, a ~tolerance~ at which we return the estimated root, and a
~max_iterations~ to break us out of a loop if we can never reach the ~tolerance~
+ Output: a ~double~ representing the estimated root
by ~a~ and ~b~: ~[a, b]~, a ~tolerance~ at which we return the estimated root once $b-a < \text{tolerance}$, and a
~max_iterations~ to break us out of a loop if we can never reach the ~tolerance~.
+ Output: a vector of size of 3, ~double~'s representing first the range ~[a,b]~ and then the midpoint,
~c~ of the range.
+ Description: recursively uses binary search to split the interval until we reach ~tolerance~. We
also assume the function ~f~ is continuous on ~[a, b]~.
#+BEGIN_SRC c
double bisect_find_root(double (*f)(double), double a, double b,
// f is continuous on [a, b]
Array_double *bisect_find_root(double (*f)(double), double a, double b,
double tolerance, size_t max_iterations) {
assert(a <= b);
// guarantee there's a root somewhere between a and b by IVT
@ -800,7 +804,8 @@ double bisect_find_root(double (*f)(double), double a, double b,
double c = (1.0 / 2) * (a + b);
if (b - a < tolerance || max_iterations == 0)
return c;
return InitArray(double, {a, b, c});
if (f(a) * f(c) < 0)
return bisect_find_root(f, a, c, tolerance, max_iterations - 1);
return bisect_find_root(f, c, b, tolerance, max_iterations - 1);
@ -810,7 +815,7 @@ double bisect_find_root(double (*f)(double), double a, double b,
+ Author: Elizabeth Hunt
+ Name: ~bisect_find_root_with_error_assumption~
+ Input: a one-ary function taking a double and producing a double, a closed interval represented
by ~a~ and ~b~: ~[a, b]~, and a ~tolerance~ at which we return the estimated root
by ~a~ and ~b~: ~[a, b]~, and a ~tolerance~ equivalent to the above definition in ~bisect_find_root~
+ Output: a ~double~ representing the estimated root
+ Description: using the bisection method we know that $e_k \le (\frac{1}{2})^k (b_0 - a_0)$. So we can
calculate $k$ at the worst possible case (that the error is exactly the tolerance) to be
@ -823,7 +828,140 @@ double bisect_find_root_with_error_assumption(double (*f)(double), double a,
uint64_t max_iterations =
(uint64_t)ceil((log(tolerance) - log(b - a)) / log(1 / 2.0));
return bisect_find_root(f, a, b, tolerance, max_iterations);
Array_double *a_b_root = bisect_find_root(f, a, b, tolerance, max_iterations);
double root = a_b_root->data[2];
free_vector(a_b_root);
return root;
}
#+END_SRC
*** ~fixed_point_iteration_method~
+ Author: Elizabeth Hunt
+ Name: ~fixed_point_iteration_method~
+ Location: ~src/roots.c~
+ Input: a pointer to a oneary function $f$ taking a double and producing a double of which we are
trying to find a root, a guess $x_0$, and a function $g$ of the same signature of $f$ at which we
"step" our guesses according to the fixed point iteration method: $x_k = g(x_{k-1})$. Additionally, a
~max_iterations~ representing the maximum number of "steps" to take before arriving at our
approximation and a ~tolerance~ to return our root if it becomes within [0 - tolerance, 0 + tolerance].
+ Assumptions: $g(x)$ must be a function such that at the point $x^*$ (the found root) the derivative
$|g'(x^*)| \lt 1$
+ Output: a double representing the found approximate root $\approx x^*$.
#+BEGIN_SRC c
double fixed_point_iteration_method(double (*f)(double), double (*g)(double),
double x_0, double tolerance,
size_t max_iterations) {
if (max_iterations <= 0)
return x_0;
double root = g(x_0);
if (tolerance >= fabs(f(root)))
return root;
return fixed_point_iteration_method(f, g, root, tolerance,
max_iterations - 1);
}
#+END_SRC
*** ~fixed_point_newton_method~
+ Author: Elizabeth Hunt
+ Name: ~fixed_point_newton_method~
+ Location: ~src/roots.c~
+ Input: a pointer to a oneary function $f$ taking a double and producing a double of which we are
trying to find a root and another pointer to a function fprime of the same signature, a guess $x_0$,
and a ~max_iterations~ and ~tolerance~ as defined in the above method are required inputs.
+ Description: continually computes elements in the sequence $x_n = x_{n-1} - \frac{f(x_{n-1})}{f'p(x_{n-1})}$
+ Output: a double representing the found approximate root $\approx x^*$ recursively applied to the sequence
given
#+BEGIN_SRC c
double fixed_point_newton_method(double (*f)(double), double (*fprime)(double),
double x_0, double tolerance,
size_t max_iterations) {
if (max_iterations <= 0)
return x_0;
double root = x_0 - f(x_0) / fprime(x_0);
if (tolerance >= fabs(f(root)))
return root;
return fixed_point_newton_method(f, fprime, root, tolerance,
max_iterations - 1);
}
#+END_SRC
*** ~fixed_point_secant_method~
+ Author: Elizabeth Hunt
+ Name: ~fixed_point_secant_method~
+ Location: ~src/roots.c~
+ Input: a pointer to a oneary function $f$ taking a double and producing a double of which we are
trying to find a root, a guess $x_0$ and $x_1$ in which a root lies between $[x_0, x_1]$; applying the
sequence $x_n = x_{n-1} - f(x_{n-1}) \frac{x_{n-1} - x_{n-2}}{f(x_{n-1}) - f(x_{n-2})}$.
Additionally, a ~max_iterations~ and ~tolerance~ as defined in the above method are required
inputs.
+ Output: a double representing the found approximate root $\approx x^*$ recursively applied to the sequence.
#+BEGIN_SRC c
double fixed_point_secant_method(double (*f)(double), double x_0, double x_1,
double tolerance, size_t max_iterations) {
if (max_iterations == 0)
return x_1;
double root = x_1 - f(x_1) * ((x_1 - x_0) / (f(x_1) - f(x_0)));
if (tolerance >= fabs(f(root)))
return root;
return fixed_point_secant_method(f, x_1, root, tolerance, max_iterations - 1);
}
#+END_SRC
*** ~fixed_point_secant_bisection_method~
+ Author: Elizabeth Hunt
+ Name: ~fixed_point_secant_method~
+ Location: ~src/roots.c~
+ Input: a pointer to a oneary function $f$ taking a double and producing a double of which we are
trying to find a root, a guess $x_0$, and a $x_1$ of which we define our first interval $[x_0, x_1]$.
Then, we perform a single iteration of the ~fixed_point_secant_method~ on this interval; if it
produces a root outside, we refresh the interval and root respectively with the given
~bisect_find_root~ method. Additionally, a ~max_iterations~ and ~tolerance~ as defined in the above method are required
inputs.
+ Output: a double representing the found approximate root $\approx x^*$ continually applied with the
constraints defined.
#+BEGIN_SRC c
double fixed_point_secant_bisection_method(double (*f)(double), double x_0,
double x_1, double tolerance,
size_t max_iterations) {
double begin = x_0;
double end = x_1;
double root = x_0;
while (tolerance < fabs(f(root)) && max_iterations > 0) {
max_iterations--;
double secant_root = fixed_point_secant_method(f, begin, end, tolerance, 1);
if (secant_root < begin || secant_root > end) {
Array_double *range_root = bisect_find_root(f, begin, end, tolerance, 1);
begin = range_root->data[0];
end = range_root->data[1];
root = range_root->data[2];
free_vector(range_root);
continue;
}
root = secant_root;
if (f(root) * f(begin) < 0)
end = secant_root; // the root exists in [begin, secant_root]
else
begin = secant_root;
}
return root;
}
#+END_SRC

BIN
doc/software_manual.pdf
View File

Binary file not shown.

333
doc/software_manual.tex
View File

@ -1,4 +1,4 @@
% Created 2023-11-01 Wed 20:52
% Created 2023-11-10 Fri 20:54
% Intended LaTeX compiler: pdflatex
\documentclass[11pt]{article}
\usepackage[utf8]{inputenc}
@ -15,13 +15,13 @@
\notindent \notag \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
\author{Elizabeth Hunt}
\date{\today}
\title{LIZFCM Software Manual (v0.2)}
\title{LIZFCM Software Manual (v0.3)}
\hypersetup{
pdfauthor={Elizabeth Hunt},
pdftitle={LIZFCM Software Manual (v0.2)},
pdftitle={LIZFCM Software Manual (v0.3)},
pdfkeywords={},
pdfsubject={},
pdfcreator={Emacs 28.2 (Org mode 9.7-pre)},
pdfcreator={Emacs 29.1 (Org mode 9.7-pre)},
pdflang={English}}
\begin{document}
@ -29,9 +29,8 @@
\tableofcontents
\setlength\parindent{0pt}
\section{Design}
\label{sec:org9458aa0}
\label{sec:orgdac8577}
The LIZFCM static library (at \url{https://github.com/Simponic/math-4610}) is a successor to my
attempt at writing codes for the Fundamentals of Computational Mathematics course in Common
Lisp, but the effort required to meet the requirement of creating a static library became
@ -47,9 +46,8 @@ the C programming language. I have a couple tenets for its design:
\item Routines are separated into "modules" that follow a form of separation of concerns
in files, and not individual files per function.
\end{itemize}
\section{Compilation}
\label{sec:orge0bab70}
\label{sec:org7755023}
A provided \texttt{Makefile} is added for convencience. It has been tested on an \texttt{arm}-based M1 machine running
MacOS as well as \texttt{x86} Arch Linux.
@ -71,15 +69,14 @@ produce an object file:
gcc -Iinc/ -lm -Wall -c src/<the_routine>.c -o build/<the_routine>.o
\end{verbatim}
Which is then bundled into a static library in \texttt{lib/lizfcm.a} which can be linked
Which is then bundled into a static library in \texttt{lib/lizfcm.a} and can be linked
in the standard method.
\section{The LIZFCM API}
\label{sec:org91f4707}
\label{sec:org940357c}
\subsection{Simple Routines}
\label{sec:orgc8c57e4}
\label{sec:org28486b0}
\subsubsection{\texttt{smaceps}}
\label{sec:orgfeb6ef6}
\label{sec:org1de3a4e}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{smaceps}
@ -103,9 +100,8 @@ float smaceps() {
return machine_epsilon;
}
\end{verbatim}
\subsubsection{\texttt{dmaceps}}
\label{sec:orgb3dc0f2}
\label{sec:org742e61e}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{dmaceps}
@ -129,11 +125,10 @@ double dmaceps() {
return machine_epsilon;
}
\end{verbatim}
\subsection{Derivative Routines}
\label{sec:orge88d677}
\label{sec:org21233d3}
\subsubsection{\texttt{central\_derivative\_at}}
\label{sec:org32a8384}
\label{sec:org6a00f6c}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{central\_derivative\_at}
@ -162,9 +157,8 @@ double central_derivative_at(double (*f)(double), double a, double h) {
return (y2 - y1) / (x2 - x1);
}
\end{verbatim}
\subsubsection{\texttt{forward\_derivative\_at}}
\label{sec:orgb6fdb9a}
\label{sec:org78f40a9}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{forward\_derivative\_at}
@ -193,9 +187,8 @@ double forward_derivative_at(double (*f)(double), double a, double h) {
return (y2 - y1) / (x2 - x1);
}
\end{verbatim}
\subsubsection{\texttt{backward\_derivative\_at}}
\label{sec:org8b6070e}
\label{sec:org888d29e}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{backward\_derivative\_at}
@ -224,11 +217,10 @@ double backward_derivative_at(double (*f)(double), double a, double h) {
return (y2 - y1) / (x2 - x1);
}
\end{verbatim}
\subsection{Vector Routines}
\label{sec:org161e049}
\label{sec:org73b87ea}
\subsubsection{Vector Arithmetic: \texttt{add\_v, minus\_v}}
\label{sec:org938756a}
\label{sec:orgf8b5da1}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name(s): \texttt{add\_v}, \texttt{minus\_v}
@ -257,9 +249,8 @@ Array_double *minus_v(Array_double *v1, Array_double *v2) {
return sub;
}
\end{verbatim}
\subsubsection{Norms: \texttt{l1\_norm}, \texttt{l2\_norm}, \texttt{linf\_norm}}
\label{sec:org53e3d42}
\label{sec:orgc5368a1}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name(s): \texttt{l1\_norm}, \texttt{l2\_norm}, \texttt{linf\_norm}
@ -291,9 +282,8 @@ double linf_norm(Array_double *v) {
return max;
}
\end{verbatim}
\subsubsection{\texttt{vector\_distance}}
\label{sec:org31d6d43}
\label{sec:org0505e0b}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{vector\_distance}
@ -312,9 +302,8 @@ double vector_distance(Array_double *v1, Array_double *v2,
return dist;
}
\end{verbatim}
\subsubsection{Distances: \texttt{l1\_distance}, \texttt{l2\_distance}, \texttt{linf\_distance}}
\label{sec:org3c2cede}
\label{sec:org1c45dae}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name(s): \texttt{l1\_distance}, \texttt{l2\_distance}, \texttt{linf\_distance}
@ -338,9 +327,8 @@ double linf_distance(Array_double *v1, Array_double *v2) {
return vector_distance(v1, v2, &linf_norm);
}
\end{verbatim}
\subsubsection{\texttt{sum\_v}}
\label{sec:orgde8ccf4}
\label{sec:org687d1bd}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{sum\_v}
@ -357,10 +345,8 @@ double sum_v(Array_double *v) {
return sum;
}
\end{verbatim}
\subsubsection{\texttt{scale\_v}}
\label{sec:orgb6465fa}
\label{sec:org5926df1}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{scale\_v}
@ -377,9 +363,8 @@ Array_double *scale_v(Array_double *v, double m) {
return copy;
}
\end{verbatim}
\subsubsection{\texttt{free\_vector}}
\label{sec:org38c1352}
\label{sec:org3458f6a}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{free\_vector}
@ -395,9 +380,8 @@ void free_vector(Array_double *v) {
free(v);
}
\end{verbatim}
\subsubsection{\texttt{add\_element}}
\label{sec:org9fa4fc9}
\label{sec:org54cba50}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{add\_element}
@ -415,9 +399,8 @@ Array_double *add_element(Array_double *v, double x) {
return pushed;
}
\end{verbatim}
\subsubsection{\texttt{slice\_element}}
\label{sec:orga743fd5}
\label{sec:org02cd40a}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{slice\_element}
@ -434,9 +417,8 @@ Array_double *slice_element(Array_double *v, size_t x) {
return sliced;
}
\end{verbatim}
\subsubsection{\texttt{copy\_vector}}
\label{sec:org8918aa7}
\label{sec:org4b0c599}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{copy\_vector}
@ -454,9 +436,8 @@ Array_double *copy_vector(Array_double *v) {
return copy;
}
\end{verbatim}
\subsubsection{\texttt{format\_vector\_into}}
\label{sec:org744df1b}
\label{sec:orgde12441}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{format\_vector\_into}
@ -484,11 +465,10 @@ void format_vector_into(Array_double *v, char *s) {
strcat(s, "\n");
}
\end{verbatim}
\subsection{Matrix Routines}
\label{sec:orge1c8a5a}
\label{sec:orgd85d8ec}
\subsubsection{\texttt{lu\_decomp}}
\label{sec:org19cc6a1}
\label{sec:org6a14cbd}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{lu\_decomp}
@ -548,7 +528,7 @@ Matrix_double **lu_decomp(Matrix_double *m) {
}
\end{verbatim}
\subsubsection{\texttt{bsubst}}
\label{sec:org786580f}
\label{sec:org8b51171}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{bsubst}
@ -573,7 +553,7 @@ Array_double *bsubst(Matrix_double *u, Array_double *b) {
}
\end{verbatim}
\subsubsection{\texttt{fsubst}}
\label{sec:org1d422c6}
\label{sec:orgf9180a0}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{fsubst}
@ -599,9 +579,8 @@ Array_double *fsubst(Matrix_double *l, Array_double *b) {
return x;
}
\end{verbatim}
\subsubsection{\texttt{solve\_matrix\_lu\_bsubst}}
\label{sec:orgbf1dbcb}
\label{sec:orgf3845f4}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{src/matrix.c}
@ -636,9 +615,8 @@ Array_double *solve_matrix_lu_bsubst(Matrix_double *m, Array_double *b) {
return x;
}
\end{verbatim}
\subsubsection{\texttt{gaussian\_elimination}}
\label{sec:orgc3ceb7b}
\label{sec:orge926b79}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{src/matrix.c}
@ -692,9 +670,8 @@ Matrix_double *gaussian_elimination(Matrix_double *m) {
return m_cp;
}
\end{verbatim}
\subsubsection{\texttt{solve\_matrix\_gaussian}}
\label{sec:orgb8fc210}
\label{sec:orgc4f0d99}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{src/matrix.c}
@ -726,10 +703,8 @@ Array_double *solve_matrix_gaussian(Matrix_double *m, Array_double *b) {
return solution;
}
\end{verbatim}
\subsubsection{\texttt{m\_dot\_v}}
\label{sec:org304f5e5}
\label{sec:orgb7015af}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{src/matrix.c}
@ -749,9 +724,8 @@ Array_double *m_dot_v(Matrix_double *m, Array_double *v) {
return product;
}
\end{verbatim}
\subsubsection{\texttt{put\_identity\_diagonal}}
\label{sec:orga145f39}
\label{sec:orge955396}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{src/matrix.c}
@ -768,9 +742,8 @@ Matrix_double *put_identity_diagonal(Matrix_double *m) {
return copy;
}
\end{verbatim}
\subsubsection{\texttt{slice\_column}}
\label{sec:org1ea6d1a}
\label{sec:org886997f}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{src/matrix.c}
@ -792,9 +765,8 @@ Matrix_double *slice_column(Matrix_double *m, size_t x) {
return sliced;
}
\end{verbatim}
\subsubsection{\texttt{add\_column}}
\label{sec:org733cc61}
\label{sec:org405e1c5}
\begin{itemize}
\item Author: Elizabet Hunt
\item Location: \texttt{src/matrix.c}
@ -816,9 +788,8 @@ Matrix_double *add_column(Matrix_double *m, Array_double *v) {
return pushed;
}
\end{verbatim}
\subsubsection{\texttt{copy\_matrix}}
\label{sec:orge8936ce}
\label{sec:org01ea984}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{src/matrix.c}
@ -836,9 +807,8 @@ Matrix_double *copy_matrix(Matrix_double *m) {
return copy;
}
\end{verbatim}
\subsubsection{\texttt{free\_matrix}}
\label{sec:orgf7b674e}
\label{sec:orgab8c2cf}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{src/matrix.c}
@ -855,9 +825,8 @@ void free_matrix(Matrix_double *m) {
free(m);
}
\end{verbatim}
\subsubsection{\texttt{format\_matrix\_into}}
\label{sec:org22902bd}
\label{sec:org9e01978}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{format\_matrix\_into}
@ -884,9 +853,9 @@ void format_matrix_into(Matrix_double *m, char *s) {
}
\end{verbatim}
\subsection{Root Finding Methods}
\label{sec:org6c22e6c}
\label{sec:org81f315b}
\subsubsection{\texttt{find\_ivt\_range}}
\label{sec:org43ba5e5}
\label{sec:orgc1dde4d}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{find\_ivt\_range}
@ -894,42 +863,45 @@ void format_matrix_into(Matrix_double *m, char *s) {
\item Input: a pointer to a oneary function taking a double and producing a double, the beginning point
in \(R\) to search for a range, a \texttt{delta} step that is taken, and a \texttt{max\_steps} number of maximum
iterations to perform.
\item Output: a pair of \texttt{double}'s representing a closed closed interval \texttt{[beginning, end]}
\item Output: a pair of \texttt{double}'s in an \texttt{Array\_double} representing a closed closed interval \texttt{[beginning, end]}
\end{itemize}
\begin{verbatim}
double *find_ivt_range(double (*f)(double), double start_x, double delta,
size_t max_steps) {
double *range = malloc(sizeof(double) * 2);
// f is well defined at start_x + delta*n for all n on the integer range [0,
// max_iterations]
Array_double *find_ivt_range(double (*f)(double), double start_x, double delta,
size_t max_iterations) {
double a = start_x;
while (f(a) * f(start_x) >= 0 && max_steps-- > 0)
while (f(a) * f(a + delta) >= 0 && max_iterations > 0) {
max_iterations--;
a += delta;
}
if (max_steps == 0 && f(a) * f(start_x) > 0)
double end = a + delta;
double begin = a - delta;
if (max_iterations == 0 && f(begin) * f(end) >= 0)
return NULL;
range[0] = start_x;
range[1] = a + delta;
return range;
return InitArray(double, {begin, end});
}
\end{verbatim}
\subsubsection{\texttt{bisect\_find\_root}}
\label{sec:orgf8a3f0e}
\label{sec:orgb42a836}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name(s): \texttt{bisect\_find\_root}
\item Input: a one-ary function taking a double and producing a double, a closed interval represented
by \texttt{a} and \texttt{b}: \texttt{[a, b]}, a \texttt{tolerance} at which we return the estimated root, and a
\texttt{max\_iterations} to break us out of a loop if we can never reach the \texttt{tolerance}
\item Output: a \texttt{double} representing the estimated root
by \texttt{a} and \texttt{b}: \texttt{[a, b]}, a \texttt{tolerance} at which we return the estimated root once \(b-a < \text{tolerance}\), and a
\texttt{max\_iterations} to break us out of a loop if we can never reach the \texttt{tolerance}.
\item Output: a vector of size of 3 \texttt{double}'s representing first the .
\item Description: recursively uses binary search to split the interval until we reach \texttt{tolerance}. We
also assume the function \texttt{f} is continuous on \texttt{[a, b]}.
\end{itemize}
\begin{verbatim}
double bisect_find_root(double (*f)(double), double a, double b,
// f is continuous on [a, b]
Array_double *bisect_find_root(double (*f)(double), double a, double b,
double tolerance, size_t max_iterations) {
assert(a <= b);
// guarantee there's a root somewhere between a and b by IVT
@ -937,19 +909,20 @@ double bisect_find_root(double (*f)(double), double a, double b,
double c = (1.0 / 2) * (a + b);
if (b - a < tolerance || max_iterations == 0)
return c;
return InitArray(double, a, b, c);
if (f(a) * f(c) < 0)
return bisect_find_root(f, a, c, tolerance, max_iterations - 1);
return bisect_find_root(f, c, b, tolerance, max_iterations - 1);
}
\end{verbatim}
\subsubsection{\texttt{bisect\_find\_root\_with\_error\_assumption}}
\label{sec:orgeb72b17}
\label{sec:org762134e}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{bisect\_find\_root\_with\_error\_assumption}
\item Input: a one-ary function taking a double and producing a double, a closed interval represented
by \texttt{a} and \texttt{b}: \texttt{[a, b]}, and a \texttt{tolerance} at which we return the estimated root
by \texttt{a} and \texttt{b}: \texttt{[a, b]}, and a \texttt{tolerance} equivalent to the above definition in \texttt{bisect\_find\_root}
\item Output: a \texttt{double} representing the estimated root
\item Description: using the bisection method we know that \(e_k \le (\frac{1}{2})^k (b_0 - a_0)\). So we can
calculate \(k\) at the worst possible case (that the error is exactly the tolerance) to be
@ -963,14 +936,160 @@ double bisect_find_root_with_error_assumption(double (*f)(double), double a,
uint64_t max_iterations =
(uint64_t)ceil((log(tolerance) - log(b - a)) / log(1 / 2.0));
return bisect_find_root(f, a, b, tolerance, max_iterations);
Array_double *a_b_root = bisect_find_root(f, a, b, tolerance, max_iterations);
double root = a_b_root->data[2];
free_vector(a_b_root);
return root;
}
\end{verbatim}
\subsubsection{\texttt{fixed\_point\_iteration\_method}}
\label{sec:org9f210ad}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{fixed\_point\_iteration\_method}
\item Location: \texttt{src/roots.c}
\item Input: a pointer to a oneary function \(f\) taking a double and producing a double of which we are
trying to find a root, a guess \(x_0\), and a function \(g\) of the same signature of \(f\) at which we
"step" our guesses according to the fixed point iteration method: \(x_k = g(x_{k-1})\). Additionally, a
\texttt{max\_iterations} representing the maximum number of "steps" to take before arriving at our
approximation and a \texttt{tolerance} to return our root if it becomes within [0 - tolerance, 0 + tolerance].
\item Assumptions: \(g(x)\) must be a function such that at the point \(x^*\) (the found root) the derivative
\(|g'(x^*)| \lt 1\)
\item Output: a double representing the found approximate root \(\approx x^*\).
\end{itemize}
\begin{verbatim}
double fixed_point_iteration_method(double (*f)(double), double (*g)(double),
double x_0, double tolerance,
size_t max_iterations) {
if (max_iterations <= 0)
return x_0;
double root = g(x_0);
if (tolerance >= fabs(f(root)))
return root;
return fixed_point_iteration_method(f, g, root, tolerance,
max_iterations - 1);
}
\end{verbatim}
\subsubsection{\texttt{fixed\_point\_newton\_method}}
\label{sec:orgedecc45}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{fixed\_point\_newton\_method}
\item Location: \texttt{src/roots.c}
\item Input: a pointer to a oneary function \(f\) taking a double and producing a double of which we are
trying to find a root and another pointer to a function fprime of the same signature, a guess \(x_0\),
and a \texttt{max\_iterations} and \texttt{tolerance} as defined in the above method are required inputs.
\item Description: continually computes elements in the sequence \(x_n = x_{n-1} - \frac{f(x_{n-1})}{f'p(x_{n-1})}\)
\item Output: a double representing the found approximate root \(\approx x^*\) recursively applied to the sequence
given
\end{itemize}
\begin{verbatim}
double fixed_point_newton_method(double (*f)(double), double (*fprime)(double),
double x_0, double tolerance,
size_t max_iterations) {
if (max_iterations <= 0)
return x_0;
double root = x_0 - f(x_0) / fprime(x_0);
if (tolerance >= fabs(f(root)))
return root;
return fixed_point_newton_method(f, fprime, root, tolerance,
max_iterations - 1);
}
\end{verbatim}
\subsubsection{\texttt{fixed\_point\_secant\_method}}
\label{sec:org63bcbe2}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{fixed\_point\_secant\_method}
\item Location: \texttt{src/roots.c}
\item Input: a pointer to a oneary function \(f\) taking a double and producing a double of which we are
trying to find a root, a guess \(x_0\), and a \(\delta\) of our first guess at which we draw the first
secant line according to the sequence \(x_n = x_{n-1} - f(x_{n-1}) \frac{x_{n-1} - x_{n-2}}{f(x_{n-1}) - f(x_{n-2})}\) which
thus simplifies to \(x_1 = (x_0 + \delta) - f(x_0 + \delta) \frac{(x_0 + \delta) - x_0}{f(x_0 + \delta) - f(x_0)} = (x_0 + \delta) - f(x_0 + \delta) \frac{\delta}{f(x_0 + \delta) - f(x_0)}\).
Additionally, a \texttt{max\_iterations} and \texttt{tolerance} as defined in the above method are required
inputs.
\item Output: a double representing the found approximate root \(\approx x^*\) recursively applied to the sequence.
\end{itemize}
\begin{verbatim}
double fixed_point_secant_method(double (*f)(double), double x_0, double delta,
double tolerance, size_t max_iterations) {
if (max_iterations <= 0)
return x_0;
double x_1 = x_0 + delta;
double root = x_1 - f(x_1) * (delta / (f(x_1) - f(x_0)));
if (tolerance >= fabs(f(root)))
return root;
double new_delta = root - x_1;
return fixed_point_secant_method(f, x_1, new_delta, tolerance,
max_iterations);
}
\end{verbatim}
\subsubsection{\texttt{fixed\_point\_secant\_bisection\_method}}
\label{sec:org72d3074}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{fixed\_point\_secant\_method}
\item Location: \texttt{src/roots.c}
\item Input: a pointer to a oneary function \(f\) taking a double and producing a double of which we are
trying to find a root, a guess \(x_0\), and a \(\delta\) of which we define our first interval \([x_0, x_0 + \delta]\).
Then, we perform a single iteration of the \texttt{fixed\_point\_secant\_method} on this interval; if it
produces a root outside, we refresh the interval and root respectively with the given
\texttt{bisect\_find\_root} method. Additionally, a \texttt{max\_iterations} and \texttt{tolerance} as defined in the above method are required
inputs.
\item Output: a double representing the found approximate root \(\approx x^*\) continually applied with the
constraints defined.
\end{itemize}
\begin{verbatim}
double fixed_point_secant_bisection_method(double (*f)(double), double x_0,
double delta, double tolerance,
size_t max_iterations) {
double begin = x_0;
double end = x_0 + delta;
double root = x_0;
while (tolerance < fabs(f(root)) && max_iterations > 0) {
max_iterations--;
double secant_root =
fixed_point_secant_method(f, begin, end - begin, tolerance, 1);
if (secant_root < begin || secant_root > end) {
Array_double *range_root = bisect_find_root(f, begin, end, tolerance, 1);
begin = range_root->data[0];
end = range_root->data[1];
root = range_root->data[2];
free_vector(range_root);
continue;
}
root = secant_root;
// the root exists in [begin, secant_root]
if (f(root) * f(begin) < 0)
end = secant_root;
else
begin = secant_root;
}
return root;
}
\end{verbatim}
\subsection{Linear Routines}
\label{sec:org4e14ee5}
\label{sec:org04f3e56}
\subsubsection{\texttt{least\_squares\_lin\_reg}}
\label{sec:orge0ed136}
\label{sec:orgbd48d8e}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Name: \texttt{least\_squares\_lin\_reg}
@ -1000,12 +1119,12 @@ Line *least_squares_lin_reg(Array_double *x, Array_double *y) {
}
\end{verbatim}
\subsection{Appendix / Miscellaneous}
\label{sec:org0130d70}
\label{sec:orgf6b30a5}
\subsubsection{Data Types}
\label{sec:org8aa1c01}
\label{sec:orgd382789}
\begin{enumerate}
\item \texttt{Line}
\label{sec:org596b0e7}
\label{sec:orgab590b9}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{inc/types.h}
@ -1018,7 +1137,7 @@ typedef struct Line {
} Line;
\end{verbatim}
\item The \texttt{Array\_<type>} and \texttt{Matrix\_<type>}
\label{sec:org9d1c7c3}
\label{sec:org5be3024}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{inc/types.h}
@ -1048,12 +1167,11 @@ typedef struct {
} Matrix_int
\end{verbatim}
\end{enumerate}
\subsubsection{Macros}
\label{sec:orgb835bfa}
\label{sec:org20a391c}
\begin{enumerate}
\item \texttt{c\_max} and \texttt{c\_min}
\label{sec:org9ca763b}
\label{sec:orgfc6117a}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{inc/macros.h}
@ -1065,9 +1183,8 @@ typedef struct {
#define c_max(x, y) (((x) >= (y)) ? (x) : (y))
#define c_min(x, y) (((x) <= (y)) ? (x) : (y))
\end{verbatim}
\item \texttt{InitArray}
\label{sec:org3454dab}
\label{sec:org472f039}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{inc/macros.h}
@ -1086,9 +1203,8 @@ typedef struct {
arr; \
})
\end{verbatim}
\item \texttt{InitArrayWithSize}
\label{sec:orga4ec165}
\label{sec:orgbe950b8}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{inc/macros.h}
@ -1107,9 +1223,8 @@ typedef struct {
arr; \
})
\end{verbatim}
\item \texttt{InitMatrixWithSize}
\label{sec:org0748f30}
\label{sec:org5965f3b}
\begin{itemize}
\item Author: Elizabeth Hunt
\item Location: \texttt{inc/macros.h}

BIN
homeworks/a.out Executable file
View File

Binary file not shown.

199
homeworks/hw-6.org Normal file
View File

@ -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

BIN
homeworks/hw-6.pdf Normal file
View File

Binary file not shown.

223
homeworks/hw-6.tex Normal file
View File

@ -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}

BIN
homeworks/hw_6_p_8 Executable file
View File

Binary file not shown.

89
homeworks/hw_6_p_8.c Normal file
View File

@ -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;
}

22
inc/lizfcm.h
View File

@ -50,12 +50,26 @@ extern int matrix_equal(Matrix_double *a, Matrix_double *b);
extern Line *least_squares_lin_reg(Array_double *x, Array_double *y);
extern double *find_ivt_range(double (*f)(double), double start_x, double delta,
size_t max_steps);
extern double bisect_find_root(double (*f)(double), double a, double b,
extern Array_double *find_ivt_range(double (*f)(double), double start_x,
double delta, size_t max_steps);
extern Array_double *bisect_find_root(double (*f)(double), double a, double b,
double tolerance, size_t max_iterations);
extern double bisect_find_root_with_error_assumption(double (*f)(double),
double a, double b,
double tolerance);
extern double fixed_point_iteration_method(double (*f)(double),
double (*g)(double), double x_0,
double tolerance,
size_t max_iterations);
extern double fixed_point_newton_method(double (*f)(double),
double (*fprime)(double), double x_0,
double tolerance,
size_t max_iterations);
extern double fixed_point_secant_method(double (*f)(double), double x_0,
double x_1, double tolerance,
size_t max_iterations);
extern double fixed_point_secant_bisection_method(double (*f)(double),
double x_0, double x_1,
double tolerance,
size_t max_iterations);
#endif // LIZFCM_H

2
notes/Oct-16.org
View File

@ -42,7 +42,7 @@ Then x^* = 0
If we construct g(x) = 10x + xe^-x
Then g'(x) = 10 + (e^-x - xe^-x) \Rightarrow g'(x) = 10 + e^0 - 0 = 11 (this wouldn't converge)n
Then g'(x) = 10 + (e^-x - xe^-x) \Rightarrow g'(x) = 10 + e^0 - 0 = 11 (this wouldn't converge)
However if g(x)) = x - (xe^-x), g'(x) = 1 - (e^-x - xe^-x) \Rightarrow g'(x^*) = 0

12
src/approx_derivative.c
View File

@ -7,8 +7,8 @@ double central_derivative_at(double (*f)(double), double a, double h) {
double x2 = a + h;
double x1 = a - h;
double y2 = (*f)(x2);
double y1 = (*f)(x1);
double y2 = f(x2);
double y1 = f(x1);
return (y2 - y1) / (x2 - x1);
}
@ -19,8 +19,8 @@ double forward_derivative_at(double (*f)(double), double a, double h) {
double x2 = a + h;
double x1 = a;
double y2 = (*f)(x2);
double y1 = (*f)(x1);
double y2 = f(x2);
double y1 = f(x1);
return (y2 - y1) / (x2 - x1);
}
@ -31,8 +31,8 @@ double backward_derivative_at(double (*f)(double), double a, double h) {
double x2 = a;
double x1 = a - h;
double y2 = (*f)(x2);
double y1 = (*f)(x1);
double y2 = f(x2);
double y1 = f(x1);
return (y2 - y1) / (x2 - x1);
}

109
src/roots.c
View File

@ -3,26 +3,26 @@
#include <math.h>
// f is well defined at start_x + delta*n for all n on the integer range [0,
// max_steps]
double *find_ivt_range(double (*f)(double), double start_x, double delta,
size_t max_steps) {
double *range = malloc(sizeof(double) * 2);
// max_iterations]
Array_double *find_ivt_range(double (*f)(double), double start_x, double delta,
size_t max_iterations) {
double a = start_x;
while (f(a) * f(start_x) >= 0 && max_steps-- > 0)
while (f(a) * f(a + delta) >= 0 && max_iterations > 0) {
max_iterations--;
a += delta;
}
if (max_steps == 0 && f(a) * f(start_x) > 0)
double end = a + delta;
double begin = a - delta;
if (max_iterations == 0 && f(begin) * f(end) >= 0)
return NULL;
range[0] = start_x;
range[1] = a + delta;
return range;
return InitArray(double, {begin, end});
}
// f is continuous on [a, b]
double bisect_find_root(double (*f)(double), double a, double b,
Array_double *bisect_find_root(double (*f)(double), double a, double b,
double tolerance, size_t max_iterations) {
assert(a <= b);
// guarantee there's a root somewhere between a and b by IVT
@ -30,7 +30,8 @@ double bisect_find_root(double (*f)(double), double a, double b,
double c = (1.0 / 2) * (a + b);
if (b - a < tolerance || max_iterations == 0)
return c;
return InitArray(double, {a, b, c});
if (f(a) * f(c) < 0)
return bisect_find_root(f, a, c, tolerance, max_iterations - 1);
return bisect_find_root(f, c, b, tolerance, max_iterations - 1);
@ -42,5 +43,85 @@ double bisect_find_root_with_error_assumption(double (*f)(double), double a,
uint64_t max_iterations =
(uint64_t)ceil((log(tolerance) - log(b - a)) / log(1 / 2.0));
return bisect_find_root(f, a, b, tolerance, max_iterations);
Array_double *a_b_root = bisect_find_root(f, a, b, tolerance, max_iterations);
double root = a_b_root->data[2];
free_vector(a_b_root);
return root;
}
double fixed_point_iteration_method(double (*f)(double), double (*g)(double),
double x_0, double tolerance,
size_t max_iterations) {
if (max_iterations <= 0)
return x_0;
double root = g(x_0);
if (tolerance >= fabs(f(root)))
return root;
return fixed_point_iteration_method(f, g, root, tolerance,
max_iterations - 1);
}
double fixed_point_newton_method(double (*f)(double), double (*fprime)(double),
double x_0, double tolerance,
size_t max_iterations) {
if (max_iterations <= 0)
return x_0;
double root = x_0 - f(x_0) / fprime(x_0);
if (tolerance >= fabs(f(root)))
return root;
return fixed_point_newton_method(f, fprime, root, tolerance,
max_iterations - 1);
}
double fixed_point_secant_method(double (*f)(double), double x_0, double x_1,
double tolerance, size_t max_iterations) {
if (max_iterations == 0)
return x_1;
double root = x_1 - f(x_1) * ((x_1 - x_0) / (f(x_1) - f(x_0)));
if (tolerance >= fabs(f(root)))
return root;
return fixed_point_secant_method(f, x_1, root, tolerance, max_iterations - 1);
}
double fixed_point_secant_bisection_method(double (*f)(double), double x_0,
double x_1, double tolerance,
size_t max_iterations) {
double begin = x_0;
double end = x_1;
double root = x_0;
while (tolerance < fabs(f(root)) && max_iterations > 0) {
max_iterations--;
double secant_root = fixed_point_secant_method(f, begin, end, tolerance, 1);
if (secant_root < begin || secant_root > end) {
Array_double *range_root = bisect_find_root(f, begin, end, tolerance, 1);
begin = range_root->data[0];
end = range_root->data[1];
root = range_root->data[2];
free_vector(range_root);
continue;
}
root = secant_root;
if (f(root) * f(begin) < 0)
end = secant_root; // the root exists in [begin, secant_root]
else
begin = secant_root;
}
return root;
}

109
test/roots.t.c
View File

@ -1,17 +1,114 @@
#include "lizfcm.test.h"
#include <math.h>
#include <stdio.h>
double g(double x) { return x * x - 9; }
double f1(double x) { return x * x - 9; }
double f2(double x) { return x * x - 5 * x + 6; }
double f2prime(double x) { return 2 * x - 5; }
double g1(double x) { return x + f2(x); }
double g2(double x) { return x - f2(x); }
UTEST(ivt, find_interval) {
double *range = find_ivt_range(&g, -100.0, 1.0, 200);
EXPECT_LT(g(range[0]) * g(range[1]), 0);
Array_double *range = find_ivt_range(&f1, -10.0, 0.10, 200);
EXPECT_LT(f1(range->data[0]) * f1(range->data[1]), 0);
free(range);
free_vector(range);
}
UTEST(root, bisection_with_error_assumption) {
double root = bisect_find_root_with_error_assumption(&g, -5, 0, 0.01);
Array_double *range = find_ivt_range(&f2, 2.5, 0.10, 200);
EXPECT_NEAR(-3, root, 0.01);
double tolerance = 0.01;
double root1 = bisect_find_root_with_error_assumption(
&f2, range->data[0], range->data[1], tolerance);
free_vector(range);
range = find_ivt_range(&f2, 0, 0.01, 500);
double root2 = bisect_find_root_with_error_assumption(
&f2, range->data[0], range->data[1], tolerance);
free_vector(range);
EXPECT_NEAR(3.0, root1, tolerance);
EXPECT_NEAR(2.0, root2, tolerance);
}
UTEST(root, fixed_point_iteration_method) {
// x^2 - 5x + 6 = (x - 3)(x - 2)
double expect_x2 = 3.0;
double expect_x1 = 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_x1, 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_x2, tolerance);
}
UTEST(root, fixed_point_newton_method) {
// x^2 - 5x + 6 = (x - 3)(x - 2)
double expect_x2 = 3.0;
double expect_x1 = 2.0;
double tolerance = 0.01;
uint64_t max_iterations = 10;
double x_0 = 1.55; // 1.5 < 1.55 < 2.5
double root1 =
fixed_point_newton_method(&f2, &f2prime, x_0, tolerance, max_iterations);
EXPECT_NEAR(root1, expect_x1, tolerance);
x_0 = 3.4; // 2.5 < 3.4 < 3.5
double root2 =
fixed_point_newton_method(&f2, &f2prime, x_0, tolerance, max_iterations);
EXPECT_NEAR(root2, expect_x2, tolerance);
}
UTEST(root, fixed_point_secant_method) {
// x^2 - 5x + 6 = (x - 3)(x - 2)
double expect_x2 = 3.0;
double expect_x1 = 2.0;
double delta = 0.01;
double tolerance = 0.01;
uint64_t max_iterations = 10;
double x_0 = 1.55; // 1.5 < 1.55 < 2.5
double root1 = fixed_point_secant_method(&f2, x_0, x_0 + delta, tolerance,
max_iterations);
EXPECT_NEAR(root1, expect_x1, tolerance);
x_0 = 3.4; // 2.5 < 3.4 < 3.5
double root2 = fixed_point_secant_method(&f2, x_0, x_0 + delta, tolerance,
max_iterations);
EXPECT_NEAR(root2, expect_x2, tolerance);
}
UTEST(root, fixed_point_hybrid_method) {
// x^2 - 5x + 6 = (x - 3)(x - 2)
double expect_x2 = 3.0;
double expect_x1 = 2.0;
double delta = 1.0;
double tolerance = 0.01;
uint64_t max_iterations = 10;
double x_0 = 1.55;
double root1 = fixed_point_secant_bisection_method(&f2, x_0, x_0 + delta,
tolerance, max_iterations);
EXPECT_NEAR(root1, expect_x1, tolerance);
x_0 = 2.5;
double root2 = fixed_point_secant_bisection_method(&f2, x_0, x_0 + delta,
tolerance, max_iterations);
EXPECT_NEAR(root2, expect_x2, tolerance);
}