В этом отношении некоторые трехслойные по времени схемы оказываются более эффективными, чем схема Кранка - Николсона.

7,2,3, Обобщенная трехслойная схема

Для уравнения диффузии можно составить обобщенную трехслойную схему, охватывающую, в частности, и уравнение (7.24). Эта схема имеет вид

Т - - а [(1 - Р) UJ1 + иГГ\ = О, (7.25)

где AT = - ТК Включение в схему дополнительного временного слоя влечет за собой повышенные требования к памяти, чтобы хранить решение. Однако современная тенденция сводится к тому, что компьютерная память становится дешевле и больше по объему. Есть и второй эффект, связанный с необходимостью выделения дополнительного времени на исполнение алгоритма, порядка 10-15%, нужного для обработки дополнительных членов.

Особенно эффективная трехслойная схема получается, если выбрать у = 0.5, р= 1.0. Эта схема имеет ошибку аппроксимации 0(Д/2, Ах), является безусловно устойчивой, может решаться с помощью алгоритма Томаса и способствует затуханию искусственно вызванных колебаний, обсуждавшихся ранее в связи с упоминанием о жестких задачах. Мы будем называть эту схему трехслойной чисто неявной (ТСЧН) и воспользуемся ею в дальнейшем при обсуждении вопроса о приближенной факторизации (см. § 8.2 и 8.3). Схема ТСЧН обладает полезным свойством Л-устойчивости (см. § 7.4).

Свойства некоторых из весьма многочисленных численных схем для моделирования уравнения диффузии (7.1) показаны в табл. 7.1. Намного больше схем дается в книге (Richtmyer, Morton, 1967].

7,2,4, Схемы повышенного порядка

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

- а [fiUJr + (1 - Р) LxJl] = О, (7.26)

W/--д]2-

Стоит отметить, что уравнение (7.26) по своей структуре подобно уравнению (7.13). При применении метода конечных элементов Mx = {j, Jy j]y так что

ATj = i. ДГу, +1 дг, +1 ATj,. (7.27)

В случае метода конечных разностей Мх = {О, 1, 0}. Параметры

Y и р могут быть подобраны так, чтобы обеспечить заданные уровни точности и/или устойчивости. В п. 7.2.2 различным схемам соответствовало значение y = О- Частный вариант выбора

Y = О, р = 0.5 приводит к методу Кранка - Николсона. Результаты как для конечно-элементной, так и для конечно-разностной форм метода Кранка - Николсона приводятся в табл. 7.3. Вариант выбора y = 0-5, Р = 1.0 кратко обсуждается в п. 7.2.3. В данном пункте параметр y будет рассматриваться как свободный, но параметр р будет считаться функцией y-

Разложение членов уравнения (7.26) в ряд Тейлора относительно узла (У, п) дает следующее выражение для главного члена ошибки усечения:

Е] 5 А;с (0.5 + Y + - р) . (7.28)

где оператор массы записывается в виде

М = {6, 1-26, 6). (7.29)

Это выражение включает в себя как конечно-разностную (б = = 0), так и конечно-элементную (6=1/6) формулировки. В выражении (7.28) для ошибки аппроксимации исключены все производные по времени, как это было сделано в табл. 7.1.

Рассмотренные ранее схемы (см., например, п. 7.2.3) соответствовали выбору р = 0.5 -f Y- Ясно, однако, что при выборе

p==0.5 + Y + -=7 (7.30)

должна быть достижима точность четвертого порядка. Это обстоятельство в свою очередь дает мотивировку выбора M;c = {-, -,-jV }, так как при таком Мх и при р = 0.5-f y получается ошибка аппроксимации четвертого порядка.

Уравнение (7.26) последовательно применяется к каждому узлу, что дат трехдиагональную систему уравнений вида

ЛуГ?:!-/ + ВуГГ + С/Г?/ = (1 + 2у) MJ } - уМхТ- +

+ (l~-p)sO?, (7.31)

Л, = Су = (1 + у)б-5р, 5, = (l + Y)(l-26) + 2sp, 0;-=Г?-1~2Г? + Г/%,.

Несмотря на то что число операций, требуемое для решения уравнения (7.31), значительно больше, чем для чисто неявной схемы или схемы Кранка - Николсона, данная схема оказывается намного более точной (см. п. 7.2.5).

7.2.5. DIFIM: численные результаты применения неявных схем

В данном пункте проводится сравнение точности схемы повышенного порядка (7.31) с точностью неявных схем более низкого порядка, а также явных схем низкого и высокого порядков, рассмотренных в § 7.1.

Программа DIFIM представляет собой обобщение программы DIFF, направленное на построение численного решения уравнения (7.1) при граничных, условиях (7.18) путем решения последовательности уравнений (7.31) так, чтобы решение продвигалось во времени для всех внутренних узлов. Для решения уравнения (7.31) используются подпрограммы BANFAC и BANSOL. Если учесть, что Л/, В/ и С/ не зависят от времени, то обращение к BANFAC необходимо сделать лишь один раз, на первом шаге по времени. Распечатка программы DIFIM приведена на рис. 7.7. Программа DIFIM может работать по пяти вариантам схем, приведенным в табл. 7.4 и соответствующим различным представлениям р как функции у, а также различным вариантам оператора Мх. Случай 1, когда у = 0, соответствует схеме Кранка - Николсона (7.26). Многие из основных параметров, используемых в программе DIFIM, совпадают с теми, которые применялись в программе DIFEX (см. табл. 7.2). Дополнительные параметры приводятся в табл. 7.5. Типовая выдача, характеризующая решение при ME = 3 и у = 0, показана на рис. 7.8.

Программа DIFIM использовалась для сравнения результатов применения различных методов, перечисленных в табл. 7.4. Сравнительные данные приводятся в табл. 7.6 и относятся к точности решения на различных сетках (Ах = 0.2, 0.1 и 0.05) при =12.0. Эти результаты были получены при а = 0.01 и