hw8 checkpoint

This commit is contained in:
Elizabeth Hunt 2023-12-08 09:21:32 -07:00
parent b8c662456f
commit b5ad184c1b
Signed by: simponic
GPG Key ID: 52B3774857EB24B1
7 changed files with 113 additions and 4 deletions

17
homeworks/hw-8.org Normal file
View File

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

View File

@ -88,4 +88,10 @@ extern Array_double *partition_find_eigenvalues(Matrix_double *m,
size_t max_iterations); size_t max_iterations);
extern Matrix_double *leslie_matrix(Array_double *age_class_surivor_ratio, extern Matrix_double *leslie_matrix(Array_double *age_class_surivor_ratio,
Array_double *age_class_offspring); 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 #endif // LIZFCM_H

View File

@ -2,9 +2,9 @@
#include <assert.h> #include <assert.h>
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
#include <string.h> n #include<string.h>
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); assert(v->size == m->cols);
Array_double *product = copy_vector(v); Array_double *product = copy_vector(v);
@ -222,6 +222,35 @@ Array_double *solve_matrix_gaussian(Matrix_double *m, Array_double *b) {
return solution; 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 *slice_column(Matrix_double *m, size_t x) {
Matrix_double *sliced = copy_matrix(m); Matrix_double *sliced = copy_matrix(m);
@ -259,7 +288,7 @@ void format_matrix_into(Matrix_double *m, char *s) {
strcpy(s, "empty"); strcpy(s, "empty");
for (size_t y = 0; y < m->rows; ++y) { for (size_t y = 0; y < m->rows; ++y) {
char row_s[256]; char row_s[5192];
strcpy(row_s, ""); strcpy(row_s, "");
format_vector_into(m->data[y], row_s); format_vector_into(m->data[y], row_s);

7
src/rand.c Normal file
View File

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

33
test/jacobi.t.c Normal file
View File

@ -0,0 +1,33 @@
#include "lizfcm.test.h"
#include <math.h>
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);
}

View File

@ -1,5 +1,12 @@
#include "lizfcm.test.h" #include "lizfcm.test.h"
#include <stdlib.h>
#include <time.h>
UTEST(basic, unit_tests) { ASSERT_TRUE(1); } 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);
}

10
test/rand.t.c Normal file
View File

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