diff --git a/doc/software_manual.org b/doc/software_manual.org index 383b0c5..2a3b347 100644 --- a/doc/software_manual.org +++ b/doc/software_manual.org @@ -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,46 +761,51 @@ 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, - double tolerance, size_t max_iterations) { +// 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 assert(f(a) * f(b) < 0); 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 @@ -914,14 +1052,14 @@ a collection of pointers of ~Array_int~'s and its dimensions. + Output: a new ~Array_type~ with the size of the given array and its data #+BEGIN_SRC c -#define InitArray(TYPE, ...) \ - ({ \ - TYPE temp[] = __VA_ARGS__; \ - Array_##TYPE *arr = malloc(sizeof(Array_##TYPE)); \ - arr->size = sizeof(temp) / sizeof(temp[0]); \ - arr->data = malloc(arr->size * sizeof(TYPE)); \ - memcpy(arr->data, temp, arr->size * sizeof(TYPE)); \ - arr; \ +#define InitArray(TYPE, ...) \ + ({ \ + TYPE temp[] = __VA_ARGS__; \ + Array_##TYPE *arr = malloc(sizeof(Array_##TYPE)); \ + arr->size = sizeof(temp) / sizeof(temp[0]); \ + arr->data = malloc(arr->size * sizeof(TYPE)); \ + memcpy(arr->data, temp, arr->size * sizeof(TYPE)); \ + arr; \ }) #+END_SRC @@ -932,14 +1070,14 @@ a collection of pointers of ~Array_int~'s and its dimensions. + Output: a new ~Array_type~ with the given size filled with the initial value #+BEGIN_SRC c -#define InitArrayWithSize(TYPE, SIZE, INIT_VALUE) \ - ({ \ - Array_##TYPE *arr = malloc(sizeof(Array_##TYPE)); \ - arr->size = SIZE; \ - arr->data = malloc(arr->size * sizeof(TYPE)); \ - for (size_t i = 0; i < arr->size; i++) \ - arr->data[i] = INIT_VALUE; \ - arr; \ +#define InitArrayWithSize(TYPE, SIZE, INIT_VALUE) \ + ({ \ + Array_##TYPE *arr = malloc(sizeof(Array_##TYPE)); \ + arr->size = SIZE; \ + arr->data = malloc(arr->size * sizeof(TYPE)); \ + for (size_t i = 0; i < arr->size; i++) \ + arr->data[i] = INIT_VALUE; \ + arr; \ }) #+END_SRC diff --git a/doc/software_manual.pdf b/doc/software_manual.pdf index d6c03b3..4355e3f 100644 Binary files a/doc/software_manual.pdf and b/doc/software_manual.pdf differ diff --git a/doc/software_manual.tex b/doc/software_manual.tex index 4d465e6..67f86fa 100644 --- a/doc/software_manual.tex +++ b/doc/software_manual.tex @@ -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/.c -o build/.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,62 +863,66 @@ 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, - double tolerance, size_t max_iterations) { +// 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 assert(f(a) * f(b) < 0); 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\_} and \texttt{Matrix\_} -\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} @@ -1076,19 +1193,18 @@ typedef struct { \end{itemize} \begin{verbatim} -#define InitArray(TYPE, ...) \ - ({ \ - TYPE temp[] = __VA_ARGS__; \ - Array_##TYPE *arr = malloc(sizeof(Array_##TYPE)); \ - arr->size = sizeof(temp) / sizeof(temp[0]); \ - arr->data = malloc(arr->size * sizeof(TYPE)); \ - memcpy(arr->data, temp, arr->size * sizeof(TYPE)); \ - arr; \ +#define InitArray(TYPE, ...) \ + ({ \ + TYPE temp[] = __VA_ARGS__; \ + Array_##TYPE *arr = malloc(sizeof(Array_##TYPE)); \ + arr->size = sizeof(temp) / sizeof(temp[0]); \ + arr->data = malloc(arr->size * sizeof(TYPE)); \ + memcpy(arr->data, temp, arr->size * sizeof(TYPE)); \ + arr; \ }) \end{verbatim} - \item \texttt{InitArrayWithSize} -\label{sec:orga4ec165} +\label{sec:orgbe950b8} \begin{itemize} \item Author: Elizabeth Hunt \item Location: \texttt{inc/macros.h} @@ -1097,19 +1213,18 @@ typedef struct { \end{itemize} \begin{verbatim} -#define InitArrayWithSize(TYPE, SIZE, INIT_VALUE) \ - ({ \ - Array_##TYPE *arr = malloc(sizeof(Array_##TYPE)); \ - arr->size = SIZE; \ - arr->data = malloc(arr->size * sizeof(TYPE)); \ - for (size_t i = 0; i < arr->size; i++) \ - arr->data[i] = INIT_VALUE; \ - arr; \ +#define InitArrayWithSize(TYPE, SIZE, INIT_VALUE) \ + ({ \ + Array_##TYPE *arr = malloc(sizeof(Array_##TYPE)); \ + arr->size = SIZE; \ + arr->data = malloc(arr->size * sizeof(TYPE)); \ + for (size_t i = 0; i < arr->size; i++) \ + arr->data[i] = INIT_VALUE; \ + arr; \ }) \end{verbatim} - \item \texttt{InitMatrixWithSize} -\label{sec:org0748f30} +\label{sec:org5965f3b} \begin{itemize} \item Author: Elizabeth Hunt \item Location: \texttt{inc/macros.h} @@ -1131,4 +1246,4 @@ value }) \end{verbatim} \end{enumerate} -\end{document} \ No newline at end of file +\end{document} diff --git a/homeworks/a.out b/homeworks/a.out new file mode 100755 index 0000000..410a7d5 Binary files /dev/null and b/homeworks/a.out differ diff --git a/homeworks/hw-6.org b/homeworks/hw-6.org new file mode 100644 index 0000000..eebc0c2 --- /dev/null +++ b/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 +#include + +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/homeworks/hw-6.pdf b/homeworks/hw-6.pdf new file mode 100644 index 0000000..c056102 Binary files /dev/null and b/homeworks/hw-6.pdf differ diff --git a/homeworks/hw-6.tex b/homeworks/hw-6.tex new file mode 100644 index 0000000..1a0ddc4 --- /dev/null +++ b/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 +#include + +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/homeworks/hw_6_p_8 b/homeworks/hw_6_p_8 new file mode 100755 index 0000000..46b58a2 Binary files /dev/null and b/homeworks/hw_6_p_8 differ diff --git a/homeworks/hw_6_p_8.c b/homeworks/hw_6_p_8.c new file mode 100644 index 0000000..56f199f --- /dev/null +++ b/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 +#include + +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/inc/lizfcm.h b/inc/lizfcm.h index 24b7fa9..2e12a50 100644 --- a/inc/lizfcm.h +++ b/inc/lizfcm.h @@ -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, - double tolerance, size_t max_iterations); +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 diff --git a/notes/Oct-16.org b/notes/Oct-16.org index 656b261..1406737 100644 --- a/notes/Oct-16.org +++ b/notes/Oct-16.org @@ -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 diff --git a/src/approx_derivative.c b/src/approx_derivative.c index b33a208..63d0b05 100644 --- a/src/approx_derivative.c +++ b/src/approx_derivative.c @@ -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); } diff --git a/src/roots.c b/src/roots.c index de16adf..d6b22af 100644 --- a/src/roots.c +++ b/src/roots.c @@ -3,34 +3,35 @@ #include // 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, - double tolerance, size_t max_iterations) { +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 assert(f(a) * f(b) < 0); 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; } diff --git a/test/roots.t.c b/test/roots.t.c index 632832a..c20f22e 100644 --- a/test/roots.t.c +++ b/test/roots.t.c @@ -1,17 +1,114 @@ #include "lizfcm.test.h" +#include #include -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); }