Метод Гаусса для вирішення систем лінійних рівнянь
На головну >>> Метод Гаусса для вирішення систем лінійних рівнянь
Тут описаний алгоритм розв'язання системи лінійних рівнянь за допомогою так званого методу Гаусса. Програму ви можете скачати розділі програми . Алгоритм реалізований на мові С.
Нехай у нас є система N лінійних рівнянь
a11x1 + a12x2 + a13x3 + ... a1NxN = b1
a21x1 + a22x2 + a23x3 + ... a2NxN = b2
a31x1 + a32x2 + a33x3 + ... a3NxN = b3
...
aN1x1 + aN2x2 + aN3x3 + ... aNNxN = bN
де xi - невідомі, aij - коефіцієнти при невідомих, bi - вільні члени в рівняннях, i, j пробігають значення від 1 до N.
Мета завдання - знаючи aij і bi знайти xi.
Суть методу Гаусса полягає в тому, що за допомогою деяких операцій вихідну систему рівнянь можна звести до більш простій системі. Ця проста система має трикутний вигляд:
a11x1 + a12x2 + a13x3 + ... a1NxN = b1 a22x2 + a23x3 + ... a2NxN = b2 a33x3 + ... a3NxN = b3 ... ... aNNxN = bNОсобливість цієї системи - в рядках з номером i все коефіцієнти aij при j <i дорівнюють нулю.
Якщо ми змогли привести нашу систему рівнянь до такого трикутного вигляду, то вирішити рівняння вже просто. З останнього рівняння знаходимо xN = bN / aNN. Далі підставляємо його в передостаннє рівняння і знаходимо з нього xN-1. Підставляємо обидва знайдених рішення в наступне з кінця рівняння і знаходимо xN-2. І так далі, поки не знайдемо x1, на чому рішення закінчується. Така процедура називається зворотної прогоном.
Тепер перейдемо до питання як же домогтися того, щоб система стала трикутної.
З лінійної алгебри (див. Наприклад, Крутицкая Н.І., Тихонравов А.В., Шишкін А.А., Аналітична геометрія і лінійна алгебра з додатками) відомо що якщо до деякої рядку системи рівнянь додати будь-яку лінійну комбінацію будь-яких інших рядків цієї системи, то рішення системи не зміниться. Під лінійною комбінацією рядків розуміється сума рядків, кожна з яких у множаться на деяке число (в принципі, будь-який).
Потрібно, щоб у другому рядку вийшло рівняння, в якій відсутній член при x1. Додамо до цього рядка перший рядок, помножену на деяке число M.
(A11x1 + a12x2 + a13x3 + ... a1NxN = b1) * M +
a21x1 + a22x2 + a23x3 + ... a2NxN = b2
отримаємо
(A11 * М + a21) x1 + ... = b1 * M + b2
Для того, щоб член при x1 дорівнював нулю, потрібно, щоб M = - a21 / a11. Виконавши цю операцію, вийшло рівняння запишемо замість другого і приступимо до третього рівняння. До нього ми додамо перше рівняння, помножене на M = - a31 / a11 і теж отримаємо нуль замість члена при x1. Таку операцію потрібно виконати над усіма іншими рівняннями. В результаті отримаємо систему такого виду:
a11x1 + a12x2 + a13x3 + ... a1NxN = b1 a22x2 + a23x3 + ... a2NxN = b2 a32x2 + a33x3 + ... a3NxN = b3 ... aN2x2 + aN3x3 + ... aNNxN = bNПісля цього будемо позбуватися від членів при x2 в третьому, четвертому, N-му рівнянні. Для цього потрібно до рівняння з j-м номером додати 2-е рівняння, помножене на M = - aj2 / a22. Виконавши цю операцію над усіма іншими рівняннями, отримаємо систему де немає членів з x2 в рівняннях з номером більше 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 () - основна функція з якої послідовно викликаються всі перераховані вище функції, і проходить процес приведення системи рівнянь до трикутного вигляду і зворотна прогін.
Програма забезпечена коментарями, що полегшує її розуміння.
Сподіваюся, що вам стало все зрозуміло.
На головну >>>
