 LU разложение. |
ВНИМАНИЕ!
Прежде чем задать вопрос, смотрите FAQ.
Рекомендуем загрузить DRKB.
  |
| Krjuger |
 1.11.2011 17:52
Сообщение
#1
|
|
Профи  Группа: Пользователи Сообщений: 652 Пол: Мужской Реальное имя: Алексей Репутация:  20   |
Собственно есть задача, полную ее суть писать не стану,но проблема заключается в следующем нужно с помощью LU разложения решить систему.Рабочий код на С++ у меня был и я его попытался перенести на Delphi.Но к несчастью моя затея обернулась неудачей.
type
TVec = array of Extended; // векторный тип
TMatr = array of TVec; // матричный тип
function LUrazl(matr: Tmatr): Tvec;
var
i, j,k, row, col: Word;
res: TVec;
_matr_: TMatr;
L,U: TMatr;
Sum: double;
begin
_matr_ := CopyMatr(matr); // Copy не работает как надо
row := Length(_matr_);
col := Length(_matr_[0]);
SetLength(L,row,col);
SetLength(U,row,col);
if row <> col-1 then
ShowMessage('Ошибка: некорректные размерности матрицы! (МГ)'); //error
for i:=0 to row-1 do
for j:= 0 to col-1 do
begin
L [i][j] := 0;
U [i][j] := 0;
if i = j then
L [i][j] := 1;
end;
//==============================================
//находим первый столбец L[][] и первую строку U[][]
for i:= 0 to row-1 do
begin
L [i][0]:= _matr_ [i][0]/ _matr_ [0][0];
U [0][i]:= _matr_ [0][i] / L [0][0];
end;
for i:= 0 to row-1 do
begin
for j:= 0 to col-1 do
begin
U [0][ i] := _matr_[0][ i];
L [i][ 0] := _matr_[i][ 0] / U[0][ 0];
sum := 0;
for k:= 0 to i do
begin
sum :=sum+ L[i][ k] * U[k][ j];
end;
U[i][ j] := _matr_[i][ j] - sum;
if (i > j) then
L[j][ i] := 0
else
begin
sum := 0;
for k:= 0 to i do
begin
sum:=sum+L[j][ k] * U[k][ i];
end;
L[j][ i] := (_matr_[j][ i] - sum) / U[i][ i];//В этой строке выдает деление на ноль так же в дебаге видно,что значение i=65000+
end;
end;
res[i]:= U[i][ i];
end;
end;
Непонятна причина такого поведения,ведь на С++ все это работало.Если надо могу добавить файл целиком. |
   |
 
|
| IUnknown |
1.11.2011 18:00
Сообщение
#2
|
 a.k.a. volvo877 Группа: Пользователи Сообщений: 1 013 Пол: Мужской Репутация: 627 |
Добавь лучше С++-ную реализацию. Посмотрим, что сделано не так, почему не работает...
|
|
|
|
| Krjuger |
1.11.2011 18:15
Сообщение
#3
|
|
Профи Группа: Пользователи Сообщений: 652 Пол: Мужской Реальное имя: Алексей Репутация: 20 |
Вот С++,так как оно мне требовалось, чтобы посчитать стопочку типовых расчетов,особого энтузиазма в оформлении нет.
#include <iostream>
using namespace std;
int main ()
{
const int n=4;
double sum = 0;
double A [ n ][ n ]= {
{ 2, -5,6,1},
{ -16, 48,-55,-5},
{ 0, 24,-25,10},
{8,21,30,4 }};
double L [ n ][ n ];
double U [ n ][ n ];
//задаем матрицу A[][] ...
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
// cout << "\na[" << i+1 << "][" << j+1 << "] = ";
// cin >> A [i][j];
L [i][j] = 0;
U [i][j] = 0;
if (i == j)
L [i][j] = 1;
}
}
//==============================================
//находим первый столбец L[][] и первую строку U[][]
for (int i = 0; i < n; i++)
{
L [i][0] = A [i][0]/ A [0][0];
U [0][i] = A [0][i] / L [0][0];
}
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
U [0][ i] = A[0][ i];
L [i][ 0] = A[i][ 0] / U[0][ 0];
double sum = 0;
for (int k = 0; k < i; k++)
{
sum += L[i][ k] * U[k][ j];
}
U[i][ j] = A[i][ j] - sum;
if (i > j)
{
L[j][ i] = 0;
}
else
{
sum = 0;
for (int k = 0; k < i; k++)
{
sum += L[j][ k] * U[k][ i];
}
L[j][ i] = (A[j][ i] - sum) / U[i][ i];
}
}
}
//====================================================
cout << "\n\n";
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
if (i == j)
L [i][j] = 1;
cout << " " << L [i][j] << " ";
}
cout << "\n\n";
}
cout << "\n\n";
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
cout << " " << U [i][j] << " ";
cout << "\n\n";
}
return 0;
}
Так как типовые были сдаты, и пара вариантов были посчитаны руками и результаты были одинаковыми,поэтому я счел,что он работает верно. Да и забыл упомянуть для текущей задачи мне требуеться,чтобы оно считалось для максимальных размерностей 101х100.А то у меня уже был случай,когда Гаусс верно считал для любой матрицы меньше 20,а для больших переставал работать,пришлось полностью пересмотреть концепцию))) Сообщение отредактировано: Krjuger - 1.11.2011 18:18 |
|
|
|
| IUnknown |
1.11.2011 18:55
Сообщение
#4
|
|
a.k.a. volvo877 Группа: Пользователи Сообщений: 1 013 Пол: Мужской Репутация: 627 |
У тебя как минимум в двух местах неправильно сделан перенос с С++... Оба - циклы по K. И там и там надо ходить до (i - 1), а не до i... Тестировал программу на той же матрице, что и Сишный вариант?
P.S. Матрицы вот такие должны получаться (это результат работы твоего чуть-чуть подпиленного дельфийского кода): L-matrix: ? Да, кстати, у тебя там балаган с возвращаемыми значениями, ты под Res память-то не выделил. Да и как собрался возвращать и L- и U-матрицы? Проще передать через параметры... Сообщение отредактировано: IUnknown - 1.11.2011 19:20 |
|
|
|
| Krjuger |
1.11.2011 19:21
Сообщение
#5
|
|
Профи Группа: Пользователи Сообщений: 652 Пол: Мужской Реальное имя: Алексей Репутация: 20 |
Нет не тестировал,я смысла не вижу,если оно даже не отрабатывает до конца,то говорить правильности решения пока рано.
Цитата У тебя как минимум в двух местах неправильно сделан перенос с С++... Оба - циклы по K. И там и там надо ходить до (i - 1), а не до i... К несчастью так нельзя,когда мы самый первый раз заходим в цикл по i у нас i=0, и i-1 в цикле по K мне даст большую беду,в Сях то такой цикл просто проскочился бы не отробатывая,а тут так не катит((( Цитата Да, кстати, у тебя там балаган с возвращаемыми значениями, ты под Res память-то не выделил. Да и как собрался возвращать и L- и U-матрицы? Проще передать через параметры... Мне не нужно возвращать эти матрицы,мне нужно вернуть вектор, который состоит из элементов диагонали матрици U это и будут решения системы,такой уж метод. Цитата Матрицы вот такие должны получаться (это результат работы твоего чуть-чуть подпиленного дельфийского кода): В Сях то он таким и получаеться. Сообщение отредактировано: Krjuger - 1.11.2011 19:26 |
|
|
|
| IUnknown |
1.11.2011 19:25
Сообщение
#6
|
|
a.k.a. volvo877 Группа: Пользователи Сообщений: 1 013 Пол: Мужской Репутация: 627 |
Цитата К несчастью так нельзя,когда мы самый первый раз заходим в цикл по i у нас i=0, и i-1 в цикле по K мне даст большую беду Бла-бла-бла... Кроме Word-а в Дельфи есть еще и Integer, который способен принимать отрицательные значения. Чуешь, откуда ошибка?В общем, вот, прогнал на FPC твой дельфийский код, ничего не вылетает, печатает то, что приведено выше: type
TVec = array of Extended;
TMatr = array of TVec;
procedure LUrazl(_matr_: Tmatr; var L, U : TMatr);
var
i, j, row, col: Word;
k : integer;
res: TVec;
Sum: double;
begin
// _matr_ := CopyMatr(matr);
row := Length(_matr_);
col := Length(_matr_[0]);
SetLength(L,row,col);
SetLength(U,row,col);
if row <> col-1 then
writeln('error');
for i:=0 to row-1 do
for j:= 0 to col-1 do
begin
L [i][j] := 0;
U [i][j] := 0;
if i = j then L [i][j] := 1;
end;
//==============================================
for i:= 0 to row-1 do
begin
L [i][0]:= _matr_ [i][0]/ _matr_ [0][0];
U [0][i]:= _matr_ [0][i] / L [0][0];
end;
for i:= 0 to row-1 do
begin
for j:= 0 to col-1 do
begin
U [0][ i] := _matr_[0][ i];
L [i][ 0] := _matr_[i][ 0] / U[0][ 0];
sum := 0;
for k:= 0 to i-1 do
begin
sum :=sum+ L[i][ k] * U[k][ j];
end;
U[i][ j] := _matr_[i][ j] - sum;
if (i > j) then
L[j][ i] := 0
else
begin
sum := 0;
for k:= 0 to i-1 do
begin
sum:=sum+L[j][ k] * U[k][ i];
end;
L[j][ i] := (_matr_[j][ i] - sum) / U[i][ i];
end;
end;
end;
end;
var
A : tmatr;
L, U : tMatr;
i, j : integer;
begin
setlength(a, 4, 4);
a[0, 0] := 2; a[0, 1] := -5; a[0, 2] := 6; a[0, 3] := 1;
a[1, 0] := -16; a[1, 1] := 48; a[1, 2] := -55; a[1, 3] := -5;
a[2, 0] := 0; a[2, 1] := 24; a[2, 2] := -25; a[2, 3] := 10;
a[3, 0] := 8; a[3, 1] := 21; a[3, 2] := 30; a[3, 3] := 4;
LURazl(A, L, U);
writeln('L-matrix:');
for i := 0 to Length(L) - 1 do
begin
for j := 0 to Length(L[0]) - 1 do
write(L[i, j]:10 :4);
writeln;
end;
writeln('U-matrix:');
for i := 0 to Length(U) - 1 do
begin
for j := 0 to Length(U[0]) - 1 do
write(U[i, j]:10 :4);
writeln;
end;
end.
Дельфи запускать лень  Сообщение отредактировано: IUnknown - 1.11.2011 19:26 |
|
|
|
| -Татьяна- |
14.06.2012 0:59
Сообщение
#7
|
|
Гость |
Цитата(Krjuger @ 1.11.2011 17:52)  Собственно есть задача, полную ее суть писать не стану,но проблема заключается в следующем нужно с помощью LU разложения решить систему.Рабочий код на С++ у меня был и я его попытался перенести на Delphi.Но к несчастью моя затея обернулась неудачей.
type
TVec = array of Extended; // векторный тип
TMatr = array of TVec; // матричный тип
function LUrazl(matr: Tmatr): Tvec;
var
i, j,k, row, col: Word;
res: TVec;
_matr_: TMatr;
L,U: TMatr;
Sum: double;
begin
_matr_ := CopyMatr(matr); // Copy не работает как надо
row := Length(_matr_);
col := Length(_matr_[0]);
SetLength(L,row,col);
SetLength(U,row,col);
if row <> col-1 then
ShowMessage('Ошибка: некорректные размерности матрицы! (МГ)'); //error
for i:=0 to row-1 do
for j:= 0 to col-1 do
begin
L [i][j] := 0;
U [i][j] := 0;
if i = j then
L [i][j] := 1;
end;
//==============================================
//находим первый столбец L[][] и первую строку U[][]
for i:= 0 to row-1 do
begin
L [i][0]:= _matr_ [i][0]/ _matr_ [0][0];
U [0][i]:= _matr_ [0][i] / L [0][0];
end;
for i:= 0 to row-1 do
begin
for j:= 0 to col-1 do
begin
U [0][ i] := _matr_[0][ i];
L [i][ 0] := _matr_[i][ 0] / U[0][ 0];
sum := 0;
for k:= 0 to i do
begin
sum :=sum+ L[i][ k] * U[k][ j];
end;
U[i][ j] := _matr_[i][ j] - sum;
if (i > j) then
L[j][ i] := 0
else
begin
sum := 0;
for k:= 0 to i do
begin
sum:=sum+L[j][ k] * U[k][ i];
end;
L[j][ i] := (_matr_[j][ i] - sum) / U[i][ i];//В этой строке выдает деление на ноль так же в дебаге видно,что значение i=65000+
end;
end;
res[i]:= U[i][ i];
end;
end;
Непонятна причина такого поведения,ведь на С++ все это работало.Если надо могу добавить файл целиком. А ВЫ не могли бы сбросить вашу задачу пожалуйста, она мне очень нужна...(lovetanya.15@mail.ru) Спасибо |
|
|
|
|
1 чел. читают эту тему (гостей: 1, скрытых пользователей: 0)
Пользователей: 0
|
Текстовая версия | 26.07.2025 11:24 |