Для решения системы линейных алгебраических уравнений (ЛАУ) вида 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;
}