Как инвертировать матрицу

Delphi , Синтаксис , Массивы

Автор: http://www.swissdelphicenter.ch

type
  RCOMat = array of array of Extended;

var
  DimMat: integer;

procedure InvertMatrix(var aa: RCOMat);
var
  numb, nula1, ipiv, indxr, indxc: array of Integer;
  i, j, l, kod, jmax, k, ll, icol, irow: Integer;
  amax, d, c, pomos, big, dum, pivinv: Double;
  ind: Boolean;
begin
  for j := 0 to Pred(DimMat) do
    ipiv[j] := 0;

  irow := 1;
  icol := 1;
  for i := 0 to Pred(DimMat) do
  begin
    big := 0;

    for j := 0 to Pred(DimMat) do
    begin
      if (ipiv[j] <> 1) then
      begin
        for k := 0 to Pred(DimMat) do
        begin
          if (ipiv[k] = 0) then
            if (Abs(aa[j, k]) >= big) then
            begin
              big := Abs(aa[j, k]);
              irow := j;
              icol := k;
            end
            else
              ;
        end;
      end;
    end;

    ipiv[icol] := ipiv[icol] + 1;
    if (irow <> icol) then
    begin
      for l := 0 to Pred(DimMat) do
      begin
        dum := aa[irow, l];
        aa[irow, l] := aa[icol, l];
        aa[icol, l] := dum;
      end;
      for l := 0 to Pred(DimMat) do
      begin
        dum := aa[irow + DimMat + 1, l];
        aa[irow + DimMat + 1, l] := aa[icol + DimMat + 1, l];
        aa[icol + DimMat + 1, l] := dum;
      end;
    end;
    indxr[i] := irow;
    indxc[i] := icol;
    if (aa[icol, icol] = 0) then
      ;
    pivinv := 1.0 / aa[icol, icol];
    aa[icol, icol] := 1.0;
    for l := 0 to Pred(DimMat) do
      aa[icol, l] := aa[icol, l] * pivinv;
    for l := 0 to Pred(DimMat) do
      aa[icol + DimMat + 1, l] :=
        aa[icol + DimMat + 1, l] * pivinv;
    for ll := 0 to Pred(DimMat) do
    begin
      if (ll <> icol) then
      begin
        dum := aa[ll, icol];
        aa[ll, icol] := 0.0;
        for l := 0 to Pred(DimMat) do
          aa[ll, l] := aa[ll, l] - aa[icol, l] * dum;
        for l := 0 to Pred(DimMat) do
          aa[ll + DimMat + 1, l] :=
            aa[ll + DimMat + 1, l] - aa[icol + DimMat + 1, l] * dum;
      end;
    end;
  end;

  for l := Pred(DimMat) downto 0 do
  begin
    if (indxr[l] <> indxc[l]) then
    begin
      for k := 0 to Pred(DimMat) do
      begin
        dum := aa[k, indxr[l]];
        aa[k, indxr[l]] := aa[k, indxc[l]];
        aa[k, indxc[l]] := dum;
      end;
    end;
  end;
end;

Перевод контента на русский язык:

Это реализация алгоритма Гаусса для инвертирования матрицы в Delphi. Алгоритм использует частичное пивотирование для предотвращения деления на ноль и улучшения численной стабильности.

Вот разблокировка работы кода:

1. Инициализация: Процедура InvertMatrix принимает 2D-массив aa типа RCOMat в качестве входного параметра, представляющего матрицу для инвертирования.
2. Перебор строк: Внешний цикл проходит по каждой строке матрицы (от 0 до Pred( DimMat ) - 1, где DimMat - целочисленная переменная, хранящая размер матрицы).
3. Поиск элемента пивота: Внутри цикла код ищет максимальное абсолютное значение в текущей колонке (icol) и обменивает строки, если необходимо, чтобы обеспечить, что элемент пивота находится на диагонали.
4. Обмен строк: Код обменивает строки irow и icol (если они не равны) для перемещения элемента пивота на диагональ.
5. Удаление: Код затем удаляет элементы в текущей колонке ниже элемента пивота, вычитая множители строки пивота из других строк.
6. Инвертирование элемента пивота: Код рассчитывает обратное значение элемента пивота и обновляет матрицу соответствующим образом.
7. Повторение процесса: Шаги 3-6 повторяются для каждой строки, пока не будут обработаны все строки.

Конечная инвертированная матрица хранится в входном массиве aa.

Обратите внимание, что эта реализация не проверяет ошибки, такие как сингулярные матрицы или матрицы не квадратной формы. Она также предполагает, что входная матрица является квадратной матрицей (т.е., имеет одинаковое количество строк и столбцов).

Вот альтернативное решение с использованием более современных языковых функций Delphi: generics:

type
  TMatrix<T> = array of array of T;

procedure InvertMatrix(var Matrix: TMatrix<Double>);
var
  I, J, K, L: Integer;
  Pivot: Double;
begin
  for I := 0 to Pred(Matrix.Length - 1) do
  begin
    Pivot := 0.0;
    for J := 0 to Pred(Matrix[I].Length - 1) do
    begin
      if Abs(Matrix[J, I]) > Pivot then
      begin
        Pivot := Abs(Matrix[J, I]);
        Matrix.J := Matrix.I;
        Matrix.I := J;
      end;
    end;

    // Удаление шага
    for K := 0 to Pred(Matrix[I].Length - 1) do
    begin
      if K <> Matrix.I then
      begin
        Pivot := Matrix[K, Matrix.I];
        Matrix[K, Matrix.I] := 0.0;
        for L := 0 to Pred(Matrix[I].Length - 1) do
          Matrix[K, L] := Matrix[K, L] - Pivot * Matrix[I, L];
      end;
    end;

    // Инвертирование элемента пивота
    Pivot := 1.0 / Matrix[Matrix.I, Matrix.I];
    for K := 0 to Pred(Matrix[I].Length - 1) do
      Matrix[Matrix.I, K] := Matrix[Matrix.I, K] * Pivot;
  end;
end;

Это реализация более компактная и легче для чтения, благодаря использованию generics и современных языковых функций Delphi. Однако она все еще использует схожий алгоритм, как оригинальный код.

Описано алгоритмическое решение задачи инвертирования матрицы в языке Delphi.

Комментарии и вопросы

Получайте свежие новости и обновления по Object Pascal, Delphi и Lazarus прямо в свой смартфон. Подпишитесь на наш Telegram-канал delphi_kansoftware и будьте в курсе последних тенденций в разработке под Linux, Windows, Android и iOS

:: Главная :: Массивы ::