функция для вычисления матричной экспонеты (С++)

[Luzuk]

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

Нашёл готовую функцию, в которой вроде бы это уже реализовано. только с ней я полностью тоже разобраться не смог, слишком много непонятных входных данных.
Исходники нашёл вот в этом документе: 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]

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

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

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

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

[Luzuk]

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

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

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

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

[Luzuk]

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

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

Ниже привожу код готовой программы с этой функцией и примером её использования, может кому ещё пригодиться ;)
Код тестировался на 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();

}