[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();

}