Обозначим буквой D точное решение этого разностного уравнения, т. е. решение, которое мы получили бы на ЭВМ при отсутствии ошибок округления, а буквой Л/ -решение, полученное на реальной ЭВМ. Если А - аналитическое решение исходного уравнения в частных производных, то можно записать

Погрешность аппроксимации = Л - Z),

Погрешность округления =N - D,

Устойчивость конечно-разностной схемы определяется изменением погрешности в процессе вычисления. ОЪрайен и др. [OBrien et al., 1950] предложили следующую классификацию устойчивости разностных схем:

1. Если полная погрешность округления растет {не растет), то разностная схема называется сильно неустойчивой {устойчивой).

2. Если отдельная погрешность округления растет {не растет), то разностная схема называется слабо неустойчивой {устойчивой).

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

3.6.1. Метод Фурье или метод Неймана

Рассмотрим разностное уравнение (3.101). Если е - погрешность округления, то численное решение разностного уравнения можно представить в виде

М = 0 + г. (3.102)

Так как численное решение должно удовлетворять разностному уравнению, то, подставляя (3.102) в (3.101), получим

Точное решение D удовлетворяет разностному уравнению (3.101), поэтому и погрешность е удовлетворяет тому же уравнению:

Так как точное решение разностного уравнения D и погрешность округления 8 удовлетворяют одному и тому же уравнению, то и растут по времени они одинаково. Если разностная схема устойчива, то рост любого возмущения, вводимого на п~м шаге по времени, ограничен; для неустойчивых конечно-разностных схем возмущениеаозрастаст.

Рассмотрим распределение погрешности на сетке в любой момент времени. Для удобства выберем момент времени =0. Схематически это распределение погрешности показано на

е(х,0)

Рис. 3.11. Начальное распределение погрешности.

рис. 3.11. Предположим, что погрешность e{x,t) можно представить в виде суммы ряда Фурье

(3.104)

причем период основной частоты (т = 1) равен 2L. Нас интересует решение в интервале длины L, поэтому волновые числа

krn = mnlL, m = 0, 1, 2, ..., М,

где М - число отрезков длины Дл:, помещающихся в отрезке длины L. Например, если интервал длины 2L разбит на отрезки пятью узлами, то Л1 = 2, а в сумму ряда входят лишь гармоники

= kJ2n = m/2L, f, = 1/2L, m = 1, /о = 0, m = 0, /2=1 , m = 2.

Напомним, что частота указывает на число волн, помещающихся в отрезке длины 2L. Если т = 0, то /о = О, а соответствующее слагаемое описывает стационарную составляющую решения.

Так как погрешность удовлетворяет линейному уравнению, то поведение каждой гармоники, входящей в (3.104), можно рассмотреть независимо. I ассмотрим член 8(х, t) = bm{t)e.

Будем искать решение уравнения в виде ze. При t = О (п = 0) оно имеет вид е". Пусть z = e\ тогда

n anM at (3.105)

причем km вещественно, но а может быть и комплексным. Подставляя (3.105) в (3.103), получим

где r = аА/(Ах)2. Разделив на Vh использовав соотношение

cosP = (P + -P)/2,

получим

«A/ = l + 2r(cosp-l), где р = kmAx, При помощи тригонометрического тождества

sin2(p/2) = (l-cos Р)/2 перепишем последнее соотношение в окончательном виде

ei = 1 ~ 4г sin2(P/2). (ЗЛ06)

Так как для каждой гармоники 8j+ = "ejf, то погрешность округления не будет возрастать на каждом шаге по маршевой координате (времени), если \е\ не превосходит единицы. Следовательно, разностная схема устойчива при

I 1 -4rsin2(p/2)Kl. (3.107)

Коэффициент 1 - 4/sin(p/2) (равный отношению назы-

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

При решении неравенства (3.107) надо рассмотреть два возможных случая:

1. Рели (1-4/sin2(p/2))>0, то 4/sin2(p/2) > 0.

" 2. Сели (1-4/sin2(p/2))<0, то 4г sin2(p/2)-1 1.

Первое неравенство выполняется для всех г > О, а второе - лишь при г 1/2. Последнее неравенство и является условием устойчивости рассматриваемой конечно-разностной схемы; оно накладывает ограничение на соотношение шагов по времени и пространственной координате. Теперь можно легко объяснить, почему в примере, приведенном в конце п. 3.3.3, получались физи-