Для решения системы линейных алгебраических уравнений (ЛАУ) вида AX=B,
где A - прямоугольная матрица размерностью m строк на n столбцов,
m>n, применение находит метод, основанный на SVD-разложении (Singular
Value Decomposition) матрицы A:
По найденному разложению матрицы A решение системы определяется
следующим образом:
Исходные тексты на языке СИ программ, реализующих SVD-разложение и решение на его базе систем ЛАУ, можно найти в литературе и в составе математических пакетов (например, http://www.netlib.org/clapack или http://alglib.sources.ru/matrixops/general).
Ниже приводится пример использования для решения системы ЛАУ функций из библиотеки GSL (GNU Scientific Library).
/* * Решение системы ЛАУ A*X=B (кол-во уравнений > кол-во неизвестных) методом SVD * (Singular Value Decomposition) * * Компиляция: gcc -o prg prg.c -L/usr/local/lib -lgsl -lgslcblas -lm */ #include <stdio.h> #include <gsl/gsl_linalg.h> int main (int argc, char *argv[]) { /* Коэффициенты матрицы A[30x4] */ double a_data[] = { 1.0, 0.1, 0.01, 0.001, 1.0, 0.3, 0.09, 0.027, 1.0, 0.5, 0.25, 0.125, 1.0, 0.7, 0.49, 0.343, 1.0, 0.9, 0.81, 0.729, 1.0, 1.1, 1.21, 1.331, 1.0, 1.3, 1.69, 2.197, 1.0, 1.5, 2.25, 3.375, 1.0, 1.7, 2.89, 4.913, 1.0, 1.9, 3.61, 6.859, 1.0, 2.1, 4.41, 9.261, 1.0, 2.3, 5.29, 12.167, 1.0, 2.5, 6.25, 15.625, 1.0, 2.7, 7.29, 19.683, 1.0, 2.9, 8.41, 24.389, 1.0, 3.1, 9.61, 29.791, 1.0, 3.3, 10.89, 35.937, 1.0, 3.5, 12.25, 42.875, 1.0, 3.7, 13.69, 50.653, 1.0, 3.9, 15.21, 59.319, 1.0, 4.1, 16.81, 68.921, 1.0, 4.3, 18.49, 79.507, 1.0, 4.5, 20.25, 91.125, 1.0, 4.7, 22.09, 103.823, 1.0, 4.9, 24.01, 117.649, 1.0, 5.1, 26.01, 132.651, 1.0, 5.3, 28.09, 148.877, 1.0, 5.5, 30.25, 166.375, 1.0, 5.7, 32.49, 185.193, 1.0, 5.9, 34.81, 205.379 }; /* Элементы вектора правых частей B[30]*/ double b_data[] = { 0.0998334, 0.29552, 0.479426, 0.644218, 0.783327, 0.891207, 0.963558, 0.997495, 0.991665, 0.9463, 0.863209, 0.745705, 0.598472, 0.42738, 0.239249, 0.0415807, -0.157746, -0.350783, -0.529836, -0.687766, -0.818277, -0.916166, -0.97753, -0.999923, -0.982453, -0.925815, -0.832267, -0.70554, -0.550686, -0.373877 }; /* Создание "представления" (view) для матрицы коэффициентов A из массива */ gsl_matrix_view a = gsl_matrix_view_array (a_data, 30, 4); /* Теперь доступ к матрице через a.matrix */ /* Создание "представления" для вектора B из массива */ gsl_vector_view b = gsl_vector_view_array (b_data, 30); /* Теперь доступ к вектору через b.vector */ gsl_matrix *v = gsl_matrix_alloc (4,4); /* Вспомогательная матрица */ gsl_vector *s = gsl_vector_alloc (4); /* Вектор сингулярных значений */ gsl_vector *work = gsl_vector_alloc (4); /* Рабочий вектор */ gsl_vector *x = gsl_vector_alloc (4); /* Вектор X */ /* Для доступа к элементам вектора можно использовать след. функции * double gsl_vector_get (const gsl_vector * v, size_t i) * void gsl_vector_set (gsl_vector * v, size_t i, double x) * * Для доступа к элементам матрицы можно использовать след. функции * double gsl_matrix_get (const gsl_matrix * m, size_t i, size_t j) * void gsl_matrix_set (gsl_matrix * m, size_t i, size_t j, double x) * * Индексация элементов - с 0 */ /* SVD - разложение матрицы A */ gsl_linalg_SV_decomp (&a.matrix, v, s, work); /* Для коэффициентов A изменены значениями !!! */ /* Решение системы ЛАУ */ gsl_linalg_SV_solve (&a.matrix, v, s, &b.vector, x); /* Печать результатов */ printf ("x = \n"); gsl_vector_fprintf (stdout, x, "%g"); gsl_vector_free (s); gsl_matrix_free (v); return 0; }