Текст подпрограммы и версий ame1r_p.zip , ame1e_p.zip , ame2r_p.zip , ame2e_p.zip |
Тексты тестовых примеров tame1r_p.zip , tame1e_p.zip , tame2r_p.zip , tame2e_p.zip |
Вычисление матричной экспоненты exp (AT) при большом по модулю значении скалярной переменной T.
Вычисляется функция от матрицы exp (AT) - матричная экспонента от произведения матрицы А и скалярной переменной T при большом по модулю значении скалярной величины T по формуле:
exp(AT) = ( exp(AT/N) )N ,
где N - некоторое целое положительное число, задаваемое при обращении к подпрограмме. Матричная экспонента exp (АТ/N) = exp (AH) аппроксимируется усеченным степенным рядом
14 F = ∑ (AH)k / k! k=0
Значение N предполагается таким, что погрешность аппроксимации exp (АН) - F является достаточно малой.
procedure AME1R(var A :Array of Real; T :Real; M :Integer; N :Integer; var R1 :Array of Real; var R :Array of Real; var E :Array of Real; var R2 :Array of Real; var IERR :Integer);
Параметры
A - | вещественный двумерный массив размера М на М, содержащий матрицу, входящую в матричную экспоненту; |
T - | скалярная величина, входящая в матричную экспоненту (тип: вещественный); |
M - | порядок матрицы А (тип: целый); |
N - | некоторое целое положительное число, участвующее в алгоритме вычисления матричной экспоненты. Значение N рекомендуется задавать таким, чтобы погрешность аппроксимации матричной экспоненты exp (AH) усеченным степенным рядом F была достаточно малой; |
E - | вещественный двумерный массив размера М на М, содержащий матрицу - результат exp (AT); |
R - | вещественный одномерный рабочий массив длины М; |
R1, R2 - | вещественные двумерные рабочие массивы размера М на М; |
IERR - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если значение параметра N меньше 1. В этом случае вычисление матричной экспоненты прекращается. |
Версии
AME2R - |
вычисление матричной экспоненты exp (AT) при большом по модулю значении скалярной переменной T по формуле exp(AT) = ( exp(АТ/N) )N в случае, когда N не задается при обращении к подпрограмме, а выбирается в ней автоматически из условия, чтобы || АТ/N || < 1, а именно: exp(AT) = ( exp(AT / [ || А || * | Т | + 1 ] ) ) [ || А || * | Т |+1 ] ,где [Х] - целая часть величины X, а || А || означает максимальную сумму модулей элементов матрицы А по строкам. Первый оператор подпрограммы имеет вид: procedure AME2R(var A :Array of Real; T :Real; M :Integer; var E :Array of Real; var R :Array of Real; var R1 :Array of Real; var R2 :Array of Real);. Число параметров подпрограммы АМЕ2R на 2 меньше числа параметров подпрограммы АМЕ1R, при этом параметры подпрограммы АМЕ2R имеют тот же смысл, что и одноименные параметры АМЕ1R. |
AME1E - AME2E | то же, что и АМЕ1R, АМЕ2R, но при этом все вычисления проводятся с удвоенным числом значащих цифр. При этом параметры А, T, E, R, R1, R2 должны иметь тип Extended. |
Вызываемые подпрограммы
UTAM10 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы АМЕ1R. |
UTAM11 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы АМЕ1E. |
Замечания по использованию
При использовании подпрограмм АМЕ1R и АМЕ1E значение
параметра N рекомендуется задавать таким, чтобы
погрешность аппроксимации матричной экспоненты exp (AH)
усеченным степенным рядом
14 ∑ (AH)k / k! k=0 то есть величина ∞ exp(AH) - F = ∑ (AH)k / k! k=15 была достаточно малой.Таким образом, выбором параметра N можно управлять точностью вычисления матричной экспоненты exp (AT). Значение параметров A, T, М, N при работе подпрограмм сохраняются. |
1) | -1 3 0 0 | A = | 4 -2 0 0 | , T = 1 | 0 0 -3 3 | | 0 0 4 -2 |
Вычисление матричной экспоненты производится с помощью подпрограммы АМЕ1E при значении параметра N = 16
Unit tame1e_p; interface uses SysUtils, Math, { Delphi } LStruct, Lfunc, UtRes_p, AME1E_p; function tame1e: String; implementation function tame1e: String; var M,N,_i,IERR :Integer; T :Extended; A :Array [0..15] of Extended; E :Array [0..15] of Extended; R :Array [0..15] of Extended; R1 :Array [0..15] of Extended; R2 :Array [0..15] of Extended; begin Result := ''; { прототип оператора DATA на FORTRANе } A[0] := -1.0e0; A[1] := 4.0e0; A[2] := 0.0e0; A[3] := 0.0e0; A[4] := 3.0e0; A[5] := -2.0e0; A[6] := 0.0e0; A[7] := 0.0e0; A[8] := 0.0e0; A[9] := 0.0e0; A[10] := -3.0e0; A[11] := 4.0e0; A[12] := 0.0e0; A[13] := 0.0e0; A[14] := 3.0e0; A[15] := -2.0e0; T := 1.e0; M := 4; N := 16; AME1E(A,T,M,N,E,R,R1,R2,IERR); Result := Result + Format('%s',[' IERR=']); Result := Result + Format('%5d ',[IERR]) + #$0D#$0A; Result := Result + #$0D#$0A; for _i:=0 to 15 do begin Result := Result + Format('%20.16f ',[E[_i]]); if ( ((_i+1) mod 2)=0 ) then Result := Result + #$0D#$0A; end; Result := Result + #$0D#$0A; UtRes('tame1e',Result); { вывод результатов в файл tame1e.res } exit; end; end. Результаты: 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) вычисляется с помощью подпрограммы АМЕ2R.
Unit tame2r_p; interface uses SysUtils, Math, { Delphi } LStruct, Lfunc, UtRes_p, AME2R_p, AM11R_p; function tame2r: String; implementation function tame2r: String; var M,_i :Integer; T :Real; E :Array [0..15] of Real; R :Array [0..3] of Real; R1 :Array [0..15] of Real; R2 :Array [0..15] of Real; E1 :Array [0..15] of Real; const A :Array [0..15] of Real = ( -1.0,4.0,0.0,0.0,3.0,-2.0,0.0,0.0,0.0,0.0,-3.0, 4.0,0.0,0.0,3.0,-2.0 ); begin Result := ''; T := 1.0; M := 4; AME2R(A,T,M,E,R,R1,R2); Result := Result + #$0D#$0A; for _i:=0 to 15 do begin Result := Result + Format('%20.16f ',[E[_i]]); if ( ((_i+1) mod 2)=0 ) then Result := Result + #$0D#$0A; end; Result := Result + #$0D#$0A; T := -1.0; AME2R(A,T,M,E1,R,R1,R2); Result := Result + #$0D#$0A; for _i:=0 to 15 do begin Result := Result + Format('%20.16f ',[E1[_i]]); if ( ((_i+1) mod 2)=0 ) then Result := Result + #$0D#$0A; end; Result := Result + #$0D#$0A; AM11R(E,E1,R,M); Result := Result + #$0D#$0A; for _i:=0 to 15 do begin Result := Result + Format('%20.16f ',[E1[_i]]); if ( ((_i+1) mod 2)=0 ) then Result := Result + #$0D#$0A; end; Result := Result + #$0D#$0A; UtRes('tame2r',Result); { вывод результатов в файл tame2r.res } exit; end; end. Результаты: 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) вычисляется с помощью подпрограммы АМЕ2R.
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