Консультация № 172614
25.09.2009, 14:16
0.00 руб.
0 10 1
Здравствуйте, уважаемые Эксперты!

Срочно нужна Ваша помощь!

Для решения одной из моих задач мне нужно найти решение системы линейных уравнений методом Гаусса.

Поискав в интернете нашел _один_ распространенный алгоритм (приведу его ниже).
НО! Он работает только для матриц с ненулевыми значениями.

Я точно знаю, что методом Гаусса можно решить и с нулевыми, но такого примера не нашел.

Если у кого есть прошу поделитесь, пожалуйста. (Матрица n*n +значения). Т около 1000. приведенный ниже алгоритм для не нулевых значений нормально отрабатывается за 3 секунды.

P.S. Это не для курсовой и т.д.

Код:

unit Unit3;

interface

uses
Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
Dialogs, Grids, StdCtrls;

const MaxDimension = 10;

type

Vector = array[1..MaxDimension] of Double;
Matrix = array[1..MaxDimension] of Vector;

TForm3 = class(TForm)
Label1: TLabel;
Edit1: TEdit;
StringGrid1: TStringGrid;
StringGrid2: TStringGrid;
Button1: TButton;
Label2: TLabel;
Label3: TLabel;
ListBox1: TListBox;
procedure Edit1Change(Sender: TObject);
procedure Button1Click(Sender: TObject);
procedure FormCreate(Sender: TObject);
private
{ Private declarations }
public
{ Public declarations }
end;

var
Form3: TForm3;

implementation

{$R *.dfm}

procedure TForm3.Button1Click(Sender: TObject);
var a: Matrix;
b,x: Vector;
h: Double;
i,j,k,n:integer;
begin
//Ввод данных
//Размерность системы
n := StrToIntDef(Text, StringGrid1.ColCount);
//Коэффициенты
for j := 0 to n - 1 do
for i := 0 to n - 1 do
a[i + 1, j + 1] := StrToFloatDef(StringGrid1.Cells[j, i], 0);
//Правая часть уравнения
for I := 0 to n - 1 do b[i + 1] := StrToFloatDef(StringGrid2.Cells[0, i], 0);
//Прямой ход - исключение переменных
for i:=1 to n-1 do
for j:=i+1 to n do
begin
a[j,i]:=-a[j,i]/a[i,i];
for k:=i+1 to n do
a[j,k]:=a[j,k]+a[j,i]*a[i,k];
b[j]:=b[j]+a[j,i]*b[i]
end;
x[n]:=b[n]/a[n,n]; //Деление на ноль
//Обратный ход - нахождение корней
for i:=n-1 downto 1 do
begin
h:=b[i];
for j:=i+1 to n do h:=h-x[j]*a[i,j];
x[i]:=h/a[i,i]
end;
//Вывод результата
for i:=1 to n do ListBox1.Items.Append('x(' + IntToStr(i) + ')=' + FloatToStr(x[i]));
end;

procedure TForm3.Edit1Change(Sender: TObject);
begin
with StringGrid1, Edit1 do
begin
ColCount := StrToIntDef(Text, 3);
RowCount := StrToIntDef(Text, 3);
end;
with StringGrid2, Edit1 do
RowCount := StrToIntDef(Text, 3)
end;

procedure TForm3.FormCreate(Sender: TObject);
var i, j: integer;
begin
//Заполнить коэф уравнения
Randomize;
for I := 0 to StrToIntDef(Text, StringGrid1.ColCount) - 1 do
for J := 0 to StrToIntDef(Text, StringGrid1.RowCount) - 1 do
StringGrid1.Cells[I, J] := IntToStr(Random(100));
for I := 0 to StrToIntDef(Text, StringGrid2.RowCount) - 1 do
StringGrid2.Cells[0, I] := IntToStr(Random(100))
end;

end.

Обсуждение

Неизвестный
25.09.2009, 14:19
общий
Т около 1000
забыл переключить раскладку естественно - N около 1000.
P.S. Просьба помещайте решение в тэгах код - в приложении переносы строк срезаются, "дешифровать" потом практически невозможно.
давно
Профессионал
304622
583
28.09.2009, 21:34
общий
это ответ
Здравствуйте, LanK.


1) Непонятно, почему у тебя написано "Деление на ноль" в строчке с делением на a[n,n]. Оно, как правило, должно происходить раньше, в цикле при делении на a[i,i]

2) Напутано с алгоритмом - правильного ответа быть не должно. После
a[j,i]:=-a[j,i]/a[i,i];
надо делать такое же деление для b
b[j]:=-b[i]/a[i,i];

А строчка
b[j]:=b[j]+a[j,i]*b[i]
должна, во первых, находиться в цикле (у тебя нет begin end), во вторых, просчитывать k-й (а не j-й) элемент, т.е.:
b[k]:=b[k]+a[j,i]*b[k]

3) Обратный ход вполне можно делать от n:
for i:=n downto 1 do
а строчку
x[n]:=b[n]/a[n,n];
преспокойно стирать. Избыточна.

4) Теперь о нулевых диагональных элементах. Это нужен алгоритм с выбором столбцов. Собственно никакого особого алгоритма тут нет. Просто у системы можно без нарушения эквивалентности менять местами столбцы. Главное при это помнить какие столбцы переставлены.

Т.е. перед главным делением
a[j,i]:=-a[j,i]/a[i,i];
надо проверить a[i,i] на равнество нулю и, если да, то найти в строке какой-нибудь (например, первый) ненулевой. (Если все нули значит система недоопределени и не имеет единственного решения.) Это можно сделать сразу
m:=i;
while (a[i,m]<>0) and (m<=n) do inc(m);
if m>n then ... // все нули
Потом цикл по этим столбцам (от начала, а не от i)
for k:=1 to n do
begin
/// меняем местами a[i,k] и a[m,k]
end;

Останется только запоминть, какие элементы поменял. Это можно держать в массиве вроде
ind:array[1..n] of integer;
Перед прямым ходом надо его заполнить
for i:=1 to n do ind[i]:=i;
А результата выводить не как x[i], а как x[ind[i]].
(В принципе, может быть в Делфе уже есть какая-нибудь удобная штука для этого. Не помню.)

Удачи!
5
спасибо
Неизвестный
29.09.2009, 10:43
общий
Сергей Бендер:
Сергей, очень неплохо Вы всё расписали.
Приведите, пожалуйста, рабочий код (в тегах code).
А то мыслей в Вашем ответе много, при "сборке" можно и ошибиться...

P.S. В Delphi "удобная штука" нет, Вы верно написали надо циклом выводить.
давно
Профессионал
304622
583
29.09.2009, 21:50
общий
Рабочий код -- это чуть позже. Сейчас некогда.

Но пока стоит IMHO самому попробовать. Я ж практически всё расписал прямым текстом ;).
Неизвестный
30.09.2009, 08:55
общий
Сергей Бендер:
Цитата: Сергей Бендер
Я ж практически всё расписал прямым текстом ;).
Так то оно так, но я СЛАУ изучал лет 20 назад, так что программировать то я могу, а вот что там делать когда строки меняются и т.д. уже тяжко "въезжать".

Цитата: Сергей Бендер
Рабочий код -- это чуть позже.
Ок. Буду ждать. Поэтому поставлю за ответ пять. (в надежде всё таки на рабочий код).

За ранее спасибо!
давно
Профессионал
304622
583
02.10.2009, 18:55
общий
LanK:
Завтра :)
давно
Профессионал
304622
583
03.10.2009, 20:45
общий
LanK:
Ну вот.

Во первых, признаю себя ослом. То, что я предлагал, по исправлению алгоритма -- ерунда. Я недостаточно вгляделся в текст твоей программы -- там несколько другой алгоритм решения, чем тот который я привык видеть.

Во вторых, лучше, конечно, менять местами строки, а не столбцы -- там не надо запоминать перестановки в x. Я думал об этом с начала, но перестраховался -- боялся, что можно упустить некоторые случаи. Похоже что зря. Будем менять местами строки.

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

В четвёртых, деление на ноль всё равно может произойти, но это для линейно зависимой матрицы системы. Тогда решения либо нет, либо много. Это уже другая задача.

Код:

//Прямой ход - исключение переменных
for i:=1 to n-1 do
begin
j:=i;
while (a[j,i]=0) and (j<=n) do inc(j);
if j<>i then
begin
for k:=1 to n do
begin
h:=a[i,k];
a[i,k]:=a[j,k];
a[j,k]:=h;
end;
h:=b[i];
b[i]:=b[j];
b[j]:=h;
end;
for j:=i+1 to n do
begin
a[j,i]:=-a[j,i]/a[i,i];
for k:=i+1 to n do
a[j,k]:=a[j,k]+a[j,i]*a[i,k];
b[j]:=b[j]+a[j,i]*b[i]
end;
end;

//Обратный ход - нахождение корней
for i:=n downto 1 do
begin
h:=b[i];
for j:=i+1 to n do
h:=h-x[j]*a[i,j];
x[i]:=h/a[i,i]
end;
Неизвестный
05.10.2009, 08:52
общий
Сергей Бендер:
Спасибо!
Правда, я уже переписал алгоритм полностью. Без деления. Хотя он получился не особо красивым, но матрицу 5000*5000 решает далее вылетает (похоже при заполнении результатов в StringGrid).
давно
Профессионал
304622
583
05.10.2009, 13:23
общий
LanK:

Правда, я уже переписал алгоритм полностью. Без деления.


А в чём там суть? Он ведь несовместную систему всё равно не может решить?
Неизвестный
06.10.2009, 09:13
общий
Сергей Бендер:
Несовместную да, конечно не может. Просто в системе получается много уравнений - как я писал около 1000. И если использовать деление получаются очень большие погрешности. Поэтому всё просто - домнажаем и вычитаем строки. Получаем "единичную" матрицу, затем прогоняем обратно по Жордану. Получаем решение. Если по Гауссу не проходит на диагональ выпадает ноль - то решений нет.
Форма ответа