|
Текст подпрограммы и версий ame1r_c.zip , ame1d_c.zip , ame2r_c.zip , ame2d_c.zip |
Тексты тестовых примеров tame1r_c.zip , tame1d_c.zip , tame2r_c.zip , tame2d_c.zip |
Вычисление матричной экспоненты exp (AT) при большом по модулю значении скалярной переменной T.
Вычисляется функция от матрицы exp (AT) - матричная экспонента от произведения мaтpицы А и скаляpнoй переменной T при большом по модулю значении скалярной величины T по формуле:
exp(AT) = ( exp(AT/N) )N ,
где N - некоторое целое положительное число, задаваемое при обращении к подпрограмме. Матричная экспонента exp (АТ/N) = exp (AH) аппроксимируется усеченным степенным рядом
14
F = ∑ (AH)k / k!
k=0
Значение N предполагается таким, что погрешность аппроксимации exp (АН) - F является достаточно малой.
int ame1r_c (real *a, real *t, integer *m, integer *n,
real *r1, real *r, real *e, real *r2, integer *ierr)
Параметры
| a - | вещественный двумерный массив размера m на m, содержащий матрицу, входящую в матричную экспоненту; |
| t - | скалярная величина, входящая в матричную экспоненту (тип: вещественный); |
| m - | порядок матрицы A (тип: целый); |
| n - | некоторое целое положительное число, участвующее в алгоритме вычисления матричной экспоненты. Значение n рекомендуется задавать таким, чтобы погрешность аппроксимации матричной экспоненты exp (AH) усеченным степенным рядом F была достаточно малой; |
| e - | вещественный двумерный массив размера m на m, содержащий матрицу - результат exp (AT); |
| r - | вещественный одномерный рабочий массив длины m; |
| r1, r2 - | вещественные двумерные рабочие массивы размера m на m; |
| ierr - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если значение параметра n меньше 1. В этом случае вычисление матричной экспоненты прекращается. |
Версии
| ame2r_c - |
вычисление матричной экспоненты exp (AT) при большом по модулю значении скалярной переменной T по формуле exp(AT) = ( exp(AT/N )N в случае, когда N не задается при обращении к подпрограмме, а выбирается в ней автоматически из условия, чтобы || AT/N || < 1, а именно: exp(AT) = ( exp(AT / [ || A || * | T | + 1 ] ) ) [ || A || * | T |+1 ] ,где [x] - целая часть величины x, а || A || означает максимальную сумму модулей элементов матрицы A по строкам. Первый оператор подпрограммы имеет вид: int ame2r_c(real *a, real *t, integer *m, real *e, real *r__, real *r1, real *r2). Число параметров подпрограммы ame2r_c на 2 меньше числа параметров подпрограммы ame1r_c, при этом параметры подпрограммы ame2r_c имеют тот же смысл, что и одноименные параметры ame1r_c. |
|
ame1d_c - ame2d_c | то же, что и ame1r_c, ame2r_c, но при этом все вычисления проводятся с удвоенным числом значащих цифр. При этом параметры a, t, e, r, r1, r2 должны иметь тип double. |
Вызываемые подпрограммы
| utam10_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы ame1r_c. |
| utam11_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы ame1d_c. |
Замечания по использованию
При использовании подпрограмм ame1r_c и ame1d_c значение
параметра n рекомендуется задавать таким, чтобы
погрешность аппроксимации матричной экспоненты exp (AH)
усеченным степенным рядом
14
∑ (AH)k / k!
k=0
то есть величина
∞
exp(AH) - F = ∑ (AH)k / k!
k=15
была достаточно малой.
Таким образом, выбором параметра n можно управлять
точностью вычисления матричной экспоненты exp (AT).
Значение параметров a, t, m, n при работе подпрограмм
сохраняются.
|
1)
| -1 3 0 0 |
A = | 4 -2 0 0 | , T = 1
| 0 0 -3 3 |
| 0 0 4 -2 |
Вычисление матричной экспоненты производится с помощью подпрограммы ame1d_c при значении параметра n = 16
int main(void)
{
/* Initialized data */
static double a[16] /* was [4][4] */ = { -1.,4.,0.,0.,3.,-2.,0.,
0.,0.,0.,-3.,4.,0.,0.,3.,-2. };
/* Local variables */
static int ierr;
extern int ame1d_c(double *, double *, int *, int *, double *,
double *, double *, double *, int *);
static double e[16] /* was [4][4] */;
static int m, n, i__;
static double r__[16] /* was [4][4] */, t,
r1[16] /* was [4][4] */,
r2[16] /* was [4][4] */;
t = 1.;
m = 4;
n = 16;
ame1d_c(a, &t, &m, &n, e, r__, r1, r2, &ierr);
printf("\n %5i \n", ierr);
for (i__ = 0; i__ <= 14; i__+= 2) {
printf("\n %24.14e %24.14e \n", e[i__], e[i__+1]);
}
return 0;
} /* main */
Результаты:
ierr = 0
первая строка
4.225205462389 + 000
3.163850636542 + 000
0.000000000000 + 000
0.000000000000 + 000
вторая строка
4.218467515389 + 000
3.170588583541 + 000
0.000000000000 + 000
0.000000000000 + 000
третья строка
0.000000000000 + 000
0.000000000000 + 000
1.166394356298 + 000
1.163915604121 + 000
четвертая строка
0.000000000000 + 000
0.000000000000 + 000
1.551887472161 + 000
1.554366224338 + 000
2) Матрица A и скалярная величина T те же, что и в примере 1. Матричная экспонента exp (AT) вычисляется с помощью подпрограммы ame2r_c.
int main(void)
{
/* Initialized data */
static float a[16] /* was [4][4] */ = { -1.f,4.f,0.f,0.f,3.f,-2.f,0.f,
0.f,0.f,0.f,-3.f,4.f,0.f,0.f,3.f,-2.f };
/* Local variables */
extern int am11r_c(float *, float *, float *, int *),
ame2r_c(float *, float *, int *, float *, float *,
float *, float *);
static float e[16] /* was [4][4] */;
static int m, i__;
static float r__[4], t, e1[16] /* was [4][4] */,
r1[16] /* was [4][4] */,
r2[16] /* was [4][4] */;
t = 1.f;
m = 4;
ame2r_c(a, &t, &m, e, r__, r1, r2);
for (i__ = 0; i__ <= 14; i__+= 2) {
printf("\n %16.7e %16.7e \n", e[i__], e[i__+1]);
}
t = -1.f;
ame2r_c(a, &t, &m, e1, r__, r1, r2);
for (i__ = 0; i__ <= 14; i__+= 2) {
printf("\n %16.7e %16.7e \n", e1[i__], e1[i__+1]);
}
am11r_c(e, e1, r__, &m);
for (i__ = 0; i__ <= 14; i__+= 2) {
printf("\n %16.7e %16.7e \n", e1[i__], e1[i__+1]);
}
return 0;
} /* main */
4.22520546235 + 00 4.21846751535 + 00 0.00000000000 + 00
3.16385063651 + 00 3.17058858350 + 00 0.00000000000 + 00
0.00000000000 + 00 0.00000000000 + 00 1.16639435628 + 00
0.00000000000 + 00 0.00000000000 + 00 1.16391560411 + 00
0.00000000000 + 00
0.00000000000 + 00
1.55188747214 + 00
1.55436622432 + 00
3) Матрица A та же, что и в первых двух примерах, T = -1. Матричная экспонента exp (AT) вычисляется с помощью подпрограммы ame2r_c.
6.36829740628 + 01 -8.47301850395 + 01 0.00000000000 + 00
-6.35476387796 + 01 8.48655203223 + 01 0.00000000000 + 00
0.00000000000 + 00 0.00000000000 + 00 2.30688401754 + 02
0.00000000000 + 00 0.00000000000 + 00 -1.72740391735 + 02
0.00000000000 + 00
0.00000000000 + 00
-2.30320522314 + 02
1.73108271176 + 02
Произведение матричных экспонент, вычисленных в примерах 2 и 3:
1.00000000000 + 00 0.00000000000 + 00 0.00000000000 + 00
-2.32830643653 - 09 9.99999997672 - 01 0.00000000000 + 00
0.00000000000 + 00 0.00000000000 + 00 9.99999998137 - 01
0.00000000000 + 00 0.00000000000 + 00 -2.32830643653 - 10
0.00000000000 + 00
0.00000000000 + 00
-1.86264514923 - 09
9.99999999069 - 01