From b5ad184c1bcef78e7c2b052bcffc71f4c7381bcb Mon Sep 17 00:00:00 2001 From: Elizabeth Hunt Date: Fri, 8 Dec 2023 09:21:32 -0700 Subject: [PATCH] hw8 checkpoint --- homeworks/hw-8.org | 17 +++++++++++++++++ inc/lizfcm.h | 6 ++++++ src/matrix.c | 35 ++++++++++++++++++++++++++++++++--- src/rand.c | 7 +++++++ test/jacobi.t.c | 33 +++++++++++++++++++++++++++++++++ test/main.c | 9 ++++++++- test/rand.t.c | 10 ++++++++++ 7 files changed, 113 insertions(+), 4 deletions(-) create mode 100644 homeworks/hw-8.org create mode 100644 src/rand.c create mode 100644 test/jacobi.t.c create mode 100644 test/rand.t.c diff --git a/homeworks/hw-8.org b/homeworks/hw-8.org new file mode 100644 index 0000000..42d9dac --- /dev/null +++ b/homeworks/hw-8.org @@ -0,0 +1,17 @@ +#+TITLE: Homework 7 +#+AUTHOR: Elizabeth Hunt +#+LATEX_HEADER: \notindent \notag \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry} +#+LATEX: \setlength\parindent{0pt} +#+OPTIONS: toc:nil + +TODO: Update LIZFCM org file with jacobi solve, format_matrix_into, rand + +* Question One +See ~UTEST(jacobi, solve_jacobi)~ in ~test/jacobi.t.c~ and the entry +~Jacobi -> solve_jacobi~ in the LIZFCM API documentation. +* Question Two +A problem arises when using the Jacobi method to solve for the previous population +distribution, $n_k$, from $Ln_{k} = n_{k+1}$, because a Leslie matrix is not diagonally +dominant and will cause a division by zero. Likewise, we cannot factor it into $L$ +and $U$ terms and apply back substitution because pivot points are zero. +* Question Three diff --git a/inc/lizfcm.h b/inc/lizfcm.h index 1bb5322..70ebd2f 100644 --- a/inc/lizfcm.h +++ b/inc/lizfcm.h @@ -88,4 +88,10 @@ extern Array_double *partition_find_eigenvalues(Matrix_double *m, size_t max_iterations); extern Matrix_double *leslie_matrix(Array_double *age_class_surivor_ratio, Array_double *age_class_offspring); + +extern double rand_from(double min, double max); + +extern Array_double *jacobi_solve(Matrix_double *m, Array_double *b, + double tolerance, size_t max_iterations); + #endif // LIZFCM_H diff --git a/src/matrix.c b/src/matrix.c index 04c5adc..5f36d12 100644 --- a/src/matrix.c +++ b/src/matrix.c @@ -2,9 +2,9 @@ #include #include #include -#include +n #include -Array_double *m_dot_v(Matrix_double *m, Array_double *v) { + Array_double *m_dot_v(Matrix_double *m, Array_double *v) { assert(v->size == m->cols); Array_double *product = copy_vector(v); @@ -222,6 +222,35 @@ Array_double *solve_matrix_gaussian(Matrix_double *m, Array_double *b) { return solution; } +Array_double *jacobi_solve(Matrix_double *m, Array_double *b, + double l2_convergence_tolerance, + size_t max_iterations) { + size_t iter = max_iterations; + + Array_double *x_k = InitArrayWithSize(double, b->size, rand_from(0.1, 10.0)); + Array_double *x_k_1 = + InitArrayWithSize(double, b->size, rand_from(0.1, 10.0)); + + while ((--iter) > 0 && l2_distance(x_k_1, x_k) > l2_convergence_tolerance) { + for (size_t i = 0; i < x_k->size; i++) { + double delta = 0.0; + for (size_t j = 0; j < x_k->size; j++) { + if (i == j) + continue; + delta += m->data[i]->data[j] * x_k->data[j]; + } + x_k_1->data[i] = (b->data[i] - delta) / m->data[i]->data[i]; + } + + Array_double *tmp = x_k; + x_k = x_k_1; + x_k_1 = tmp; + } + + free_vector(x_k); + return x_k_1; +} + Matrix_double *slice_column(Matrix_double *m, size_t x) { Matrix_double *sliced = copy_matrix(m); @@ -259,7 +288,7 @@ void format_matrix_into(Matrix_double *m, char *s) { strcpy(s, "empty"); for (size_t y = 0; y < m->rows; ++y) { - char row_s[256]; + char row_s[5192]; strcpy(row_s, ""); format_vector_into(m->data[y], row_s); diff --git a/src/rand.c b/src/rand.c new file mode 100644 index 0000000..ff375ed --- /dev/null +++ b/src/rand.c @@ -0,0 +1,7 @@ +#include "lizfcm.h" + +double rand_from(double min, double max) { + double range = (max - min); + double div = RAND_MAX / range; + return min + (rand() / div); +} diff --git a/test/jacobi.t.c b/test/jacobi.t.c new file mode 100644 index 0000000..dc13d6e --- /dev/null +++ b/test/jacobi.t.c @@ -0,0 +1,33 @@ +#include "lizfcm.test.h" +#include + +Matrix_double *generate_ddm(size_t n) { + Matrix_double *m = InitMatrixWithSize(double, n, n, rand_from(0.0, 1.0)); + + for (size_t y = 0; y < m->rows; y++) { + m->data[y]->data[y] += sum_v(m->data[y]); + } + + return m; +} + +UTEST(jacobi, jacobi_solve) { + Matrix_double *m = generate_ddm(2); + + Array_double *b_1 = InitArrayWithSize(double, m->rows, 1.0); + Array_double *b = m_dot_v(m, b_1); + + double tolerance = 0.001; + size_t max_iter = 400; + Array_double *solution = jacobi_solve(m, b, tolerance, max_iter); + + for (size_t y = 0; y < m->rows; y++) { + double dot = v_dot_v(m->data[y], solution); + EXPECT_NEAR(b->data[y], dot, 0.1); + } + + free_matrix(m); + free_vector(b_1); + free_vector(b); + free_vector(solution); +} diff --git a/test/main.c b/test/main.c index 5d16fe7..6a92d4a 100644 --- a/test/main.c +++ b/test/main.c @@ -1,5 +1,12 @@ #include "lizfcm.test.h" +#include +#include UTEST(basic, unit_tests) { ASSERT_TRUE(1); } -UTEST_MAIN(); +UTEST_STATE(); +int main(int argc, const char *const argv[]) { + srand(time(NULL)); + + return utest_main(argc, argv); +} diff --git a/test/rand.t.c b/test/rand.t.c new file mode 100644 index 0000000..9ed2e9c --- /dev/null +++ b/test/rand.t.c @@ -0,0 +1,10 @@ +#include "lizfcm.test.h" + +UTEST(rand, rand_from) { + double min = -2.0; + double max = 5.0; + for (size_t i = 0; i < 1000; i++) { + double r = rand_from(min, max); + ASSERT_TRUE(min <= r && r <= max); + } +}