1626435587-55f52a4de97976f3c6215fa7c103f544 (Смирнов 2017 - Основы вычислительной физики ч2), страница 4

Выражение в правой части (30) есть многочлен степени от , называемый.Как решать характеристическое уравнение (30)? Даже вычислениекоэффициентов характеристического полинома с использованием хорошо известных формул для определителя требует (!) арифметических операций, что делает данный метод совершенно непригодным при ≥ 15 (и крайне медленным при меньших ). Вместе с тем, как ужеотмечалось ранее, в физических расчётах возникают матрицы порядка > 103 или даже 104 .

Ниже мы познакомимся с идеей нескольких наиболее простых методов решения спектральной задачи и свойственнымиим ограничениями.характеристическим полиномом1.7.1. Метод интерполяцииПростейшее эффективное решение задачи (29) известно как методинтерполяции. Выберем произвольным образом + 1 точку 0, . .

. , и вычислим в них характеристический полином: ≡ ( ), =0, . . . , . Затем построим интерполяционный полином ˜ степени ,проходящий через заданные точки (0 , 0 ), . . . , ( , ) [1, гл. 5]. Поскольку через + 1 точку можно провести один и только один полиномстепени , имеем () ≡ ˜ () ∀. Следовательно, описанный вышеалгоритм позволяет вычислять коэффициенты полинома ˜ , отличающегося от только погрешностью вычислений. Зная коэффициентыполинома , можно достаточно эффективно вычислить корни этогополинома — например, используя метод парабол [2, с.

146].Оценим количество арифметических операций, необходимых дляиспользования метода интерполяции. Поскольку вычисление ( )18есть нахождение определителя матрицы − с известными коэффициентами, оно требует ≈ 23 3 операций (см. п. 1.4 на с. 13). Вычисление значений полинома в + 1 точке потребует ≈ 23 4 операций. Построение интерполяционного полинома требует всего (2 ) операций;считая, что метод парабол сходится за 10 итераций [2, с. 165], вычислительной сложностью поиска корней полинома также можно пренебречьпри > 20.Поскольку метод интерполяции основан на вычислении конечныхразностей [1, п.

5.2], его применимость ограничена матрицами сравнительно небольшого порядка ( ≤ 10 . . . 20). Кроме того, к числунедостатков данного метода можно отнести относительно высокуюсложность его программирования, поскольку метод включает в себянесколько подзадач: вычисление определителей, построение интерполяционного полинома, поиск корней.1.7.2. Степенной методЗачастую интерес представляют не все решения (29), но лишь одно или несколько наибольших (или, наоборот, наименьших) собственных значений. Например, в телекоммуникациях и нелинейной оптике наибольший интерес представляет основная или несколько первыхмод оптических волокон, поскольку высоковозбуждённые моды имеют, как правило, высокие потери. Аналогично, в задачах квантовоймеханики зачастую бывает необходимо найти лишь основное и первыевозбуждённые квантовые состояния.

Кроме того, при дискретизацииˆспектральной задачи ()= () и переходе к её конечномерному (матричному) аналогу (29) погрешность ответа быстро возрастаетс числом нулей волновой функции (). Вследствие этого лишь относительно небольшая часть численных решений — основное и несколькопервых возбуждённых состояний — имеют отношение к исходной задаче и несут физический смысл. В таких случаях вместо полного решения проблемы собственных значений оказывается выгоднее использовать тот или иной итерационный процесс, сходящийся к одному собственному значению и одному собственному вектору.

Один из наиболеепростых примеров таких итерационных процессов, позволяющих найти решение, известен какили.Перенумеруем собственные числа симметричной матрицы в порядке убывания модуля:частичной проблемы собственных значенийстепенной метод счёт на установление|1 | > |2 | ≥ |3 | ≥ . . . ≥ | |.19(31)В качестве начального приближения для итерационного процесса выберем произвольным образом векторu(0) = 1 u1 + 2 u2 + .

. . ,(32)где uj — -й собственный вектор матрицы: uj = uj . На каждой итерации будем умножать вектор на матрицу:u() = u(−1) .(33)Покажем, что итерационный процесс (33) сходится к собственному вектору u1 , соответствующему наибольшему по модулю собственному значению 1 . Используя начальные условия (32), на -й итерации процесса(33) будем иметь:⎛⎞(︂ )︂∑︁∑︁u() = uj = 1 ⎝1 u1 +uj ⎠ .(34)1=1=2Здесь мы вынесли 1 , с тем чтобы коэффициент при u1 в разложенииu() был равен 1 — так же, как и в начальных условиях (32). Приэтом видно, что коэффициенты при uj ( = 2, .

. . , ) содержат малыевеличины вида ( /1 ) , которые стремятся к нулю при → ∞ в силунеравенств (31). Следовательно, в пределе → ∞ направление вектораu() будет совпадать с направлением искомого вектора u1 . Обрываяпроцесс (33) на -м шаге, получаем в результате искомый вектор u1с погрешностью ((2 /1 ) ). Из (33) видно, что для вычисления 1достаточно найти отношение произвольной компоненты вектора u надвух последовательных итерациях:()()(−1)1; = /+ ((2 /1 ) ).(35)Итерационный процесс (33) можно останавливать, когда несколько значений 1 , вычисленных по формуле (35) с использованием разных компонент вектора u() , совпадают в требуемом числе знаков.Используя в (35) норму векторов вместо фиксированной ( -й) компоненты u(k) , можно вдвое ускорить сходимость метода:√︃(u(k) , u(k) )()+ ((2 /1 )2 ).(36)1 =(k−1)(u, u(k−1) )Заметим, что если начальное приближение (32) было выбрано так,что 1 ≈ 0, то, несмотря на малость множителей ( /1 ) , основной20∑︀вклад в сумму uj (34) при недостаточно больших значениях может давать не первое, а второе или последующие слагаемые4 .

С одной стороны, данное обстоятельство может приводить к ошибочномурезультату; для предотвращения ошибки следует выбирать вектор u(0)несимметричным, либо его симметрия должна соответствовать симметрии искомого ответа5 . С другой стороны, указанное обстоятельство может быть использовано для нахождения второго и последующих собственных значений и соответствующих им собственных векторов. Дляэтого нужно выбрать начальное приближение u(0) ортогональным найденным ранее векторам u1 , u2 , .

. ., вычитая из произвольно выбранногоначального приближения его проекцию на первые найденных ранеерешений:∑︁(uj , u(0) )u(0) ← u(0) −uj.(37)‖uj ‖2=1Поскольку векторы u1 , u2 , . . . , um были найдены в результате численного решения и известны приближённо, процесс ортогонализации (37)имеет смысл проводить не только перед началом итерационного процесса (33), но также и после его завершения, и в нескольких промежуточных точках.1.7.3. Метод обратных итераций со сдвигомРассмотренный выше степенной метод позволяет искать наибольшее по модулю собственное значение матрицы. Для решения обратнойзадачи — поискапо модулю собственного значения — очевидно, могут быть использованы итерации с заменой матрицы на−1 .

Действительно, если 1 , . . . , — собственные числа матрицы ,−1то матрица −1 будет иметь собственные числа, равные −11 , . . . , .−1При выполнении условия (31), степенной метод сойдётся к , что позволит найти наименьшее по модулю собственное значение матрицы.Выполним теперь ещё одно преобразование матрицы , что позволит нам весьма эффективно (быстро) найти любое собственное значе˜ ≈ .

А именние , для которого известно начальное приближение ˜ и применим степенной метод дляно,матрицу на величину ˜ −1 :( − )наименьшегосдвинемВ пределе → ∞ направление вектора u(k) может совпадать с направлениемдажепри 1 = 0 из-за наличия ошибок округления.5 Например, при поиске основного уровня энергии квантовой частицы в симметричной потенциальной яме () = (−) начальное приближение (0) () должнобыть либо чётной функцией , либо функцией без определённой чётности.4u121(︀)︀˜ −1 u() .u(+1) = − (38)сдвинутойИдея перехода к обратнойматрице в (38) состоит в том,˜ −1 матрицы ( − )˜ −1 будет больши́мчто собственное число ( − )по величине, что обеспечит быструю сходимость степенного метода.В соответствии с (35), (36), сходимость метода (38) будет тем быст˜ , т.

е. чем более точным былорее, чем меньше модуль разности | − |˜начальное приближение . Следовательно, скорость сходимости можноповысить (сделать квадратичной), если в (38) вместо фиксированного˜ использовать приближение () , найденное на предыдущейсдвига ˜итерации. Заметим, однако, что использование переменного сдвига в (38) предполагает вычисление обратной матрицы на каждом шагеитерационного процесса. В этой связи выгоднее переписать (38), до˜ , и вычислять u(+1) ,множив обе части равенства на матрицу ( − )решая систему линейных уравнений:(︀)︀ − () u(+1) = u() ,(+1) = () +‖u() ‖.‖u(+1) ‖(39)Для начала итерационного процесса (39) необходимо выбрать произвольный вектор u(0) и начальное приближение к искомому собственному значению (0) .

Применяя (39) при = 0, находим вектор u(1) изрешения системы линейных уравнений и вычисляем (1) . Далее повторяем процесс (39) в цикле ( = 1, 2, . . .), завершая его по достижениитребуемого уровня точности.1.7.4. Метод НьютонаПосмотрим на спектральную задачу u−u = 0 как на систему из уравнений на +1 неизвестную величину 1 , . . . , , . Система является нелинейной ввиду наличия члена u. Линеаризуем её подстановкамиu = u() +u(+1) , = () +(+1) . Здесь u() и () — известные на -йитерации приближения к искомым собственному вектору и собственному значению, u(+1) и (+1) — малые поправки к ним, которые будутнайдены на ( + 1)-й итерации.