Текст подпрограммы и версий agj0r_p.zip |
Тексты тестовых примеров tagj0r_p.zip |
Вычисление всех собственных значений и собственных вектоpов в обобщенной проблеме AX = λBX, где A - матрица Якоби, а B - диагональная положительная матрица.
Обобщенная проблема Ax = λBx с якобиевой матрицей А и диагональной положительной матрицей B последовательностью преобразований сводится к стандартной задаче Qz = λz. Здесь Q - симметричная трехдиагональная матрица, собственные значения которой совпадают с собственными значениями исходной обобщенной проблемы, а собственные векторы z связаны с собственными векторами x обобщенной проблемы соотношением x = Cz, где матрица C определяется последовательностью выполненных преобразований.
Для решения стандартной проблемы применяется неявный QL - алгоритм со сдвигом. Собственные векторы z нормированы таким образом, что zTz = E, где E - единичная матрица.
Дж.Х. Уилкинсон, Алгебраическая проблема собственных значений, "Hаука", M., 1970.
procedure AGJ0R(var A :Array of Real; var B :Array of Real; var V :Array of Real; var EV :Array of Real; var RAB :Array of Real; N :Integer; var IERR :Integer);
Параметры
A - | вещественный двумерный массив размера N на 3, содержащий якобиеву матрицу A в компактной форме записи для ленточных матриц; |
B - | вещественный вектоp длины N, содержащий элементы диагональной матрицы B; |
V - | вещественный двумерный массив размера N на N, в столбцах которого содержатся вычисленные оpтонормированные собственные векторы; |
EV - | вещественный вектоp длины N, содержащий вычисленные в возрастающем порядке собственные значения; |
RAB - | вещественный вектоp длины 2 на N, используемый как рабочий; |
N - | порядок исходных матриц (тип: целый); |
IERR - | целая переменная служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; значение IERR полагается равным: |
номеpу собственного значения, для вычисления которого потребовалось более 30 итераций, при этом собственные значения с индексами 1, 2, ..., IERR - 1 и соответствующие собственные векторы вычислены правильно, но собственные значения неупорядочены; IERR = 65 - когда исходная диагональная матрица B не является положительно определенной; IERR = 66 - когда в исходной матрице Якоби не все попарные произведения соответствующих произведений побочных диагоналей неотрицательны; IERR = 67 - когда в матрице Якоби имеется pавное нулю произведение соответствующих элементов побочных диагоналей, причем сомножители равны нулю неодновpеменно; в этом случае нет симметризации, необходимой для правильного вычисления собственных вектоpов; IERR = 68 - когда матрица Q получена, но матрица преобразования C не вычислена, так как она алгоритмически выражена; в этом случае собственные векторы не вычисляются. |
Версии : нет
Вызываемые подпрограммы
AFJ0R - | приведение матрицы Якоби к симметрической трехдиагональной матрице преобразованием подобия. |
AFE1R - | одновременное приведение вещественной симметрической трехдиагональной матрицы и диагональной положительно определенной матрицы двусторонним диагональным преобразованием соответственно к вещественной симметрической трехдиагональной и единичной матрицам. |
UTAG10 - UTAG11 | подпрограммы выдачи диагностических сообщений при работе подпрограммы AGJ0R. |
Замечания по использованию
Подпрограмма AGJ0R сохраняет матрицу A. |
Unit TAGJ0R_p; interface uses SysUtils, Math, { Delphi } Lstruct, Lfunc, UtRes_p, AGJ0R_p; function TAGJ0R: String; implementation function TAGJ0R: String; var J,I,N,_i,IERR :Integer; EV :Array [0..4] of Real; V :Array [0..24] of Real; RАВ :Array [0..9] of Real; const A :Array [0..14] of Real = ( 0.0,1.0,1.0,1.0,1.0,0.0,0.0,0.0,0.0,0.0,4.0,4.0, 4.0,4.0,0.0 ); B :Array [0..4] of Real = ( 2.0,2.0,2.0,2.0,2.0 ); begin Result := ''; { результат функции } Result := Result + #$0D#$0A; Result := Result + Format('%s', [#$0D#$0A + ' ВСЕ СОБСТВЕННЫЕ ЗНАЧЕНИЯ И СОБСТВЕННЫЕ ВЕКТОРЫ ' + #$0D#$0A + ' УРАВНЕНИЯ BИДA:AX=LAMBDA*BX, ГДЕ A - ВЕЩЕСТВЕННАЯ MATPИЦA' + #$0D#$0A + ' ЯКОБИ, B - ВЕЩЕСТВЕННАЯ ДИАГОНАЛЬНАЯ МАТРИЦА ' + #$0D#$0A + ' A' + #$0D#$0A]); for I:=1 to 5 do begin for J:=1 to 3 do begin Result := Result + Format(' %20.16f ',[A[(I-1)+(J-1)*5]]) + #$0D#$0A; end; end; Result := Result + #$0D#$0A; Result := Result + #$0D#$0A; Result := Result + Format('%s',[#$0D#$0A + ' B' + #$0D#$0A]); for I:=1 to 5 do begin Result := Result + Format(' %20.16f ',[B[I-1]]) + #$0D#$0A; end; Result := Result + #$0D#$0A; N := 5; AGJ0R(A,B,V,EV,RAB,N,IERR); Result := Result + Format('%s', [#$0D#$0A + ' PEЗУЛЬTAT' + #$0D#$0A]) + #$0D#$0A; Result := Result + #$0D#$0A; Result := Result + Format('%s',[#$0D#$0A + ' A' + #$0D#$0A]); for I:=1 to 5 do begin for J:=1 to 3 do begin Result := Result + Format(' %20.16f ',[A[(I-1)+(J-1)*5]]) + #$0D#$0A; end; end; Result := Result + #$0D#$0A; Result := Result + #$0D#$0A; Result := Result + Format('%s',[#$0D#$0A + ' B' + #$0D#$0A]); for I:=1 to 5 do begin Result := Result + Format(' %20.16f ',[B[I-1]]) + #$0D#$0A; end; Result := Result + #$0D#$0A; Result := Result + Format('%s',[#$0D#$0A + ' EV' + #$0D#$0A]); Result := Result + #$0D#$0A; for _i:=0 to 4 do begin Result := Result + Format('%20.16f ',[EV[_i]]); if ( ((_i+1) mod 3)=0 ) then Result := Result + #$0D#$0A; end; Result := Result + #$0D#$0A; Result := Result + #$0D#$0A; Result := Result + Format('%s',[#$0D#$0A + ' V' + #$0D#$0A]); for I:=1 to 5 do begin for J:=1 to 5 do begin Result := Result + Format(' %20.16f ',[V[(I-1)+(J-1)*5]]) + #$0D#$0A; end; end; Result := Result + #$0D#$0A; Result := Result + Format('%s',[#$0D#$0A + ' IERR' + #$0D#$0A]); Result := Result + Format('%8d ',[IERR]) + #$0D#$0A; UtRes('TAGJ0R',Result); { вывод результатов в файл TAGJ0R.res } exit; end; end. Результаты: | -1.7320508 | | 0.2041241 | | -1.0000000 | | -0.1767767 | EV = | 0.0000000 | , V1 = | 0.1020621 | , | 1.0000000 | | -0.0441942 | | 1.7320508 | | 0.0127578 | | -0.3535534 | | -0.4082483 | | 0.1767767 | | -0.0000000 | V2 = | 0.0000000 | , V3 = | 0.1020621 | , | -0.0441942 | | 0.0000000 | | 0.0220971 | | -0.0255155 | | -0.3535534 | | -0.2041241 | | -0.1767767 | | -0.1767767 | V4 = | -0.0000000 | , V5 = | -0.1020621 | , | 0.0441942 | | -0.0441942 | | 0.0220971 | | -0.0127578 | IERR = 0.