Как и прежде, для обеспечения устойчивости требуется, чтобы

I(V-b)J<1-0-

В общем случае, когда нельзя получить явных формул для собственных значений матрицы А-В, (4.29) удобно переписать в виде

(4.30>

где D = A~B.

Максимальные собственные значения матрицы D можно найти, пользуясь степенным методом [Carnahan et al., 1969]. Это потребует расчетов по формуле

Jг+l DX (4.31>

повторяемых до тех пор, пока Х + не совпадет с X . Если сходимость достигнута, то Х + является собственным вектором, соответствующим максимальному собственному значению которое определяется максимальным элементом вектора X +L Обычно на каждой итерации масштаб подбирается так, чтобы максимальный элемент вектора Х равнялся единице. Стартовое значение Х выбирается произвольно. Карнаган, Лютер и Уилкс предлагают программу такого расчета и обсуждают возможные трудности, связанные с применением метода.

Применительно к двухслойной схеме общего вида (4.23) максимальные собственные значения матрицы D = А-В для граничных условий Дирихле и Неймана, поставленных в точке л: = О, приводятся в табл. 4.2 для некоторого набора значений Р и 5. Эти результаты были получены с помощью степеннбга*

Таблица 4.2. Данные о численной устойчивости для двухслойной схемы общего вида

метода. В случае граничного условия Дирихле в точке л: = О ]результаты согласуются с требованием

<-ПГ2р- р<0.5,

обеспечивающим устойчивость (Jimax 1). Ясно, что необходимость выполнения граничного условия Неймана (4.27) несколько уменьшает устойчивость схемы при одинаковых значениях р и S.

4,3,4. Метод Неймана: схема ВВЦП

При определении критериев устойчивости наиболее общепринятым является метод Неймана, так как в общем случае его легче всего применять и так как он самый непосредственный и надежный. К сожалению, однако, этот метод можно применять только для установления необходимого и достаточного условия для линейных задач с начальными условиями и постоянными коэффициентами.

Как правило, практические задачи охватывают и переменные коэффициенты, и различного рода нелинейности, и усложненные варианты граничных условий. А если что-то из этого имеет место, то метод может применяться только локально и при замораживании нелинейных членов. Для такой более общей ситуации метод Неймана обеспечивает необходимые, но не всегда достаточные условия устойчивости. Строго говоря, метод применим только для внутренних точек. Однако метод Неймана может дать эвристическую информацию о том, как влияют граничные условия для производных на численную устойчивость, если этот метод применять по отдельности к тем алгоритмам, которые используются на границах [Тгарр, Ramshaw, 1976].

Метод Неймана состоит в том, что ошибки, распределенные вдоль линий сетки на одном временном слое, разлагаются в конечный ряд Фурье. После этого устойчивость или неустойчивость вычислительного алгоритма определяется в зависимости от того, будут ли отдельные компоненты Фурье для распределения ошибок затухать или возрастать при переходе на следующий временной слой. Итак, вектор начальной ошибки выражается с помощью конечного комплексного ряда Фурье, -так что в точке xj ошибка определяется формулой

tf = Z ciemf / = 2,3,...,/-!, (4.32)

где 1=(-1)/2 и Qm = mnAx. Ряд Фурье (4.32) по виду совпадает с тем, что было в выражении (3.30), если не считать того что здесь распределение ошибок на интересующем нас интервале О л: 1 предполагается периодическим. Если учесть линейность вычислительного алгоритма, достаточно изучить распространение ошибки, обусловленное единственным членом ехр(Шт/) представления с помощью ряда Фурье согласно формуле (4.32). Поэтому в дальнейшем обозначение Вт мы будем писать без индекса т.

С учетом вида начального распределения ошибок решение-ошибочного уравнения (4.17), соответствующего схеме ВВЦП, ищется по методу разделения переменных в форме

J==Gfe, (4.33)

причем в данной формуле зависимость данной компоненты Фурье от времени включается в комплексный коэффициент (G) верхний индекс в символе которого подразумевает, что величина G возведена в степень я Подстановка выражения (4.33) в уравнение (4.17) позволяет получить уравнение

-е/ е(/-1) + (1 2s) (Of е + s [Gf е

После некоторых алгебраических манипуляций получим

G=l-4ssin2(e/2). (4.34)

Величину G можно рассматривать как коэффициент усиления моды Фурье с индексом т для распределения ошибок, соответствующий продвижению на один шаг вперед по времени;, это подтверждается тем, что, как следует из (4.33),

rVl/=G. (4.35)

Можно отметить, что G является функцией 5 и 0. Иначе говоря, G(5, 0) зависит от размера элементов сетки, а также от того, какая именно мода Фурье рассматривается, поскольку 5 = аМ/Ах и 0 = тпАх. Ошибки будут оставаться ограниченными, если абсолютная величина G, называемая также избытком, никогда не будет превышать единицу ни для каких мод Фурье 0. Таким образом, общее требование устойчивости сводится к условию

IGK1 при любых 0. (4.36)

Как следует из (4.34), требование устойчивости для схемы ВВЦП гласит, что

- 1 < 1 - 4s sin2 (0/2) < 1 при любых 0, (4.37)