add least dominant eigenvalue

This commit is contained in:
Elizabeth Hunt 2023-11-27 13:27:49 -07:00
parent f50a1c5a4d
commit 2a35f68ac4
Signed by: simponic
GPG Key ID: 52B3774857EB24B1
3 changed files with 62 additions and 0 deletions

View File

@ -76,6 +76,9 @@ extern double fixed_point_secant_bisection_method(double (*f)(double),
extern double dominant_eigenvalue(Matrix_double *m, Array_double *v, extern double dominant_eigenvalue(Matrix_double *m, Array_double *v,
double tolerance, size_t max_iterations); double tolerance, size_t max_iterations);
extern double least_dominant_eigenvalue(Matrix_double *m, Array_double *v,
double tolerance,
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);
#endif // LIZFCM_H #endif // LIZFCM_H

View File

@ -48,3 +48,38 @@ double dominant_eigenvalue(Matrix_double *m, Array_double *v, double tolerance,
return lambda; return lambda;
} }
double least_dominant_eigenvalue(Matrix_double *m, Array_double *v,
double tolerance, size_t max_iterations) {
assert(m->rows == m->cols);
assert(m->rows == v->size);
double shift = 0.0;
Matrix_double *m_c = copy_matrix(m);
for (size_t y = 0; y < m_c->rows; ++y)
m_c->data[y]->data[y] = m_c->data[y]->data[y] - shift;
double error = tolerance;
size_t iter = max_iterations;
double lambda = shift;
Array_double *eigenvector_1 = copy_vector(v);
while (error >= tolerance && (--iter) > 0) {
Array_double *eigenvector_2 = solve_matrix_lu_bsubst(m_c, eigenvector_1);
Array_double *normalized_eigenvector_2 =
scale_v(eigenvector_2, 1.0 / linf_norm(eigenvector_2));
free_vector(eigenvector_2);
eigenvector_2 = normalized_eigenvector_2;
Array_double *mx = m_dot_v(m, eigenvector_2);
double new_lambda =
v_dot_v(mx, eigenvector_2) / v_dot_v(eigenvector_2, eigenvector_2);
error = fabs(new_lambda - lambda);
lambda = new_lambda;
free_vector(eigenvector_1);
eigenvector_1 = eigenvector_2;
}
return lambda;
}

View File

@ -43,6 +43,30 @@ UTEST(eigen, leslie_matrix_dominant_eigenvalue) {
free_matrix(leslie); free_matrix(leslie);
} }
UTEST(eigen, least_dominant_eigenvalue) {
Matrix_double *m = InitMatrixWithSize(double, 3, 3, 0.0);
m->data[0]->data[0] = 2.0;
m->data[0]->data[1] = 2.0;
m->data[0]->data[2] = 4.0;
m->data[1]->data[0] = 1.0;
m->data[1]->data[1] = 4.0;
m->data[1]->data[2] = 7.0;
m->data[2]->data[1] = 2.0;
m->data[2]->data[2] = 6.0;
double expected_least_dominant_eigenvalue = 0.87689; // 5 - sqrt(17)
double tolerance = 0.0001;
uint64_t max_iterations = 64;
Array_double *v_guess = InitArrayWithSize(double, 3, 2.0);
double approx_least_dominant_eigenvalue =
least_dominant_eigenvalue(m, v_guess, tolerance, max_iterations);
EXPECT_NEAR(expected_least_dominant_eigenvalue,
approx_least_dominant_eigenvalue, tolerance);
}
UTEST(eigen, dominant_eigenvalue) { UTEST(eigen, dominant_eigenvalue) {
Matrix_double *m = InitMatrixWithSize(double, 2, 2, 0.0); Matrix_double *m = InitMatrixWithSize(double, 2, 2, 0.0);
m->data[0]->data[0] = 2.0; m->data[0]->data[0] = 2.0;