Здесь описан алгоритм решения системы линейных уравнений с помощью так называемого метода Гаусса. Программу вы можете скачать разделе программы. Алгоритм реализован на языке С.

Пусть у нас есть система N линейных уравнений

a₁₁x₁ + a₁₂x₂ + a₁₃x₃ + ... a_1Nx_N = b₁
a₂₁x₁ + a₂₂x₂ + a₂₃x₃ + ... a_2Nx_N = b₂
a₃₁x₁ + a₃₂x₂ + a₃₃x₃ + ... a_3Nx_N = b₃
...
a_N1x₁ + a_N2x₂ + a_N3x₃ + ... a_NNx_N = b_N

где x_i - неизвестные, a_ij - коэффициенты при неизвестных, b_i - свободные члены в уравнениях, i,j пробегают значения от 1 до N.

Цель задачи - зная a_ij и b_i найти x_i.

Суть метода Гаусса состоит в том, что с помощью некоторых операций исходную систему уравнений можно свести к более простой системе. Эта простая система имеет треугольный вид:

a11x1 +	a12x2 +	a13x3 +	...	a1NxN = b1
	a22x2 +	a23x3 +	...	a2NxN = b2
		a33x3 +	...	a3NxN = b3
...
			...	aNNxN = bN

Особенность этой системы - в строках с номером i все коэффициенты a_ij при j<i равны нулю.

Если мы смогли привести нашу систему уравнений к такому треугольному виду, то решить уравнения уже просто. Из последнего уравнения находим x_N= b_N / a_NN. Дальше подставляем его в предпоследнее уравнение и находим из него x_N-1. Подставляем оба найденных решения в следующее с конца уравнение и находим x_N-2. И так далее, пока не найдем x₁, на чем решение заканчивается. Такая процедура называется обратной прогонкой.

Теперь перейдем к вопросу как же добиться того, чтобы система стала треугольной.

Из линейной алгебры (см. например, Крутицкая Н.И., Тихонравов А.В., Шишкин А.А., Аналитическая геометрия и линейная алгебра с приложениями) известно что если к некоторой строке системы уравнений прибавить любую линейную комбинацию любых других строк этой системы, то решение системы не изменится. Под линейной комбинацией строк понимается сумма строк, каждая из которых у множается на некоторое число (в принципе, любое).

Нужно, чтобы во второй строке получилось уравнение, в которой отсутствует член при x₁. Прибавим к этой строке первую строку, умноженную на некоторое число M.

(a₁₁x₁ + a₁₂x₂ + a₁₃x₃ + ... a_1Nx_N = b₁)*M +
a₂₁x₁ + a₂₂x₂ + a₂₃x₃ + ... a_2Nx_N = b₂

Получим

(a₁₁*М + a₂₁) x₁ + ... = b₁*M + b₂

Для того, чтобы член при x₁ равнялся нулю, нужно, чтобы M = - a₂₁ / a₁₁. Проделав эту операцию, получившееся уравнение запишем вместо второго и приступим к третьему уравнению. К нему мы прибавим первое уравнение, умноженное на M = - a₃₁ / a₁₁ и тоже получим ноль вместо члена при x₁. Такую операцию нужно проделать над всеми остальными уравнениями. В результате получим систему такого вида:

a11x1 +	a12x2 +	a13x3 +	...	a1NxN = b1
	a22x2 +	a23x3 +	...	a2NxN = b2
	a32x2 +	a33x3 +	...	a3NxN = b3
...
	aN2x2 +	aN3x3 +	...	aNNxN = bN

После этого будем избавляться от членов при x₂ в третьем, четвертом, N-ом уравнении. Для этого нужно к уравнению с j-м номером прибавить 2-ое уравнение, умноженное на M = - a_j2 / a₂₂. Проделав эту операцию над всеми остальными уравнениями, получим систему где нет членов с x₂ в уравнениях с номером больше 2.

И так далее... Проделав это для третьего члена, четвертого... до тех пор, пока не кончатся уравнения, получим в итоге систему треугольного вида.

Ниже приведен текст программы.

//HOWTO solve system of linear equations
//ax=b
//where a[N][N] is square matrix,
//b[N] and x[N] are columns

//data are stored in input text file:

//input file format:
//a_11 a_12 a_13... a_1N b_1
//a_21 a_22 a_23... a_2N b_2
//...
//a_N1 a_N2 a_N3... a_NN b_N

#include <stdio.h>
#include <process.h>
float **a, *b, *x;
int N;
char filename[256];
FILE* InFile=NULL;
void count_num_lines(){
   //count number of lines in input file - number of equations
   int nelf=0;       //non empty line flag
   do{
       nelf = 0;
       while(fgetc(InFile)!='\n' && !feof(InFile)) nelf=1;
       if(nelf) N++;
   }while(!feof(InFile));
}

void freematrix(){
   //free memory for matrixes
   int i;
   for(i=0; i<N; i++){
       delete [] a[i];
   }
   delete [] a;
   delete [] b;
   delete [] x;
}

void allocmatrix(){
   //allocate memory for matrixes
   int i,j;
   x = new float[N];
   b = new float[N];
   a = new float*[N];
   if(x==NULL || b==NULL || a==NULL){
       printf("\nNot enough memory to allocate for %d equations.\n", N);
       exit(-1);
   }
   for(i=0; i<N; i++){
       a[i] = new float[N];
       if(a[i]==NULL){
	   printf("\nNot enough memory to allocate for %d equations.\n", N);
       }
   }
   for(i=0; i<N; i++){
       for(j=0; j<N; j++){
	   a[i][j]=0;
       }
       b[i]=0;
       x[i]=0;
   }
}

void readmatrix(){
   int i=0,j=0;
   //read matrixes a and b from input file
   for(i=0; i<N; i++){
       for(j=0; j<N; j++){
	   fscanf(InFile, "%f", &a[i][j]);
       }
       fscanf(InFile, "%f", &b[i]);
   }
}

void printmatrix(){
   //print matrix "a"
   int i=0,j=0;
   printf("\n");
   for(i=0; i<N; i++){
       for(j=0; j<N; j++){
	   printf(" %+f*X%d", a[i][j], j);
       }
       printf(" =%f\n", b[i]);
   }
}

void testsolve(){
   //test that ax=b
   int i=0,j=0;
   printf("\n");
   for(i=0; i<N; i++){
       float s = 0;
       for(j=0; j<N; j++){
	   s += a[i][j]*x[j];
       }
       printf("%f\t%f\n", s, b[i]);
   }
}
void printresult(){
   int i=0;
   printf("\n");
   printf("Result\n");
   for(i=0; i<N; i++){
       printf("X%d = %f\n", i, x[i]);
   }
}
void diagonal(){
   int i, j, k;
   float temp=0;
   for(i=0; i<N; i++){
       if(a[i][i]==0){
	   for(j=0; j<N; j++){
	       if(j==i) continue;
	       if(a[j][i] !=0 && a[i][j]!=0){
		   for(k=0; k<N; k++){
		       temp = a[j][k];
		       a[j][k] = a[i][k];
		       a[i][k] = temp;
		   }
		   temp = b[j];
		   b[j] = b[i];
		   b[i] = temp;
		   break;
	       }
	   }
       }
   }
}
void cls(){
   for(int i=0; i<25; i++) printf("\n");
}
void main(){
   int i=0,j=0, k=0;
   cls();
   do{
       printf("\nInput filename: ");
       scanf("%s", filename);
       InFile = fopen(filename, "rt");
   }while(InFile==NULL);
   count_num_lines();
   allocmatrix();
   rewind(InFile);
   //read data from file
   readmatrix();
   //check if there are 0 on main diagonal and exchange rows in that case
   diagonal();
   fclose(InFile);
   printmatrix();
   //process rows
   for(k=0; k<N; k++){
       for(i=k+1; i<N; i++){
	   if(a[k][k]==0){
	       printf("\nSolution is not exist.\n");
	       return;
	   }
	   float M = a[i][k] / a[k][k];
	   for(j=k; j<N; j++){
	       a[i][j] -= M * a[k][j];
	   }
	   b[i] -= M*b[k];
       }
   }
   printmatrix();
   for(i=N-1; i>=0; i--){
       float s = 0;
       for(j = i; j<N; j++){
	   s = s+a[i][j]*x[j];
       }
       x[i] = (b[i] - s) / a[i][i];
   }

   InFile = fopen(filename, "rt");
   readmatrix();
   fclose(InFile);
   printmatrix();
   testsolve();
   printresult();
   freematrix();
}

Поскольку весь метод уже описан, коротко опишу, как эта программа получает уравнения на вход, что выдает в результате решения.

Чтобы воспользоваться этой программой, нужно запустить скомпилированный исполняемый файл. В первую очередь программа спросит, откуда ей брать коэффициенты для уравнений. Создайте в любом текстовом редакторе (но только не в Word-е а, например в notepad-е) файл, где напишите коэффициенты уравнений построчно через пробел, приблизительно так:

a11	a12	a13	...	a1N	b1
a21	a22	a23	...	a2N	b2
a31	a32	a33	...	a3N	b3
aN1	aN2	aN3	...	aNN	bN

Например:

2	-3	1	5
-1	-1	2	0
1	2	-1	3

Этот файл необходимо создать в той директории, где лежит программа, иначе она его не найдет. В результате работы программы, она выдаст результат, нечто вроде:

X0 = 3.000000
X1 = 1.000000
X2 = 2.000000

Это и есть решение системы уравнений, т.е. набор неизвестных Х.

Коротко опишем, для чего служит такое большое количество подпрограмм в данной программе:

• void count_num_lines() - подсчитывает количество уравнений в системе
• void allocmatrix() - выделяет память для массивов для хранения коэффициентов уравнений, свободных членов и решения
• void readmatrix() - прочитывает из файла коэффициенты и свободные члены в массивы
• void printmatrix() - распечатывает систему уравнений
• void diagonal() - делает так, чтобы на главной диагонали не было нулей, чтобы не пришлось однажды делить на ноль в процессе приведения матрицы к треугольному виду
• void testsolve() - подставляет решение в систему и распечатывает рядом то, что получается в левой части уравнения и для сравнения печатает рядом свободные члены, те, что в правой части уравнения; получившиеся два столбца должны совпадать, если уравнения решены правильно
• void printresult() - распечатывает получившийся столбец решений
• void freematrix() - освобождает память, которая была выделена ранее
• void cls() - стирает экран в начале работы программы
• void main() - основная функция из которой последовательно вызываются все вышеперечисленные функции, и проходит процесс приведения системы уравнений к треугольному виду и обратная прогонка.

Программа снабжена комментариями, что облегчает ее понимание.

Надеюсь, что вам стало все понятно.