VPF::Нахождение обратной матрицы - Форум программистов

Mpak

Дата 19.12.2006, 15:33 (ссылка)

(нет голосов)

Загрузка ...

Опытный

Профиль
Группа: Участник
Сообщений: 445
Регистрация: 12.11.2006

Репутация: нет
Всего: 7

Не могу догадаться КАК....
Методом обратной матрицы

Методом Крамера сделал...
Подайте идею smile

Это сообщение отредактировал(а) Mpak - 19.12.2006, 15:33

--------------------

В любой откомпилированной программе есть, по крайней мере, одна ошибка...
P.S. А у меня их минимум две...

SoWa

Дата 19.12.2006, 16:01 (ссылка)

(нет голосов)

Загрузка ...

Харекришна

Профиль
Группа: Комодератор
Сообщений: 2422
Регистрация: 18.10.2004

Репутация: 6
Всего: 74

Два метода:
1) Методом Гаусса. Строишь расширенную матрицу (A|E), приводишь её слева к единичному виду. Получаешь (E|B) Матрица В- обратная для A.
2) Не знаю названия.
Обратная матрица- это матрица вида ((Av)т)/(|A|)
Где Av- матрица алгебраических дополнений, буква "т"- транспонирование, В знаменателе определитель исходной матрицы.
Думаю, не надо рассказывать, что такое алгебраическое дополнение?

--------------------

Всем добра smile

Mpak

Дата 19.12.2006, 16:17 (ссылка)

(нет голосов)

Загрузка ...

Опытный

Профиль
Группа: Участник
Сообщений: 445
Регистрация: 12.11.2006

Репутация: нет
Всего: 7

SoWa, я это все знаю....не за даром на второй курс был переведен smile

Но все равно спасибо...

Я не могу перевести это дело в программный код smile

--------------------

В любой откомпилированной программе есть, по крайней мере, одна ошибка...
P.S. А у меня их минимум две...

SoWa

Дата 19.12.2006, 16:33 (ссылка)

(нет голосов)

Загрузка ...

Харекришна

Профиль
Группа: Комодератор
Сообщений: 2422
Регистрация: 18.10.2004

Репутация: 6
Всего: 74

Второй метод очень легко. Достаточно сделать код нахождения определителя. Не сложно. Руккурентно.
PS Я только на первом курсе smile

PPS Ты на twilight??? Зареган?

--------------------

Всем добра smile

Mpak

Дата 19.12.2006, 16:46 (ссылка)

(нет голосов)

Загрузка ...

Опытный

Профиль
Группа: Участник
Сообщений: 445
Регистрация: 12.11.2006

Репутация: нет
Всего: 7

>PPS Ты на twilight??? Зареган?

Это цитата из песни

Нет, а что? Есть такое?

--------------------

В любой откомпилированной программе есть, по крайней мере, одна ошибка...
P.S. А у меня их минимум две...

esperant0

Дата 19.12.2006, 17:36 (ссылка)

(нет голосов)

Загрузка ...

Опытный

Профиль
Группа: Участник
Сообщений: 714
Регистрация: 20.5.2005

Репутация: 4
Всего: 14

Рекурентный поиск определителя имеет экспоненциальную сложность а потому не эффективен.

--------------------

Student->Teacher Assistant ->Research assistant->Microsoft Software Development Engineer

Пользователь получил наказание за то, что проигнорировал замечание которое было написано модератором а затем стерто и которое он - пользователь не мог видеть.

SoWa

Дата 19.12.2006, 18:01 (ссылка)

(нет голосов)

Загрузка ...

Харекришна

Профиль
Группа: Комодератор
Сообщений: 2422
Регистрация: 18.10.2004

Репутация: 6
Всего: 74

Есть быстрее?

--------------------

Всем добра smile

esperant0

Дата 19.12.2006, 18:03 (ссылка)

(нет голосов)

Загрузка ...

Опытный

Профиль
Группа: Участник
Сообщений: 714
Регистрация: 20.5.2005

Репутация: 4
Всего: 14

Цитата(SoWa @ 19.12.2006, 18:01)

Есть быстрее?

Да, метод приведения к единичной матрице, имеет кубическую сложность

--------------------

maxim1000

Дата 19.12.2006, 18:07 (ссылка)

(нет голосов)

Загрузка ...

Эксперт

Профиль
Группа: Участник
Сообщений: 3334
Регистрация: 11.1.2003
Где: Киев

Репутация: 33
Всего: 110

Цитата(SoWa @ 19.12.2006, 17:01

)

Есть быстрее?

Гаусса уже упоминали
он меет кубическую сложность

--------------------

qqq

SoWa

Дата 19.12.2006, 18:53 (ссылка)

(нет голосов)

Загрузка ...

Харекришна

Профиль
Группа: Комодератор
Сообщений: 2422
Регистрация: 18.10.2004

Репутация: 6
Всего: 74

Так. Запутались вы.
Я спросил про алгоритм вычисления определителя, более быстрый, чем рекурентный.

А то, что метод Гаусса быстрее метода алгебраических дополнений я не спорил smile

Теперь человеку в раздел "Центр помощи"

Это сообщение отредактировал(а) SoWa - 19.12.2006, 20:09

--------------------

Всем добра smile

maxim1000

Дата 20.12.2006, 03:10 (ссылка)

(нет голосов)

Загрузка ...

Эксперт

Профиль
Группа: Участник
Сообщений: 3334
Регистрация: 11.1.2003
Где: Киев

Репутация: 33
Всего: 110

Цитата(SoWa @ 19.12.2006, 17:53

)

Я спросил про алгоритм вычисления определителя, более быстрый, чем рекурентный.

А то, что метод Гаусса быстрее метода алгебраических дополнений я не спорил

ну так метод Гаусса вполне может быть использован для вычисления определителя...

Добавлено @ 03:11
кстати, тут: http://forum.vingrad.ru/topic-37776/view-all.html
по-моему, есть его реализация, только для обращения матрицы
если его подправить, можно попутно и определитель считать...

--------------------

qqq

SoWa

Дата 20.12.2006, 04:52 (ссылка)

(нет голосов)

Загрузка ...

Харекришна

Профиль
Группа: Комодератор
Сообщений: 2422
Регистрация: 18.10.2004

Репутация: 6
Всего: 74

Хм. Ну все равно Метод гаусса почти реккурентный:
Приводим первый столбец к базисному виду. Определитель равен прозведению единицы на минор n-1 порядка. Ищем минор далее...
Хотя можно привести к диагональной единичной иатрице и все без проблемм рассчитается.
Разобрались.

--------------------

Всем добра smile

V.A.KeRneL

Дата 21.12.2006, 02:05 (ссылка)

(нет голосов)

Загрузка ...

Vadim A. Kazantsev

Профиль
Группа: Участник
Сообщений: 291
Регистрация: 3.12.2006
Где: Moscow, Russia

Репутация: 1
Всего: 14

Вот тут довольно эффективная итерационная реализация метода нахождения определителя матрицы из стандартной библиотеки Ruby (файл "matrix.rb"):

Код


# ...

class Matrix
  @RCS_ID='-$Id: matrix.rb,v 1.11 1999/10/06 11:01:53 keiju Exp keiju $-'

# ...

  #
  # Returns the determinant of the matrix.  If the matrix is not square, the
  # result is 0.
  #   Matrix[[7,6], [3,9]].determinant
  #     => 63
  #
  def determinant
    return 0 unless square?
    
    size = row_size - 1
    a = to_a
    
    det = 1
    k = 0
    begin 
      if (akk = a[k][k]) == 0
        i = k
        begin
          return 0 if (i += 1) > size
        end while a[i][k] == 0
        a[i], a[k] = a[k], a[i]
        akk = a[k][k]
        det *= -1
      end
      (k + 1).upto(size) do
        |i|
        q = a[i][k] / akk
        (k + 1).upto(size) do
          |j|
          a[i][j] -= a[k][j] * q
        end
      end
      det *= akk
    end while (k += 1) <= size
    det
  end
  alias det determinant

# ...

end  # class Matrix

# ...

Тут обязательно надо заметить, что в комментарии приведёно неправильное значение (63) определителя матрицы Matrix[[7,6], [3,9]]. Правильное: Matrix[[7,6], [3,9]].determinant = 7*9 - 6*3 = 63 - 18 = 45.

Такой результат получится, если в программе, использующий данный класс (Matrix), подключить (require) только файл "matrix.rb", без "mathn.rb". Подключение последнего решает проблемму.

Код


(E:\__tmp__\c)
[[email protected]]$ irb
irb(main):001:0> require "matrix"
=> true
irb(main):002:0> Matrix[[7, 6], [3, 9]].determinant
=> 63
irb(main):003:0> require "mathn"
=> true
irb(main):004:0> Matrix[[7, 6], [3, 9]].determinant
=> 45
irb(main):005:0>

Вся проблема тут, насколько я понимаю, заключается в использовании в строчке «q = a[i][k] / akk» операции [целочисленного] деления для матрицы целых чисел.

Определитель матрицы вещественных (действительных) чисел всегда вычисляетя корректно, даже без подключения файлика "mathn.rb":

Код


(E:\__tmp__\c)
[[email protected]]$ irb
irb(main):001:0> require "matrix"
=> true
irb(main):002:0> Matrix[[7.0, 6.0], [3.0, 9.0]].determinant
=> 45.0
irb(main):003:0> require "mathn"
=> true
irb(main):004:0> Matrix[[7.0, 6.0], [3.0, 9.0]].determinant
=> 45.0
irb(main):005:0>

Вот переписал код на C для более широких масс:

Код


#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>

void inline SWAP_ROWS(double **, int, int, int);
double determinant(const double **, int, int);

void 
inline 
SWAP_ROWS(double ** arr2d, int row1, int row2, int cols) 
{ 
    double t;
    int j;
    
    for (j = 0; j < cols; ++j) { 
        t = arr2d[row1][j];
        arr2d[row1][j] = arr2d[row2][j];
        arr2d[row2][j] = t;
    }
}

double 
determinant(const double ** matrix, int rows, int cols) 
{ 
    int size;
    double ** a;
    double det;
    int i, j, k;
    double akk;
    double q;
    
    if (rows != cols) 
        return (0.0);
    
    size = rows - 1;
    a = (double **) matrix;
    
    det = 1.0;
    k = 0;
    do { 
        if ((akk = a[k][k]) == 0.0) { 
            i = k;
            do { 
                if (++i > size) 
                    return (0.0);
            } while (a[i][k] == 0.0);
            SWAP_ROWS(a, i, k, cols);
            akk = a[k][k];
            det *= -1.0;
        }
        for (i = k + 1; i <= size; ++i) { 
            q = a[i][k] / akk;
            for (j = k + 1; j <= size; ++j) { 
                a[i][j] -= a[k][j] * q;
            }
        }
        det *= akk;
    } while (++k <= size);
    
    return (det);
}

int 
main(void) 
{ 
    int rows = 3;
    int cols = rows;
    double ** matrix;
    int i, j;
    
    matrix = (double **) malloc(rows * sizeof(double *));
    if (matrix == NULL) 
        exit(EXIT_FAILURE);
    for (i = 0; i < rows; ++i) { 
        matrix[i] = (double *) calloc(cols, sizeof(double));
        if (matrix[i] == NULL) 
            exit(EXIT_FAILURE);
    }
    matrix[0][0] = 1; matrix[0][1] = 2; matrix[0][2] = 3;
    matrix[1][0] = 2; matrix[1][1] = 3; matrix[1][2] = 1;
    matrix[2][0] = 3; matrix[2][1] = 1; matrix[2][2] = 2;
    /* matrix = {{1, 2, 3},  */
    /*           {2, 3, 1},  */
    /*           {3, 1, 2}}; */
    /* determinant(matrix)  //=> -18.0 */
    
    printf("determinant(matrix) = %.3lf\n", 
            determinant((const double **) matrix, rows, cols));  //=> -18.000
    
    return (0);
}

З.Ы. Сорри за не самое читабельное заполнение матрицы, у меня GCC под виндой глючит...

З.З.Ы. Я не очень силён в математике (пока), но мне кажется использован именно метод Гаусса, просто оптимизирован. Если я не прав, то обязательно поправьте меня, пожалуйста!

Это сообщение отредактировал(а) V_A_KeRneL - 21.12.2006, 02:09

--------------------

«C'est un pense-creux d'ici. C'est le meilleur et le plus irascible homme du monde...» © Ф.М. Достоевский, «Бесы»
---/)/)---(\.../)---(\(\
--(':'=)---(=';'=)---(=':')
(")(")..)-(").--.(")-(..(")(")

V.A.KeRneL

Дата 21.12.2006, 02:27 (ссылка)

(нет голосов)

Загрузка ...

Vadim A. Kazantsev

Профиль
Группа: Участник
Сообщений: 291
Регистрация: 3.12.2006
Где: Moscow, Russia

Репутация: 1
Всего: 14

А вообще, коллекцию реализаций by Johna Smith таких стандартных (классических) алгоритмов как поиск обратной матрицы любезно выложил podval вот здесь: http://forum.vingrad.ru/topic-37776/view-all.html
(В разделе Матрицы, есессно.)

Код


//////////////////////////////////////////////////////////////////////////////
//
//  Matrix operations (inversion)
//  (c) Johna Smith, 1996
//
//  Given: A (N x N), det A != 0                 -1
//  This algorithm inverts matrix A and returns A  .
//      -1
//   A*A  = E
//  We simply convert matrix A into matrix E, and do the same operations
//  over matrix B (initially B=E). Finally A=E, B=A^(-1).
//
//////////////////////////////////////////////////////////////////////////////
#define N               3
#include <stdio.h>
float a[N][N]={{4,8,0},{8,8,8},{2,0,1}};
void Invert(float *matrix)
{
  float e[N][N];
  // initializing matrix e
  for (int i=0;i<N;i++)
    for (int j=0;j<N;j++)
      e[i][j]=(i==j?1:0);
  // converting matrix to e
  for(i=0;i<N;i++)
  {
    // normalizing row (making first element =1)
    float tmp=*(matrix+i*N+i);
    for(int j=N-1;j>=0;j--)
    {
      e[i][j]/=tmp;
      *(matrix+i*N+j)/=tmp;
    }
    // excluding i-th element from each row except i-th one
    for(j=0;j<N;j++)
      if (j!=i)
      {
        tmp=*(matrix+j*N+i);
        for(int k=N-1;k>=0;k--)
        {
          e[j][k]-=e[i][k]*tmp;
          *(matrix+j*N+k)-=*(matrix+i*N+k)*tmp;
        }
      }
  }
  // now e contains inverted matrix so we need only to copy e to matrix
  for(i=0;i<N;i++)
    for(int j=0;j<N;j++)
      *(matrix+i*N+j)=e[i][j];
}
void show_matrix(float *matrix) // this functions displays matrix
{
 for(int i=0;i<N;i++)
 {
   for(int j=0;j<N;j++)
     printf("%f ",*(matrix+i*N+j));
   printf("\n");
 }
 printf("\n");
}
void main(void)
{
  // display matrix A
  show_matrix(a[0]);
  // Invert it
  Invert(a[0]);
  // display the inverted matrix
  show_matrix(a[0]);
  // Invert it again
  Invert(a[0]);
  // display the inversion of the inverted matrix
  show_matrix(a[0]);
}

Ещё обязательно стоит посмотреть «Прежде чем изобретать велосипед (алгоритмы решения многих задач)» от того же podval'а.

З.Ы. smile

После этого поста я стал «Шустрым», и у меня появилась первая звёздочка! smile

Это сообщение отредактировал(а) V_A_KeRneL - 29.12.2006, 19:31

--------------------

popovda

Дата 28.12.2006, 16:29 (ссылка)

(нет голосов)

Загрузка ...

Опытный

Профиль
Группа: Участник
Сообщений: 290
Регистрация: 9.6.2006
Где: Москва

Репутация: 1
Всего: 6

Ребята. Вообще лучше быть осторожнее с методом Гаусса, т.к. он ведет к накоплению больших ошибок округления. А уж если говорить о поиске обратной матрицы (и определителя), то надо использовать или LU-разложение,метод вращений, метод ортогональных преобразований, или итерационные методы, которые не ведут к накоплению ошибки округления. Кстати, по LU-разложению на форуме мною была выложена методичка - можете найти. Или пользуйтесь Фортраном smile

. Как Вам такое:

antiA = .i. A;
DetA = det(A);

evals = eig(A); - собственные числа.
evals = eig(B,evects); - собств. числа и собств. вектора. (Можно правда и из др. языков IMSL библиотеки использовать)

Там алгоритмы отработаны.

--------------------

С уважением, Попов Д.А.

0 Пользователей читают эту тему (0 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема \| Алгоритмы \| Следующая тема »