На втором полушаге используется уравнение (8.13), представляемое в виде

-0.5s,rf.v-i + (1 + Sy)Ttk - 0.bSyTtk\x =

=o.5s/; + (1 - s,) r;, -f 0.5s/; (8. i5)

в течение второго полушага решение на временном слое п -f 1 неизвестно, однако известно решение на промежуточном временном слое, обозначенном звездочкой. Система уравнений, <:вязанных со всеми узлами, расположенными на одной сеточной линии вдоль направления у (фиксированное /) решается

относительно Тл k = 2, NY-1. Процесс повторяется

для каждой из упомянутых сеточных линий, т. е. для / = 2, ... NX - 1.

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

С СТУ l-2s,sin(9,/2) 1-2S, sin (9,/2) и -и и - 1 + sin2 (х12) \ + 2sy sin2 (9/2) *

Как показывает изучение формулы (8.16), при любом значении Sx, Sy, Qx и Qy получим G 1. Однако, рассматривая по отдельности G и G, найдем, что если полный шаг по времени является безусловно устойчивым, то каждый из полушагов устойчив лишь условно. Впрочем, интерес представляет лишь полный шаг по времени.

Составная схема (8.12) и (8.13) согласуется с уравнением (8.1) и обладает ошибкой аппроксимации порядка О (АР, Дл:, Ау). Второй порядок точности является следствием симметрии схемы в точной аналогии с тем, как второй порядок точности схемы Кранка - Николсона был следствием симметрии по отношению к временному слою (п + 1/2).

Однако для достижения суммарной ошибки аппроксимации порядка O(At) необходимо для промежуточного решения T*j вводить граничные значения так, чтобы они были совместимы с алгоритмами для внутренних точек (8.14) и (8.15). Например, если ставятся граничные условия Дирихле, то задание Тх, k = bk при х=1 приводит к алгоритму с ошибкой аппроксимации порядка О (ДО- Чтобы ошибка аппроксимации лолучила порядок 0(At)y необходимо задать Tlx,k помощью

соотношений (8.12) и (8.13), т. е. по формуле

Tnx. ife = 0.5 {Ы + ЬТ) - 0.25 AL - ЬЦ (8.17)

Аналогичная проблема возникает и при использовании схем с приближенной факторизацией; надлежащие формы задания граничных условий обсуждаются в п. 8,3.2 и § 8.4.

Таким образом, можно прийти к выводу, что схема НПН для двумерного случая обладает требуемыми свойствами: безусловно устойчива, имеет второй порядок точности и допускает экономичный метод решения. Схема НПН распространяется и на трехмерное уравнение диффузии, при решении которого вместо алгоритмов (8.14) и (8.15) необходимо вести счет на трех дробных шагах, каждый из которых занимает время Д 3. В случае трех измерений схема НПН оказывается экономичной, имеет второй порядок точности по пространству, но является лишь условно устойчивой. Для ее устойчивости необходимо, чтобы Sx. Sy, Sz 1.5, где Sz = azM/Az.

8.2.2, Обобщенная двухслойная схема

Попытаемся теперь дать обобщение идеи расщепления. Неявная двухслойная конечно-разностная схема общего вида применительно к уравнению (8.1) может быть записана в виде

- (1 - Р) (axLxxTl , + ayLyyTl О -- р [aUJt + aLyyTl = О, (8.18)

Ai /. Л = i /, Л - i /. /fe.

Здесь величину АГ/, л можно рассматривать как поправку к решению на временном слое п, требуемую для продвижения на временной слой (п+1). Для минимизации накапливающейся ошибки округления полезно, чтобы величина KTk фигурировала в явном виде в тексте компьютерной программы. Роль коэффициента р в уравнении (8.18) состоит в задании определенных весовых множителей при членах, соответствующих временным слоям ли (/z-f 1); та же идея была использована и при выводе уравнения (7.24).

С помощью перестановки членов уравнения (8.18) можно получить неявный алгоритм для определения АГ/, . Вначале член Tk разлагается в ряд Тейлора в окрестности п-то временного слоя в соответствии с формулой (3.16), а именно

= П. + МЩ + 0.5 + ....

Подстановка формулы (8.19) в уравнение (8.18) дает

(8.20)

или после некоторой перегруппировки

[1 - Р А/ (aUx + ayLyy)] ATk = А/ + %L/i/) П (8.21)

В левой части соотношения (8.21) фигурируют алгебраические операторы, действующие как в одном, так и в другом направлениях. Чтобы можно было воспользоваться алгоритмом Томаса, соотношение (8.21) заменяется другим соотношением, получаемым в результате приближенной факторизации, а именно

<1 - Р ШхЬхх) (1 - Р ШуЬуу) ATtk = (axLxx + ayLyy) Г?. (8.22)

Как показывает сравнение соотношений (8.21) и (8.22), в левой части последнего содержится дополнительный член

fMaxayLxxLyyATlX.

Отсюда следует, что соотношение (8.22) аппроксимирует (8.21) с точностью 0(Д2).

На каждом шаге по времени выполнение соотношения (8.22) может быть реализовано в форме двухэтапного алгоритма. Первый этап сводится к решению следующей системы уравнений, справедливых на каждой из сеточных линий, направленных параллельно оси х (линии постоянного k на рис. 8.2):

(1 - Р MaLxx) ДГд k = А/ {axLxx + ayLyy) Т], (8.23)

Это уравнение служит для определения АГ/, . что можно рассматривать как промежуточное приближение для АГ/,. Если Lxx представляет собой трехточечный оператор с центральной разностью, то уравнения (8.23) будут иметь форму трехдиагональной системы, которая может быть эффективно решена с помощью алгоритма Томаса (п. 6.2.2). На втором этапе решается следующая система уравнений:

(1 ~ Р MayLyy) ATtk = АГ;. ь (8.24)

что можно аппроксимировать формулой

-J-fO(A/2). (8.19)