CC BY
12
УДК 519.614, 519.615
Методы коррекции последовательной линейной интерполяции для нахождения нулей функций и характеристического полинома матриц специального вида
В.И. Иордан
Алтайский государственный университет (Барнаул, Россия)
Sequential Linear Interpolation Correction Methods Used for Finding Zeros of Functions and Characteristic Polynomials of Matrices of a Special Form V.I. Jordan
Altai State University (Barnaul, Russia)
Рассматриваются вычислительные схемы метода коррекции последовательной линейной интерполяции (МКПЛИ) для нахождения нулей нелинейных (в т.ч. трансцедентных) функций, а также нулей характеристического полинома таких матриц специального вида, как почти треугольные (форма Хессенберга), трехдиагональные и другие виды матриц, получаемые, например, методами Гивенса или Хаусхолдера из матриц общего вида. Предлагаемые вычислительные схемы МКПЛИ для случаев простого и кратного корней (в т.ч. для патологически близких корней) имеют структурно-функциональную общность. Схемы МКПЛИ, предназначенные для локализации и уточнения кратного корня, можно использовать и для локализации группы близких между собой корней, состоящей из простых корней и корней различной кратности (в т.ч. и патологически близких корней). Схемы МКПЛИ обладают устойчивостью вычислений и высокой скоростью сходимости (порядок скорости сходимости приблизительно равен двум). По результатам вычислительных экспериментов для МКПЛИ и других эффективных методов получены зависимости времени диагонализации матриц специального вида от порядка этих матриц.
Ключевые слова: нули функций, характеристический полином, матрицы специального вида, простые и кратные корни, патологические близкие корни.
DOI 10.14258/izvasu(2018)1-16
Computational schemes of the method of "correction of sequential linear interpolation (MCSLI)" are considered in this paper. MCSLI are used for finding zeros of nonlinear (including transcendental) functions, as well as zeros of characteristic polynomials of matrices of a special form, such as almost triangular (Hessenberg form), tridiagonal and others forms of matrices obtained, for example, by Givens or Householder methods from matrices of general form. The proposed computational schemes of MCSLI for cases of simple and multiple roots (including pathologically close roots) have a structural and functional similarity. MCSLI schemes designed to localize and improve multiple roots can also be used to localize a group of closely related roots consisting of simple roots and roots of different multiplicity (including pathologically close roots). The schemes of MCSLI have computational stability and a high convergence rate (the order of the convergence rate is approximately equal to two). Based on the results of computational experiments for MCSLI and other effective methods, the dependences of diagonalization time of matrices of a special form on the order of these matrices are obtained.
Keywords: zeros of functions, characteristic polynomial, matrices of a special form, simple and multiple roots, the pathologically close roots.
Введение
К матрицам специального (компактного) вида можно отнести почти треугольные (форма Хессенберга), трехдиагональные, ленточные и другие виды матриц, получаемых, например, преобразованиями Гивенса или Хаусхолдера из матриц общего вида [1-9]. Характеристическое уравнение матрицы
компактного вида А записано ниже (Е — единичная диагональная матрица)
/(х) = ае1(А - х • Е) = 0, (1)
Подразумевая под /(х) характеристический полином матрицы А либо какую-либо нелинейную функ-
цию, нуль/(х), другими словами, корень с (простой или кратный) уравнения (1) можно находить последовательной полиномиальной интерполяцией т-го порядка: для линейной интерполяции (метода секущих) т = 1, для метода Мюллера (метода парабол) т = 2 [1, с. 385]. А именно, на каждом шаге итерационного процесса нахождения (п + 1)-го приближения х к корню функцию /(х) аппроксимируем полино-
1. Математическое обоснование метода коррекции последовательной линейной интерполяции
Идея предлагаемых схем коррекции метода последовательной линейной интерполяции (метода секущих) [10] состоит в том, что приближение хп+1 метода секущих с помощью схемы коррекции смещается в сторону точного значения корня путем введения в итерационную схему метода секущих некоторых па-
мом т-й степени по (т + 1) приближениям в точках раметров, обеспечивающих сохранение линейности
(х , /(х )), (х ,, /(х ,)), ..., (х , /(х )). Такие методы
4 п^4 п// ч п-1 ^ 4 п-1'' 4 п-т^4 п-т' '
особенно эффективны при определении с высокой
параметрического метода коррекции последовательной линейной интерполяции (МКПЛИ). Поясним это,
точностью определенной части корней (в пределах обозначив простой корень как с.
до 40% от их полного количества, для матриц характеристические корни являются их собственными значениями [1, 2]).
В окрестности простого корня, /(с) = 0 и / '(с) ^ 0, разложим в ряды Тейлора значения /(хп-1) и /(хп), которые удобнее переобозначить через уп-1 и уп:
Уп-1 = У(1) (с)-(хп_1 - с)+ у(2) (с)-(хп_1 - с)2 / 2!+... + у(к) (с)-(хп_1 - с)* / к!+.. Уп = У(1)(с)-(хп- с)+у(2)(с).(хп- с)2 / 2!+... + у(к)(с).(хп - с)к / к!+...
(2)
где у(к)(с) — значение к-й производной в точке с. чив а1 = у(2)(с)/[2!-у(1)(с)] и а2 = у(3)(с)/[3!-у(1)(с)] , Отбросив в (2) члены, начиная с к > 4, обозна- имеем:
уп/ уп-1 —
(хп - с) [у(1) (с)+У2' (с) • (хп - с) / 2!+ у(3) (с) • (хп - с)2 / 3!]
(хп-1 - с) [у(1) (с)+у(2) (с)-(хп-1 - с) / 2!+ у(3) (с)-(хп-1 - с)2 / 3!]
(хп - с) _ [1 + а1 • (хп - с) + а2 • (хп - с)2 ] (хп-1 - с) [1 + а1 • (хп-1 - с) + а2 • (хп-1 - с)2 ]
(3)
Разности (хп-1 - с) и (хп - с) малы, и для слу- 1/(1 + 2) = 1- 2 + 22 - 23 +.....+ (-1)к • 2к +.....,
чая простого корня (у(1)(с)^0) величина или при отбрасывании малых членов, начиная со сте-
2=а1 -(хп-1 -с) +а2 • (хп-1 -с)2 также достаточно мала. пени (хп-1 - с)3 и выше, имеем равенство
Тогда в (3) разложение члена 1/(1 + г) в ряд Маклорена (в точке 2 = 0) имеет вид:
1 / (1+2) — 1-2 + 22 — 1 - а • (хп-1 - с) - а - а12) • (хп-1 - с)2,
позволяющее от (3) перейти к выражению
(хп - с)
уп /уп-1-/
[1 + а1 • (хп -хп-1)+а2 • (хп-хп-1)^(хп + хп-1 -2^с) +
"(хп-1 - с)
+а1 {х- -с)2 -а2 •(хп -с)(хп-1 -с)-а1 • а2 •(хп -с)-^- -с)^ +хп-1 -2с) ]
Отбросив в (5) степени (хп-1 - с)2, используя вспомогательные равенства
|х ,- с = (х - с) - (х - х ,),
"-1 4 п ' 4 п п-1''
хп + хп-1 - 2 • с = 2 • (хп - с) - (хп - хп-1)
обозначив а'— а1 + (2 • а2 - а^) • (хп - с), а' —-а2 +а^ + (а3 - а1 а2) • (хп - с), получим
(хп- с)
уп /уп-1-
(хп-1 - с)
•[1+а' (хп-хп-1)+а^ (хп-хп-1)2].
(4)
(5)
(6)
Если в (6) принять а' = а' = 0, тогда уп/уп-1 — (хп - с)/(хп-1 - с) преобразуется в итерационную схему метода секущих:
(х - х ,)
4 п_п-1'
(уп - уп-1)
уп.
(7)
Соответственно выражение (6) преобразуется в итерационную схему
(х — х ,) • у
4 п п—1' У п
Учитывая в (2) члены более высокого порядка, получим аналогично схеме (8) общую схему (9) с дополнительными членами (хп — хп—1)', где 3 < г < М и а'= 1:
Уп — Уп—1 •[1 + а' • (хп — хп—1) + а2 • (хп — хп—1)2 ]
(х , . — х ) • у ,
4 п—1—1 п' У п—1-
(8)
хп+1 х п—1—1
Уп—1—) — Уп '£а-М) • (хп—1—) — хп )г (10)
г=0
а0М) = 1; 1 = 0,1,..., М.
(х — х ,) • у
4 п п—1' У п
У — У , а'• (х — х ,)'
Уп У п—1 г 4 п п—1'
(9)
Последовательно исключая из каждого уравнения системы (10) неизвестные параметры а(М) , где г = 1,2,...,М, получим полиномиальное уравнение степени (М+1) с одним неизвестным хп+1. Т.е., система Как видно из (3), параметры а, следовательно, (10) эквивалентна методам полиномиальной интерпо-и {а/: г = 1,2,...,М}, определяются через производные ляции [1]. Подтвердим это утверждение для М = 1. различных порядков в точке с. Однако можно избе- А именно система (10) будет состоять из двух уравнений: жать вычисления производных, если использовать уже известные значения уп,уп—1,уп—2,...,уп—1—М и хп, хп—1, хп—2,..., хп—1—М. Определим оценки параметров {аг : г=1,2,...,М}, обозначенные как {а(М): г=1,2,...,М}, следующим образом.
В первом подходе, соблюдая некоторую аналогию с (9), запишем
(х , — х ) • у ,
4 п—1 п' У п—1
Уп—1 — Уп [1 + а1(1) • (хп—1 — хп)]
(хп—2 — хп ) • Уп—2 Уп—2 — Уп [1 + а1(1) • (хп—2 — хп)]
Выразим <51(1) из первого и второго уравнений и приравняем их между собой, последующие преобразования позволяют записать:
[Уп—1 • (хп — хп+1 ) — Уп • (хп—1 — хп+1 )] • (хп—2 — хп + 1 ) • (хп—2 — хп ) = = [Уп—2 • (хп — хп+1 ) — Уп • (хп—2 — хп + 1 )] • (хп—1 — хп+1 ) • (хп—1 — хп )
Последнее уравнение легко преобразуется в уравнение
[Уп • (хп—1 —хп—2 ) + Уп—1 • (хп—2 —хп ) + Уп—2 • (хп —хп—1 )] • х2+1 — —[Уп •(х2— 1 —хп2—2)+ Уп—1 •(хп2—2 —х2)+ Уп—2-(хп2 — ^1 )]хп + 1 +
+[Уп -хп—1 'хп—2 •(хп—1 —хп—2)+ Уп—1 ^п хп—2 ^п—2 —хп) +
+Уп—2 'хп—1 •(хп —хп—1)] = 0
(11)
Положив К = х{ — хг—1, = К / к—1, 8{ = 1 + , [1, с. 385], представляющего собой метод квадратич-Як = /к—2 • К — /к—1 • + /к • Л + ¿к) , легко прове- ной интерполяции: рить, что (11) сводится к уравнению метода Мюллера
Ак2+1 Л А—1■ •( /к—2 Л — /к—1 ¿к + /)+Лс+1 А—1 •Як + / = 0.
Во втором подходе общности с (9) еще больше. Запишем систему:
(х — х , ) • У
4 п п—1— у У п
У — У , .•'^а(М) • (х — х , У
Уп Уп—1—1 г 4 п п—1—у
г=0
а0М) = 1; 1 = 0,1,..., М
Произведя исключение хп+1 из всех (М+1) уравнений системы (12), получим систему уравнений, линейную по отношению к {й(М): ¿=1,2,...,М}:
Уп—У-• £а(М) -(хп -хп-1)'' Уп-Уп—1—,' £а(М) -(хп Хп—1—,У
Хп Хп—1—1
а0м) = 1; ) = 0,1,...,М.
(13)
От системы (13) можно перейти к матричной форме В Л = С :
£Ь, • а(М) = с,; ¿,] = 1,2,...,М.
где
Ь.. = у (х —х , )• (х —х ,)¿ — у , • (х —х ,)• (х —х , .)',
}¿ у п— 1 4 п п— 1 — ]' 4 п п — 1 У п— 1 —} 4 п п — 1 4 п п —1—} ' 7
с1 = Уп • (хп—1—хп—1—,) + Уп—1 • (хп—1—1 —хп) + Уп—1—1 • (хп—хп—1),
¿, 1 = 1,2,..., М
(14)
(15)
Система (14) является линейной и имеет единственное решение, т.к. ёе^В) ^ 0 в силу линейной независимости строк (столбцов) матрицы В.
Определив решение Л = {<а((М): ¿=1,2,..., М} из (14) и (15), достаточно просто вычислить новое приближение хп+1 из схемы (9), которая определяет основу метода МКПЛИ. Схема (9) не является интерполяционной в чистом виде, т.к. лишь одна точка (хп, уп ) оказывается интерполяционной (рис.). Введем скорректированное для уп—1 значение
Уп—1 = Уп—1 £
¿=0
М
М)
а
соединяющей точки (хп—1, у>п—1) и (хп,уп), обеспечивает более эффективное сближение к корню с нового приближения с = хп+1 по сравнению с приближением с' = хп+1, получаемым по «обычной» схеме (7) метода секущих (рис.). В каждой итерации уп и уп выбираются с противоположными знаками (рис.), и тогда корень с будет лежать между хп—1 и хп, тем самым обеспечивая устойчивость сходимости схемы (9).
• (хп — хп—1)'. Тогда наклон секущей,
Графическая схема МКПЛИ
Для М = 1 решение системы (14) с учетом (15) определяется просто:
ГД1)" гУп Уп—1 Уп Уп—2Ч/(Уп—1 — Уп—2)
а = I
х„ — х„
х„ — х„
п—1 п п—2
Конкретной итерационной схемой (М = 1), использующей (16), будет схема
(хп — хп—1) • Уп
Уп — Уп—1 •[! + а1(1) • (хп — хп—1)]
(16)
Для М = 2 искомыми будут два параметра: а\ и а2 . Введя обозначения
(Дх
= х , — х , Ду
я—к п—1 , У я
= у,
п—к У п—1 у
1 у' , , = Ду , , / Дх , ,; к = 0,1,2,3; 1 = 0,1,2,3
п — к ,п—1 У п—к ,п — 1 п—к ,п—1 ' > > > :> :»:»:»
и получив решение системы (14) и (15), запишем
И2) =а(1) +[ Уп—2/ ^п—1,п—2 — Дхп,п—11-й(2)
а(2) =
ДУп—1,п—2 [ У ] — ДУп—1,п—3 [У
у • у . - Дх , . + у ■ у _ - Дх _ , + у • у _ - Дх . _
✓ п—1 у п—2 п—1,п—2 ✓ п—1 У п—3 п—3,п—1 ✓ п—2 У п—3 п—2,п—3
(18)
где а ) определяется по формуле (16). Следовательно, вторая итерационная схема с параметрами а( и а (М = 2) имеет вид:
У2(2)
(х — х ,) - у
4 я я—1' ✓я
уя—у—1 -[1+а;2)-х—хя—1)+а22) -(хя —хя—1 )2 ]
(19)
В методе «секущих», а также при выводе итераци- сти т. Для кратных корней предлагается использо-онных схем (17) и (19) метода МКПЛИ производится вать схемы, подобные схемам (9). Для корня кратно-
уточнение простого локализованного корня.
Допустим же, что в процессе деления «спектра» корней (локализации корней) имеем корень кратно-
сти т верно соотношение
Уп / Уя—1
(х„ — с)т [у(т) (с)/ т!+ у(т+1) (с) - (хп — с) / (т +1)!+......]
(хп—1 — с)т [у(т) (с)/ т!+у(т+1) (с)-(хп—1 — с)/(т +1)!+.....]
так как у(1) (с) = у(2) (с) =.......= у(т 1) (с) = 0 для корня
кратности т. Или
уп / уп—1
(хп — с)т [1 + Д - (хп — с) +......]
(хп—1 — с)т [1 + Д -(хп—1 — с) +.....]
где Д = [у(т+0(с)-т!]/[у(т)(с)-(т + г)!] и I = 1,2,.....
Пренебрегая степенями при коэффициентах Д , получим
уп /уп—1 = (хп — с)т/(хп—1 — с)т. (21)
Если т — нечетное, тогда, обозначив уп = (уп )1/т и Уп—1 = (уп—1)1/т , получим:
Уп / Уп—1 = (хп — с)/(хя—1 — с), (22)
(х — х ,)
4 я п—1'
Уп
(23)
" (^п — ^п—1)
для уточнения кратного корня. Учитывая в (20)
(20) степени при Д , получим схему
(х — х ,) - у
4 я п—1' У я
г—\~ (М)
уя — Уп—1 вг - (хп — хя—1)г
(24)
которая вместе с (23) дает схемы метода коррекции последовательной линейной интерполяции (МКПЛИ) для кратного корня. Схемы (23) и (24) можно использовать не только для уточнения кратного корня, а также для разделения и последующей локализации группы близких простых корней, в т.ч. «патологически близких» корней. Покажем это.
Допустим, что в интервале между хп и хп—1 находится р корней с1,с2,...,ср с кратностью соответствен-
Если т — четное, то уп и уп—1 одного но к1,к2,...,к и выполняется условие
к1 + к2 +.... + к = т ,
(25)
знака. Учитывая знаки разностей (хп — с)
и (хп—1 — с), определим уп = |уп |"т - sign(xn — с)
и уп—1 = |уп—11' - sign(xn—1 — с), чтобы выполнялось (22).
Тогда аналогично вышеприведенным рассуждениям Из (22) получаем подобную методу секущих схему и выражениям (20), (21) можно записать
т.е. у = /(х) = (х — с^1 -(х — с^2 ---(х — с )р -^(х)
уп
уп—1
хя — с1 к1 хп — с2 к2 хя —ср кр
. хя—1 — с1 . . хп—1 — с2 . хп— —ср
1 + 71 - (хп — с) +...
1 + 71 - (хп—1 — с) +..
хп — с
. хя—1 — с ,
где y = • i!] и
4>(xn) 1 + jl • (xn — c) +...
^(x„-1) 1 + Yi • (xn-i - c) +..
. Из (26) следует
[ xn — ci k1 xn — c2 k2 xn — cp kp m 1 + (yJ m) • (xn — c) +..
1 xn—1 — ci . . xn—1 — c2 . xn—1 — cp . 1 + (yJ m) • (xn—1 — c) +
(27)
Л Л-i
x — c
Значение с в выражении (27), аналогичном выражению (22), является следующим приближением с = xn+1 к какому-либо из корней с1,с2,...,ср , определяемого из схем (23), (24). Определим смысл величины с = х ,,.
Если в (27) пренебречь малыми величинами / m) • (хп — с)' и (7. / т) • (хп-1 — с)' в последнем сомножителе левой части уравнения, тогда
xn — c1 k1 xn — c2 k'2 xn —cp
. xn—1 — c1 . . xn—1 — c2 . xn— —V
(28)
С учетом (25) правая часть уравнения (28) представляет собой «среднее геометрическое» значение для величин (xn — c)/(xn—1 — c), где г = 1,2,...,p . Следовательно, значение c = xn+1, соответствующее выражениям (26) - (28) и определяемое по схеме (23) или (24), будет разделять множество корней на два подмножества. Определенное количество итераций МКПЛИ позволяет довести до конца процесс разделения корней, т.е. каждое подмножество будет состоять из одиночного корня, который будет уточняться по схемам (17) и (19). Патологически близкие корни, оказавшиеся в интервале погрешности вычислений на ЭВМ, могут оказаться неразделимыми (их можно считать корнем соответствующей кратности).
2. Результаты вычислительных экспериментов
Погрешность приближения хя+1 к корню с определяется [10] соотношением |хя+1 — c| = G'' • |хя — c|1,
где G — константа, связанная с производными от f(x), а параметр у связан с порядком скорости сходимости р. Значение р получено в [10], его оценка приблизительно равна 2. Разработаны программы на различных языках программирования (Pascal, C/C++ и т.д.),
использующие схемы для уточнения простого корня (17), (19) и схемы (24) для кратного корня, а также свойство последовательности Штурма по разделению корней [1], которые позволяют находить часть корней характеристического полинома матриц специального вида (собственные значения этих матриц) в определенном интервале либо все собственные значения матриц.
Численные эксперименты показали эффективность вычислительных схем МКПЛИ [10]. А именно зависимость времени вычисления всех корней трехдиагональных матриц от порядка матриц N определяется зависимостью t = с • [10] показатель в ~ 1.585, а для наиболее быстродействующего QR-метода [8] показатель в = 2 (время счета больше).
Заключение
В задачах диагонализации матриц специального вида и нахождения корней нелинейных (в т.ч. и транс-цедентных) уравнений могут быть использованы разработанные на основе МКПЛИ программы автора, сохраняющие инвариантной форму матриц специального вида и обеспечивающие абсолютную сходимость и устойчивость вычислений, а также высокое быстродействие, сравнимое с быстродействием QR-метода.
Библиографический список
1. Уилкинсон Дж.Х. Алгебраическая проблема собственных значений / пер. с англ. В.В. Воеводина и В.Н. Фадеевой. — М., 1970.
2. Парлетт Б. Симметрическая проблема собственных значений. Численные методы // пер. с англ. Х.Д. Икрамова и Ю.А. Кузнецова. — М., 1983.
3. Watkins D.S. The Matrix Eigenvalue Problem: GR and Krylov Methods // D.S. Watkins. — SIAM. — 2007.
4. Prodi G. Eigenvalues of non-linear problems // G. Prodi (ed.). — Berlin, 2010.
5. Новиков М.А. Одновременная диагонализация трех вещественных симметричных матриц // Известия вузов. Математика. — 2014. — № 12.
6. Кочура А.Е., Подкользина Л.В., Ивакин Я.А., Нид-зиев И.И. Сингулярные матричные пучки в обобщенной симметричной проблеме собственных значений // Труды СПИИРАН. — 2013. — Вып. 3 (26).
7. Кузнецов Ю.И. Проблема собственных значений симметричной теплициевой матрицы // Сибирский журнал вычислительной математики. — 2009. — Т. 12, № 4.
8. Иордан В.И. Быстродействующие алгоритмы диа-гонализации трехдиагональных симметричных матриц
на основе элементарных плоских вращений // Изв. Алт. гос. ун-та. — 2017. — № 1 (93). 001 10.14258/ izvasu(2017)1-15.
9. Калинина Е.А. Кратные собственные числа матрицы с элементами, полиномиально зависящими от параметра // Вестник СПбГУ Сер. 10. Прикладная математика. Информатика. Процессы управления. — 2016. — Вып. 2. Э01 10.21638/1170ШрЬи10.2016.203.
10. Иордан В.И. Эффективные методы определения энергетического спектра матриц большой размерности в задачах экспериментальной физики : дисс. ... канд. физ.-мат. наук. — Барнаул, 2003.
