Собственно есть задача, полную ее суть писать не стану,но проблема заключается в следующем нужно с помощью 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;
Добавь лучше С++-ную реализацию. Посмотрим, что сделано не так, почему не работает...
Вот С++,так как оно мне требовалось, чтобы посчитать стопочку типовых расчетов,особого энтузиазма в оформлении нет.
#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;
}
У тебя как минимум в двух местах неправильно сделан перенос с С++... Оба - циклы по K. И там и там надо ходить до (i - 1), а не до i... Тестировал программу на той же матрице, что и Сишный вариант?
P.S.
Матрицы вот такие должны получаться (это результат работы твоего чуть-чуть подпиленного дельфийского кода):L-matrix:
1.0000 0.0000 0.0000 0.0000
-8.0000 1.0000 0.0000 0.0000
0.0000 3.0000 1.0000 0.0000
4.0000 5.1250 -10.4688 1.0000
U-matrix:
2.0000 -5.0000 6.0000 1.0000
0.0000 8.0000 -7.0000 3.0000
0.0000 0.0000 -4.0000 1.0000
0.0000 0.0000 0.0000 -4.9063
?
Да, кстати, у тебя там балаган с возвращаемыми значениями, ты под Res память-то не выделил. Да и как собрался возвращать и L- и U-матрицы? Проще передать через параметры...
Нет не тестировал,я смысла не вижу,если оно даже не отрабатывает до конца,то говорить правильности решения пока рано.
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.
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;