Разностная схема третьего порядка точности для решения жесткого обыкновенного дифференциального уравнения с линейными коэффициентами Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 11, № 2, 2006

РАЗНОСТНАЯ СХЕМА ТРЕТЬЕГО ПОРЯДКА

ТОЧНОСТИ ДЛЯ РЕШЕНИЯ ЖЕСТКОГО ОБЫКНОВЕННОГО ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ С ЛИНЕЙНЫМИ КОЭФФИЦИЕНТАМИ

В. Г. Зверев Томский государственный университет, Россия e-mail: zverev@niipmm.tsu.ru

A new implicit difference scheme for numerical analysis of the stiff Cauchy problem for a first order ordinary linear differential equation is proposed. This scheme is constructed using the increased accuracy Taylor expansion of function in the vicinity of the desired solution and the integration of the differential equation. In the case of linear coefficients the obtained difference scheme has the third order of approximation. Good rate of convergence to the exact solutions is shown on test examples in a wide range of a grid parameter. Comparison with other known one-step methods is carried out.

Введение

Задача Коши для обыкновенного дифференциального уравнения (ОДУ) первого порядка является математической моделью широкого круга прикладных проблем в различных областях науки. В качестве примера могут быть названы механика жидкости и газа, химическая кинетика, биология и т. д. Классы уравнений, для которых возможны точные решения в элементарных функциях, крайне узки и не охватывают всего многообразия возникающих на практике ситуаций. Поэтому наиболее общим подходом к получению решения ОДУ является применение численных методов [1, 2].

Несмотря на простоту математической формулировки, задачи с начальными данными обладают свойством жесткости, которое проявляется при наличии в решении быстро и медленно меняющихся компонент [3]. Математически это связано с появлением малого параметра при производной, что приводит к особенностям типа пограничного слоя. Последнее характерно для задач химической кинетики [4], где скорости реакций различаются на несколько порядков, а также для задач неравновесной газовой динамики [5]. В указанных условиях классические численные методы интегрирования сталкиваются со значительными трудностями, что приводит к неудовлетворительным результатам по точности и экономичности алгоритмов.

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

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

параметру сходимостью. Это прежде всего методы экспоненциальной подгонки, использование аналитических решений при построении разностных аналогов, сгущающихся по определенному закону сеток и др. Данной проблеме посвящены работы [6-13] отечественных и зарубежных авторов.

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

В качестве примера можно назвать известную в практике газодинамических расчетов неравновесных течений полидисперсных сред полностью неявную разностную схему второго порядка точности [14, 15]. Особенностью ее построения является использование отрезка ряда Тейлора для вычисления зависимой переменной в полуцелом узле через искомое решение с заменой первой производной соответствующим выражением из дифференциального уравнения. Использование этого приема позволило обеспечить полностью неявную форму разностного аналога, второй порядок аппроксимации и правильную асимптотику численного решения в случае малых значений параметра при производной, что очень важно для задач с пограничными слоями. Учитывая эти обстоятельства, представляют интерес дальнейшее повышение точности и расширение возможностей схем такого класса.

Цель настоящей работы — в развитие исследований [16, 17] разработка неявного экономичного А-устойчивого метода повышенного порядка точности для численного интегрирования жесткой задачи Коши для обыкновенного линейного дифференциального уравнения первого порядка.

1. Математическая постановка задачи и анализ предельных случаев

Рассмотрим задачу Коши для обыкновенного линейного дифференциального уравнения первого порядка следующего вида:

¿и

Ьп(х) = е-—+ а(х)и(х) = f (х), и(0) = и0, х Е (0,Х). (1)

ах

Здесь е Е (0,1] — малый параметр; а(х), f (х) — гладкие функции, причем а(х) > ао > 0 отделена от нуля для всех х > 0. Задача (1) имеет единственное решение [8, 9], согласно принципу максимума для него имеет место априорная оценка, не зависящая от е:

||и(х)||те < |по| + /ао, 0 < х < X, ||и(х)||те = тах |п(х)|. (2)

При больших и умеренных значениях е ~ 1 решение уравнения (1) не содержит каких-либо особенностей, определяется видом функций а(х), f (х) и может быть легко получено численно с помощью явных и неявных классических схем [8].

При стремлении параметра е ^ 0 решение (1) приобретает особенность типа пограничного слоя в силу того, что редуцированная часть задачи не удовлетворяет начальному

условию f (х)/а(х) = щ- В качестве иллюстрации рассмотрим задачу Коши

еи'(х) + и(х) = х, и(0) = 1, х € (0,1),

(3)

с точным решением, содержащим два слагаемых:

и(х) = (х — е) + (1 + е) ехр(—х/е).

(4)

Первое слагаемое в (4) представляет собой медленно меняющуюся функцию при любом значении е € (0,1]. Второе слагаемое определяется отношением р = х/е и быстро уменьшается в начале участка интегрирования. На рис. 1 показана зависимость и(х) для различных значений е. При уменьшении е образуется область резкого изменения функции — пограничный слой толщиной ~ е (р ~ 1), так как начальное условие и(0) = 1 не удовлетворяет решению редуцированной задачи

Классические схемы численного интегрирования в этом случае наталкиваются на серьезные трудности. Явные разностные методы, например схема Эйлера, в силу условия устойчивости при малых значениях е имеют жесткое ограничение вида р^ = к/е < 1 на шаг к интегрирования. Необходимость его выполнения распространяется и на участок медленного роста функции, протяженность которого не связана и несоизмерима с толщиной пограничного слоя. Данное обстоятельство делает невозможным применение явных схем к решению задач такого класса.

В неявных разностных схемах ограничение по устойчивости на шаг интегрирования формально отсутствует. Однако схемы первого порядка точности из-за большой ошибки аппроксимации сильно размазывают зоны пограничного слоя, что нежелательно в численных расчетах. Как отмечается в [8], для неявной схемы Эйлера ошибка сеточного решения мала при малых и больших значениях р^ = к/е< 1 и становится значительной,

и(х) = х, и(0) = 0.

(5)

0.0

0.0 0.2 0.4 0.6 0.8

1.0

X

Рис. 1. Аналитическое (1)—(3) и численное решения задачи Коши. Кривые 1-3 соответствуют е = 1; 0.1; 0.01; расчеты (е = 0.1, Н = 0.25): △ — метод Эйлера (21), * — схема (8) [14, 15], □ — схема (20), о — схема (19).

когда е и к одного порядка. На рис. 1 показано численное решение по этой схеме при е = 0.1 и постоянном грубом шаге к = 0.25 (р^ = 2.5) с явно выраженными сглаживающими свойствами. Относительная ошибка 8 = | (и — и^)/и^| на первом шаге составляет ~ 48 %. Лучше качественно и количественно при этих параметрах выглядит ситуация для неявной А-устойчивой схемы второго порядка [14, 15], где уровень ошибки падает до ~ 24%.

Таким образом, практические вопросы повышения точности численного решения наряду со схемами специального вида требуют разработки и применения экономичных неявных схем повышенного порядка аппроксимации.

2. Разностная схема второго порядка точности

Следуя источникам [14, 15], кратко рассмотрим основные этапы построения неявной схемы второго порядка точности. Разностный аналог дифференциального уравнения вытекает из применения к уравнению (1) формулы прямоугольников

иг+1 = и + кге-1 [/¿+1/2 — ат/2ит/2], (6)

где г — номер слоя интегрирования; кг = х^+1 — хг. Разлагая неизвестную функцию иг+1/2 в ряд Тейлора до членов О(к2) в окрестности искомого решения и^+1 и заменяя первую производную соответствующим выражением из дифференциального уравнения (1), получим

ит/2 = иг+1 — кг <+1/2 + 0(к2) = ит(1 + кге-1ат/2) — кге-1/т/2 + О(кг2). (7) Подставляя (7) в (6) и преобразуя, получим итоговый вид разностного аналога

иг + кге-1[/г+1/2 + /г+1^г+1/2/2] , -1

иг+1 =--.-77;-, гг+1/2 = аг+1/2кге , (8)

1 + ^г+1/2 + ^г+1/2^г+1/2

где г — сеточный параметр. При табличном (численном) способе задания функций а(х), /(х) их значения в полуцелых узлах определяются со вторым порядком точности:

аг+1/2 = (аг + аг+1)/2, /г+1/2 = (/г + /г+1)/2. (9)

Из (8) следует, что при е ^ 0 иг+1 ^ (//а)г+1, при е ^ то иг+1 ^ иг. Таким образом, в предельных случаях разностный аналог (8) аппроксимирует предельные решения уравнения (1). В случае кусочно-постоянных /(£) = /, а(£) = а схема (8) упрощается и принимает вид

иг + кге-1/[1 + г/2] -1

иг+1 = [1 + г + г2/2] , г = акге . (10)

Точное решение (1) в этом случае запишется как

иг+1 = иге(г ) + (кге-1)/в(г), в(г) = в-*, в (г) = (1 — в-* )/г, (11)

где в(г), в (г) — сеточные функции [16, 17]. В [12] оно лежит в основе разностной схемы, теоретическая оценка равномерной по е сходимости которой составляет О (к). Графики в (г) и однородного решения в (г) уравнения (1) показаны на рис. 2 и 3. Сравнивая (10)

Рис. 2. Сеточная функция ехр(—г) и ее ап- Рис. 3. Сеточная функция в (г) и ее аппрок-проксимации: 1 и 1', 2, 3 — первый, второй и симации: в1 (г), в2(г), вз(г) — первый, второй третий порядки точности. и третий порядки точности.

и (11), видим, что (10) вытекает из точного решения (11) при следующем разложении экспоненты:

ехр(г) « 1 + г + г2/2 + 0(г3).

Близость используемых в (10) аппроксимаций сеточных функций к своим оригиналам (рис. 2 и 3) обусловливает хорошую точность схемы. Фактически (10) представляет собой рациональную аппроксимацию второго порядка точного решения (11). Таким образом, дальнейшее повышение возможностей расчетной схемы связано с улучшением точности численного решения однородного уравнения.

3. Методика построения разностной схемы повышенного порядка аппроксимации

Как следует из (7), (8), линейный член разложения иг+1 в ряд Тейлора в окрестности (г+1)-го слоя обеспечивает второй порядок аппроксимации схемы (8). Однако формальный учет последующего члена разложения в (7) не дает желаемого результата в улучшении точности однородного уравнения и схемы в целом. Причиной является погрешность формулы прямоугольников. Поэтому необходимо применение квадратурной формулы более высокого порядка или использование другого подхода. Последний случай рассматривается ниже. Проинтегрировав уравнение (1) от хг до получим

«г+1 = «г + £

-1

Хг+1

Ег+1

/т - а(0и(£ж

(12)

X

X

г

Представим «(£) при £ € [£¿,£¿+1 ] в виде отрезка ряда Тейлора в окрестности (г + 1)-го слоя с погрешностью 0(Л3):

«(£) = «¿+1 + «¿+1(£ - £¿+1) + «¿'+1 (£ - £¿+1)72 + 0(£ - £г+1)3. (13)

Заменим производные соответствующими выражениями из уравнения (1), обозначив

^¿+1 = «¿+1 = е-1(/ - аи^+ь (14)

N¿+1 = «¿'+1 = е-1[/' - ае-1 / + и(а2е-1 - а')^+ь

Подставляя (13) в (12), с учетом обозначений (14) запишем

Ег+1 х1+1

«¿+1 = « + е-1 / /(£)а£ - «¿+1 / а(£)а£ -

, ж

Жг+1 жг+1

- M¿+^ J а(£)(£ - £¿+1 )а£ - ^ / а(£)(£ - £¿+1)^ | + О(^). (15)

хг хг

Необходимо заметить, что выражение (15) получено без каких-либо допущений и ограничений на вид функций а(£), /(£) и имеет погрешность 0(Л4), связанную с усечением в (13) ряда Тейлора. Поэтому оно может являться основой для получения различных разностных схем.

Конкретизируем вид зависимостей коэффициентов уравнения а(£), /(£). Для одно-шагового метода и табличной формы использования данных имеет смысл рассматривать линейную зависимость

а(£) = a¿+l + а^+1 (£ - £¿+1), a¿+l = (a¿+l - a¿£ € [£¿,£¿+1], (16) /(£) = /¿+1 + Л+1(£ - £¿+l), /¿+1 = (/"¿+1-

Используя (16), точно вычислим интегралы в (15). В результате получим

«¿+1 = « + Ле-1 (/¿+1/2 - a¿+l/2M¿+l + M¿+lah¿/2 - N¿+1^/6) + 0(Л4), (17)

где a¿+l/2, /¿+1/2 соответствуют (9), другим обозначениям отвечают формулы

а = ^¿+1 + 2a¿)/3, а = (a¿+l + 3a¿)/4.

Подставляя выражения для производных М^+1 и N¿+1 (14) и приводя подобные, получим

= « + ^¿е-1 [/¿+1/2 + /¿+1 (¿/2) + (а/ - Л/'Ы(а/6)] (1§)

«¿+1 1 + а+1/2 + ¿¿+1^/2 +(г2 - а'е-1Л|^+1(а/6) , ( )

где г = аЛе-1, а = аЛе-1.

Разностная схема (18) при линейных коэффициентах а(£), /(£) имеет порядок аппроксимации 0(Л3). Преобразуем (18), используя выражения (16) для а', /', найдем

= « + ^¿е-1 [/¿+1 (1 + 25/3 + а+1а/3)/2 + /¿(1 + а/3)/2] (19)

«¿+1 1 + ¿¿+1/2 + (¿¿+125/3 + ¿¿а/3)/2 + а2+1а/6 , ( )

X

где г = ак^г-1, а = (3а^+1 + 5а^)/8. Из принципа максимума и положительности коэффициента а(х) (сеточного параметра г) следует А-устойчивость схемы для любых значений к, а(х).

В предельном случае при г — 0 из (19) следует ^¿+1 — (//а)^+1, что отвечает правильной асимптотике.

Если в (13) ограничиться разложением второго порядка, то данный способ построения приводит к схеме вида

иг + кг^-1 [/г+1/2 + /г+1 (г/2)] ^¿+1 =----. (20)

¿+1 1 + £¿+1/2 + ¿¿+1^/2 У У

Сравнивая (20) с (8), видим, что при одинаковой структуре записи схема (20) имеет отличие в способе усреднения коэффициента а при старшей степени к2. Любопытно, что разностный аналог (8) также может быть получен на основе выражения (15). Для этого в интеграле при М^+1 функцию а(£) на [х^,х^+1 ] следует заменить приближенным значением

аг+1/2 •

Наконец, если в (15) ограничиться двумя первыми интегралами, для вычисления которых использовать формулу правых прямоугольников, то получим неявную схему Эйлера первого порядка

иг+1 = (и + кгг-1/г+1)/[1 + ¿¿+1]. (21)

4. Анализ упрощенных вариантов схемы

Рассмотрим случай кусочно-постоянного коэффициента а(£) = а и линейного источника /(£). Точное решение задачи Коши в этом случае можно представить в виде

иг+1 = иге(г) + кгг 1[/г+1^(г) + /¿п(г)],

(22)

Рис. 4. Сеточные функции и их аппроксимации: £1^), — первый, £2^), —

второй, пэ(^) — третий порядки точности.

е(а) = exp(—а), £(а) = [а - 1 + ехр(-а)]/а2 п(а) = [1 - ехр(-а)(а + 1)]/а2

-1

£(а) + п(а) = в (а) = [1 - ехр(-а)]/а, а = аЛ^е

где е(а), £(а), п(а), в (а) — сеточные функции [16, 17]. Их графики показаны на рис. 2-4. Функция е(а) является решением однородного уравнения (1), £(а), п(а) играют роль весовых множителей при источнике на (г + 1)-м и г-м слоях.

Нетрудно видеть, что схема (19) соответствует записи (22) при следующем представлении сеточных функций:

ез(а)

2

аа

-1

£з (а) = 0.5 ез(а).

2 а2

1 + 3а + У

ез(а^

1+

Пз(а) = 0.5

Отметим, что для схем второго порядка (8), (20) эти выражения имеют вид б2(а)= [1 + а + а2/2]-1, £2(а) = 0.5[1 + а] б2(а), П2(а) = 0.5в2(а).

(23)

(24)

На рис. 2-4 хорошо видно, что приближение (23) значительно точнее, чем (24), воспроизводит сеточные функции точного решения (22) во всем диапазоне значений аргумента, особенно при малых и больших а. Асимптотический анализ (22) показывает, что аппроксимации (23) следуют из точного решения (22) при использовании разложения экспоненты в ряд Тейлора до третьего порядка точности.

Для кусочно-постоянного источника точное решение (22) переходит в (11) и разностная схема (19) принимает вид

«¿+1

« + ^¿е-1/[1 + а/2 + а2/6]

1 + а + а2/2 + а з/6

что соответствует следующей аппроксимации сеточных функций е(а), в (а): ез(а) = [1 + а + а 2/2 + аз/6]-1, вз(а) = [1 + а/2 + а2/6]вз (а).

(25)

(26)

5. Результаты численных расчетов и их анализ

Результаты численного решения рассмотренной задачи 3 с помощью предложенных схем (20) и (19) при е = 0.1 и грубом шаге Л = 0.25 показаны на рис. 1. Для сравнения приведены данные по неявной схеме Эйлера (21) первого порядка. С теоретической точки зрения схемы (8) и (20) дают одинаковые результаты, качественно и количественно превышающие возможности неявной схемы Эйлера. Однако они уступают результатам, полученным по схеме (19) третьего порядка аппроксимации. Ее применение позволило сократить относительную ошибку численного решения на первом грубом расчетном шаге до ~ 12 % и с приемлемой для практики точностью воспроизвести особенности формирования (кривая 2) пограничного слоя. В числе других преимуществ схемы (19) — высокая экономичность, так как используются только рациональные выражения. Необходимо заметить, что более трудоемкая разностная схема специального вида, предложенная в [16, 17], решает эту задачу точно для любых значений Л и е.

1.0

и(х)

0.8-

0.6-

0.4-

0.2-

0.0

0.0

0.5

1.0

1.5

2.0

X

Рис. 5. Решение задачи (27) при различных е: кривые 1-3 соответствуют е = 1; 0.1; 0.01; расчеты: Д — метод Эйлера [12], V — схема [12], * — схема (8) [14, 15], □ — схема (20), о — схема (19); е = 1, Н = 1; е = 0.1, Н = 0.25.

Для более детальной оценки возможностей предложенных разностных схем (20) и (19) рассмотрим тестовую задачу Коши с переменным коэффициентом а(х) [12]:

Аналитическое решение п(х) при различных значениях е показано на рис. 5. Для сравнения здесь в соответствии с работой [12] приведены результаты по неявной схеме Эйлера и равномерной по е схеме (11) первого порядка. Видно, что полученные при грубом шаге по схеме (19) данные графически близки к кривым 1, 2 точного решения. Более подробно ошибка Ди численного решения пи при различных значениях к и е указана в таблице:

Ее анализ показывает, что, несмотря на иной способ усреднения коэффициента а(х) в (20), ошибки по схемам (8) и (20) второго порядка незначительно отличаются друг от друга. Схема (19) третьего порядка по сравнению со схемами (8) и (20) имеет гораздо более низкий уровень ошибки во всем рассматриваемом диапазоне значений к и е.

Заключение

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

В основе построения схемы лежат разложение функции в ряд Тейлора в окрестности искомого решения и прямое интегрирование дифференциального уравнения.

ей'(х) + (1 + х)п(х) = (1 + х), п(0) = 0, х € [0, 2], п(х) = 1 - ехр(—(2х + х2)/(2е)).

(27)

Ди = тах |п(х) — пи(х)|, х € [0, 2].

Погрешность численного решения задачи (27) при различных значениях Н и е

Схема (8) [14, 15] Схема (20) Схема (19)

h £ = 1 £ = 0.1 £ = 0.01 £ = 1 £ = 0.1 £ = 0.01 £ = 1 £ = 0.1 £ = 0.01

1 2.7e-2 6.0e-3 6.6e-5 3.8e-2 6.7e-3 7.4e-5 4.1e-3 1.0e-3 1.2e-6

0.1 6.2e-4 3.1e-2 1.4e-2 8.1e-4 3.2e-2 1.5e-2 2.0e-5 6.2e-3 3.6e-3

0.01 6.8e-6 5.4e-4 3.2e-2 8.9e-6 5.7e-4 3.2e-2 2.3e-8 1.2e-5 7.0e-3

0.001 6.9e-8 5.8e-6 5.7e-4 9.0e-8 6.1e-6 5.7e-4 2.4e-11 1.3e-8 1.4e-5

0.0001 6.9e-10 5.9e-8 6.1e-6 9.0e-10 6.2e-8 6.1e-6 2.5e-14 1.3e-11 1.5e-8

Тестовыми расчетами подтверждена хорошая для практики сходимость численных результатов к точным решениям, в том числе и при малых значениях параметра при производной.

Список литературы

[1] Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989.

[2] КАлиткин Н.Н. Численные методы. М.: Наука, 1978.

[3] РАкитский Ю.В., Устинов С.М., Черноруцкий И.Т. Численные методы решения жестких систем. М.: Наука, 1979.

[4] Евсеев Г.А., КАлюжный В.В. Экономичный метод численного интегрирования уравнений химической кинетики // Численные методы механики сплошной среды: Сб. науч. тр. АН СССР. ВЦ; ИТПМ. 1974. Т. 5, № 3. С. 21-28.

[5] Пирумов У.Г., Росляков Г.С. Газовая динамика сопел. М.: Наука, 1990.

[6] Васильева А.Б., Бутузов В.Ф. Асимптотические разложения сингулярно возмущенных уравнений. М.: Наука, 1973.

[7] Ломов С.А. Введение в общую теорию сингулярных возмущений. М.: Наука, 1981.

[8] Дулан Э., Миллер Дж., Шилдерс У. Равномерные численные методы решения задач с пограничным слоем. М.: Мир, 1983.

[9] Багаев Б.М., ШАйдуров В.В. Сеточные методы решения задач с пограничным слоем. Ч. 1. Новосибирск: Наука, 1998.

[10] Hairer E., Wanner G. Solving Ordinary Differential Equations II: Stiff and Differential-algebraic Problems. Berlin: Springer-Verlag, 1991.

[11] Roos H.G., Stynes M., Tobiska L. Numerical Methods for Singularly Perturbed Differential Equations. Berlin: Springer-Verlag, 1996.

[12] Титов В.А., Шишкин Г.И. Численное решение задачи Коши для обыкновенного дифференциального уравнения с малым параметром при производной // Численные методы механики сплошной среды: Сб. науч. тр. АН СССР. ВЦ; ИТПМ. 1978. Т. 9, № 7. С. 112-121.

[13] Боглаев И.П. О численном интегрировании сингулярно возмущенной задачи Коши для обыкновенного дифференциального уравнения // Журн. вычисл. математики и мат. физики. 1985. Т. 25, № 7. С. 1009-1022.

[14] ВАСЕНИН И.М., Архипов В.А., Бутов В.Г. и др. Газовая динамика двухфазных течений в соплах. Томск: Изд-во Том. ун-та, 1986.

[15] Рычков А.Д. Математическое моделирование газодинамических процессов в каналах и соплах. Новосибирск: Наука, 1988.

[16] Зверев В.Г. О численном решении сингулярно возмущенной задачи Коши для дифференциального уравнения первого порядка // Вычисл. гидродинамика. Томск: Изд-во Том. ун-та, 1999. С. 73-82.

[17] Зверев В.Г., Гольдин В.Д. Об одной разностной схеме для решения обыкновенного дифференциального уравнения первого порядка с малым параметром при производной // Сопряженные задачи механики и экологии. Томск: Изд-во Том. ун-та, 2000. С. 166-174.

Поступила в редакцию 10 августа 2005 г., в переработанном виде —13 января 2006 г.