 |
|
функция для вычисления матричной экспонеты (С++)
Помогите пожалуйста написать функция для вычисления матричной экспоненты.
Ни как не могу полностью разобраться, как она вообще считается и по этому и сам алгоритм её расчёта написать не могу.
Нашёл готовую функцию, в которой вроде бы это уже реализовано. только с ней я полностью тоже разобраться не смог, слишком много непонятных входных данных.
Исходники нашёл вот в этом документе: 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;
}
|
Ещё у меня есть проверочный вариант расчёта матричной экспоненты.
Матрица A
|1 2|
|3 4|
Её экспонента
|51.9690 74.7366 |
|112.1048 164.0738|
Вот только как это получается, не знаю. Было бы неплохо если бы мне по шагам объяснили, как этот пример решается. А там думаю я уже и сам функцию смогу написать.
|
Ну все.. я сам разобрался.
Входной матрицей является аргумент double **a
А экспонента матрицы а, выходит через аргумент double **ea
При этом
int n - это размерность матрицы
double dt - число на которое умножаются все элементы матрицы (это вроде бы не нужный для моей задачи аргумент)
Остальное, это промежуточные параметры, которые можно спрятать внутри функции ;)
Вроде бы так...
|
В общем, я функцию переработал в нормальный вид и убрал лишнее.
На вход функции подается указатель на массив с матрицей и её размерность.
На выходе получаем указатель на массив с матричной экспонентой.
Ниже привожу код готовой программы с этой функцией и примером её использования, может кому ещё пригодиться ;)
Код тестировался на 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();
}
|
Время: 10:13.
© OSzone.net 2001-