Компьютерный форум OSzone.net  

Компьютерный форум OSzone.net (http://forum.oszone.net/index.php)
-   Программирование и базы данных (http://forum.oszone.net/forumdisplay.php?f=21)
-   -   функция для вычисления матричной экспонеты (С++) (http://forum.oszone.net/showthread.php?t=118967)

Luzuk 03-10-2008 06:24 914715

функция для вычисления матричной экспонеты (С++)
 
Помогите пожалуйста написать функция для вычисления матричной экспоненты.
Ни как не могу полностью разобраться, как она вообще считается и по этому и сам алгоритм её расчёта написать не могу.

Нашёл готовую функцию, в которой вроде бы это уже реализовано. только с ней я полностью тоже разобраться не смог, слишком много непонятных входных данных.
Исходники нашёл вот в этом документе: http://www.rosavtodor.ru/doc/STO/STO_MADI.doc

Ниже привожу функцию, из этого файла, которая судя по описанию вычисляет матричную экспоненту и матрицу Аа

Код:

#define EPS 0.1e-15
#define Max(A,B) (((A)>(B))?(A):(B))
#define Min(A,B) (((A)<(B))?(A):(B))
#define Abs(A) (((A) > 0.)?(A): -(A))

expmm(double **a, double dt, int n, double **ea, double **aa, double **wm, double *w)
 {
  double am, em, emi;
  int i, j, k, ii;
    em = 0.;
    for( i = 0; i < n; i++ ){
      for( j = 0; j < n; j++ ){
        ea[i][j] = 0.;
        aa[i][j] = 0.;
        wm[i][j] = 0.;
        a[i][j] *= dt;
        am = Abs(a[i][j]);
        em = Max(am,em);
      }
      ea[i][i] = 1.;
      aa[i][i] = 1.;
      wm[i][i] = 1.;
    }
    emi = 1.;
    ii = 0;
    while( emi > EPS ){
      ii++;
      if( ii >= 40 ) break;
      emi = 0.;
      for( j = 0; j < n; j++ ){
        for( i = 0; i < n; i++ )
          w[i] = wm[i][j];
        for( i = 0; i < n; i++ ){
          wm[i][j] = 0.;
          for( k = 0; k < n; k++ )
            wm[i][j] += a[i][k] * w[k];
        }
      }
      for( i = 0; i < n; i++ )
        for( j = 0; j < n; j++ ){
          wm[i][j] /= (double)ii;
          ea[i][j] += wm[i][j];
          aa[i][j] += wm[i][j] / (double)(ii + 1);
          am = Abs(wm[i][j]);
          emi = Max(am,emi);
        }
    }
  return ii;
 }


Luzuk 03-10-2008 07:07 914724

Ещё у меня есть проверочный вариант расчёта матричной экспоненты.

Матрица A
|1 2|
|3 4|

Её экспонента
|51.9690 74.7366 |
|112.1048 164.0738|

Вот только как это получается, не знаю. Было бы неплохо если бы мне по шагам объяснили, как этот пример решается. А там думаю я уже и сам функцию смогу написать.

Luzuk 03-10-2008 09:23 914763

Ну все.. я сам разобрался.

Входной матрицей является аргумент double **a
А экспонента матрицы а, выходит через аргумент double **ea

При этом
int n - это размерность матрицы
double dt - число на которое умножаются все элементы матрицы (это вроде бы не нужный для моей задачи аргумент)
Остальное, это промежуточные параметры, которые можно спрятать внутри функции ;)

Вроде бы так...

Luzuk 03-10-2008 11:42 914893

В общем, я функцию переработал в нормальный вид и убрал лишнее.

На вход функции подается указатель на массив с матрицей и её размерность.
На выходе получаем указатель на массив с матричной экспонентой.

Ниже привожу код готовой программы с этой функцией и примером её использования, может кому ещё пригодиться ;)
Код тестировался на Visual C++ 6
Код:

#include <stdio.h>
#include <math.h>
#include <conio.h >


#define EPS 0.1e-15
#define Max(A,B) (((A)>(B))?(A):(B))
#define Min(A,B) (((A)<(B))?(A):(B))
#define Abs(A) (((A) > 0.)?(A): -(A))

double *expmm(double *a, int n)
 {
  int i, j, k, ii;
  double am, em, emi;

  double *w = new double[n*n];

  double **wm = new double*[n];
  for(i=0;i<n;i++)
          wm[i] = new double[n];

  double *ea = new double[n*n];

    em = 0.;
    for( i = 0; i < n; i++ ){
      for( j = 0; j < n; j++ ){
        ea[i*n+j] = 0.;
        wm[i][j] = 0.;
        am = Abs(a[i*n+j]);
        em = Max(am,em);
      }
      ea[i*n+i] = 1.;
      wm[i][i] = 1.;
    }
    emi = 1.;
    ii = 0;
    while( emi > EPS ){
      ii++;
      if( ii >= 40 ) break;
      emi = 0.;
      for( j = 0; j < n; j++ ){
        for( i = 0; i < n; i++ )
          w[i] = wm[i][j];
        for( i = 0; i < n; i++ ){
          wm[i][j] = 0.;
          for( k = 0; k < n; k++ )
            wm[i][j] += a[i*n+k] * w[k];
        }
      }
      for( i = 0; i < n; i++ )
        for( j = 0; j < n; j++ ){
          wm[i][j] /= (double)ii;
          ea[i*n+j] += wm[i][j];
          am = Abs(wm[i][j]);
          emi = Max(am,emi);
        }
    }
 
  delete w;
  delete []wm;

  return ea;
 }

void main()
{
        int n=2;
        int i,j;


        double *a = new double[n];
        double *ea;

        for(i=0;i<n*n;i++) a[i]=i+1;

ea = expmm(a,n);

printf("MATRIX A:\n");
for( i=0; i<n; i++){
        for( j=0; j<n; j++){
                printf("%.1f ",a[i*n+j]);
        }printf("\n");
}
printf("\n\nEXP MATRIX A:\n");
for( i=0; i<n; i++){
        for( j=0; j<n; j++){
                printf("%.3f ",ea[i*n+j]);
        }printf("\n");
}

delete []ea;

printf("\n\nPress any key to quit\n");
getch();

}



Время: 02:23.

Время: 02:23.
© OSzone.net 2001-