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

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

Том 13, № 3, 2008

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

уравнения

В. Г. Зверев, В. Д. Гольдин Томский государственный университет, Россия e-mail: [email protected]

A special one-step finite difference scheme for numerical analysis of the stiff Cauchy problem for a first order ordinary differential equation is proposed. It has uniform convergence on small parameter and is based on application of the asymptotic Laplace method for estimation of the integral solution of the Cauchy problem. A second order rational approximation of this special difference scheme is given. The comparison of the proposed difference scheme as well as its approximation efficiency with other known one-step difference methods is carried out on test problems.

Введение

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

Наличие быстро и медленно меняющихся компонент решения, а также несоизмеримость масштаба процесса и общего промежутка интегрирования делают задачу жесткой и обусловливают трудности численного интегрирования [5]. Жесткие задачи возникают почти во всех прикладных областях и в наибольшей степени характерны для задач химической кинетики, горения, неравновесной газовой динамики и т.д. [6, 7].

Математически жесткость связана с появлением в уравнении малого параметра при производной, стремление которого к нулю приводит к сингулярно возмущенной задаче Коши и обусловливает появление узкой области значительного изменения искомой функции — пограничного слоя (ПС) [8, 9].

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

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

малости погрешности аппроксимации. В конечном счете, это приводит к росту ошибки численного решения. В [10] рассматривались причины, сводящие на нет усилия классических схем при решении сингулярно возмущенной задачи Коши.

Практический интерес и необходимость создания алгоритмов, с приемлемой точностью работающих во всем диапазоне изменения коэффициентов и в области определения ОДУ, дали толчок к развитию специальных вычислительных технологий, обладающих равномерной по малому параметру сходимостью. Данной проблеме посвящена обширная библиография, например, [10-17] и цитируемая там литература.

Подходы, на основе которых строятся специальные схемы для решения жесткого ОДУ, весьма разнообразны. В [10] для этой цели используется метод экспоненциальной подгонки. В основе построения схемы [18] точности О(к) лежит аналитическое решение однородного ОДУ с постоянными коэффициентами. В [13] используются решения задачи с кусочно-постоянными коэффициентами, специальные формулы интегрирования, сгущающиеся по определенному закону сетки [17]. В [19], основываясь на методе построения точных разностных схем [20], получены специальные схемы первого и второго порядка для решения линейной и нелинейной сингулярно возмущенной задачи Коши. Следуя [19], в [21] построена равномерно сходящаяся схема, в которой используется аналитическое решение неоднородного уравнения с линейной правой частью. В [22] для этой цели применяются асимптотические методы вычисления интегралов.

При решении жестких задач практическое распространение получили также неявные абсолютно устойчивые одношаговые разностные схемы повышенного порядка точности. Они не содержат специальных функций, экономичны, обеспечивают асимптотику решения и приемлемую точность расчета задач с пограничными слоями. В качестве примера можно назвать схему второго порядка [23-25]. В [26] в развитие этого направления предложена схема третьего порядка для решения задачи Коши с линейными коэффициентами, а в [27] дано ее обобщение до пятого порядка включительно.

Цель настоящей работы — дальнейшее развитие исследований [21, 22] по разработке и апробации специальной схемы для численного интегрирования сингулярно возмущенной задачи Коши. Для ее построения используется интегральная форма решения и применяется асимптотический метод Лапласа для приближенной оценки интегралов.

1. Математическая постановка задачи и ее анализ

Рассмотрим линейную задачу Коши для жесткого ОДУ первого порядка с малым параметром е при производной:

¿и

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

ах

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

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

0<х<Х

При е ^ 0 решение уравнения (1) имеет особенность типа ПС, так как редуцированная часть задачи не удовлетворяет начальному условию, f(0)/а(0) = и0.

Для иллюстрации рассмотрим следующий пример [13] с известным точным решением, которое при различных е показано на рис. 1:

еи'(х)+ и(х) = 1 - х, и(0) = 0, х е (0,1), (3)

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

Из рисунка видно, что по мере уменьшения е на левом краю области интегрирования формируется ПС толщиной е, обусловленный наличием в (4) быстро изменяющейся компоненты решения (второе слагаемое), в то время как на остальной части решения доминирует медленная компонента (первое слагаемое).

Классические схемы численного интегрирования в этом случае наталкиваются на серьезные трудности. В [10] показано, что ошибки сеточного решения по явной и неявной схемам Эйлера уже на первом шаге зависят от отношения Н/е . Это означает отсутствие равномерной по параметру е сходимости. Жесткое ограничение по устойчивости Н/е < 1 на шаг интегрирования делает непригодным применение явной схемы при е ^ 0 в практических расчетах, так как необходимость его выполнения распространяется и на участок медленного хода решения. Неявная схема первого порядка сильно размазывает ПС, что нежелательно в численных расчетах. Согласно [10], ошибка решения велика, когда Н и е одного порядка.

Результаты расчета задачи при грубом шаге интегрирования показаны на рис. 1 (при е = 0.01, Н = 0.05, чтобы не загромождать рис. 1, б, они приведены в области пограничного слоя, х < 0.2). Явная схема Эйлера при таких параметрах из-за условия устойчивости не позволяет получить решение, качественно соответствующее точному. Неявная схема из-за схемной вязкости размазывает погранслойные особенности решения. В противоположность этому гораздо точнее выглядят результаты по неявной схеме второго порядка [23-25].

а б

Рис. 1. Аналитическое и численное решение задачи (3): кривые 1-3 — точное решение при £ = 1 (а); 0.1; 0.01 (б). Расчеты: Д — неявная схема Эйлера; * —схема второго порядка [23-25]; о, □ — специальная схема (11) и ее аппроксимация (13) второго порядка; £ = 1, Н = 0.5; £ = 0.1, Н = 0.25; £ = 0.01, Н = 0.05 (при х < 0.2)

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

2. Методика построения специальной разностной схемы

На отрезке [0,Х] определим неравномерную сетку так, чтобы точки разрыва коэффициентов а(ж), f (ж) служили узлами сетки:

{'i } ,

0,1,

, N; xo — 0, Xn — X; hi — Xi+1 — Xi.

(5)

Точное решение задачи Коши (1) на сеточном отрезке [ж^ж^] имеет следующее интегральное представление:

u(x) — u(xi) exp [—е ^(x)] + £

A(x) —

-i

exp

-£

-i

a(n)dn

f(£)d£,

(6)

a(e ж.

Конкретизация зависимости и (ж) происходит при подстановке выражений для а(ж), f (ж). Однако учет а(ж) уже в линейном приближении приводит к тому, что (6) не выражается в элементарных функциях. Поэтому необходимо получение расчетного выражения.

Для оценки интеграла в (6) воспользуемся асимптотическим методом Лапласа [28-30]. При малых £ ^ 0 основной вклад в интеграл вносит окрестность точки максимума показателя экспоненты. Так как а (ж) > 0, то он достигается на правом конце отрезка [ж^ж]. Применим к (6) интегрирование по частям, получим

u(x) — u(xi) exp [—e-1A(x)] + f(x) — f(x)exp [—e-1A(x)] —

17 К)!

I_a(£ )J exp -е

a(x) -i

a(xj) a(n)dn

d£-

(7)

Видно, что в (7) присутствует асимптотическое решение жесткой задачи (1) при £ ^ 0, а именно а(ж)и(ж) = f(ж). Оставшийся интеграл вносит малую поправку в решение. Возьмем в нем точно ту часть, которая связана с экспонентой, а комплекс ^ (£)/а(£)] /а(£) с погрешностью О (ж — ж^) заменим значением в средней точке ж = 0.5 (жi + ж). В результате получим

ш

a(£)

exp

—£

-1

a(n)dn

d£

f (x) f (xi)

a(x) a(xi)

= + O(x — xi) > е (1 — exp (—е 1A(

(^c 'X i) a

где a — 0.5 (a(xi) + a(x)). Для сомножителя имеет место оценка

е (1 — exp (—e-1A(x))) < A(x) < C1(x — xi), C1 — const

x

(8)

x

x

i

x

г

x

i

г

x

x

i

1

Поэтому отбрасывание второго слагаемого в (8) вносит погрешность О(х — х^)2. В итоге приближенную зависимость и(х) на сеточном интервале можно записать в виде

и(х) = и(хЛ ехр [—е-1А(х)1 + ^^ — ^^ехр [—е-1А(х)! —

а(х) а(х^)

/(х) /(хгП 1 — ехР (—е-1А(х)) (9)

— < - —

а(х) а(х^) ] ае-1(х — х^)

Положим в (9) х = х^+1 и вычислим А(х^+1) с погрешностью О (Н3) по формуле трапеций, которая точна для линейной функции:

Г Хг+1

А(хг+1) = а(п)(^)п ~ аН. (10)

./х;

Отметим, что при х = х^+1 отброшенное в (8) слагаемое имеет порядок О (Н3). В результате получим одношаговую разностную схему:

U+i = Uj exp(—z) + /j— — exp(—z) —

/¿+1

в(z), (11)

z = as hj, в (z)

которую также можно записать в виде

^ \ — exP(—z)

Uj+i = Ui exp(—z) + ( /j+i [1 — в(z)] + /j [в(z) — exp(—z)] 1 .

I a+i a J

Здесь z — сеточный параметр, в(z) — сеточная функция, введенная в [21, 22, 26].

z

3. Свойства специальной схемы и ее рациональная аппроксимация

Рассмотрим свойства схемы в случае предельных значений сеточного параметра г. Так как метод Лапласа является асимптотическим при е ^ 0, то уменьшение е не приводит к каким-либо трудностям. При больших г ^ то(е ^ 0) функция в(г) ~ 1/г и разностное выражение (11) сразу выходит на асимптотику и+1 ^ (/г+1/а^+1). Это означает, что для жесткой неоднородной задачи при а(х) > 0 преобладающим является весовой множитель у источника на (г + 1)-м слое. При малых г ^ 0 функция в (г) ~ 1 — г/2 + О(г2) и схема (11) принимает вид аналога формулы трапеций:

z

Uj+i = Uj + 2

/j+i + —j

При a(x) = const, /(x) — линейной функции, а также в другом важном практическом случае: /(x)/a(x) = const, a(x) — линейная функция, схема (11) дает точное решение задачи. Для произвольных a(x),/(x) она имеет второй порядок и в случае жесткой задачи (z ^ то, s ^ 0) обладает равномерной сходимостью. Это следует из того факта, что все остаточные члены в аппроксимационных выражениях (8), (10) от s не зависят. Технология построения разностной схемы (11) позволяет получить приближенное решение не только в узлах расчетной сетки (5), но и, опираясь на выражение

(9), на сеточном интервале (0,к^, а значит, и на всем промежутке интегрирования [0, X]. Следует отметить, что схема (11) обобщается на нелинейный случай, когда коэффициенты зависят от искомой функции. При этом (11) становится нелинейным уравнением относительно ^+1, которое решается методом итераций.

Из (11) путем преобразований можно получить разностный аналог, который приведен в [21] и построен на основе другого подхода. Для этого введем сеточные функции е(г), £(¿0, п(г), играющие роль однородного решения и весовых множителей при источнике на (г + 1)-м и г-м слоях, в результате получим следующий вид схемы [21]:

и+1 = и^е(г) + £(г) + ^фА , (12)

I ат а J

1 — в (г) ехр(—г) + г — 1 е(г) = ехр(—£ (г) =-=---,

•у 2

n(z)

z z2

в (z) — exp(—z) 1 — exp(—z)(z + 1)

z2

Сеточные функции е(г), в(г), £(г), п(г) и их аппроксимации подробно рассматривались в [21, 22, 26].

Представляет практический интерес рациональная аппроксимация специальной схемы (12). С этой целью выделим в £(г), п(г) функцию ехр(—г) и используем разложение второго порядка, получим следующие аппроксимации сеточных функций и вид схемы:

62(z)— [1 + z + z2/2 , ^2 (z) —0, 5(1 + z)e2(z), П2 (z) —0, 5(1 + z)e2(z),

Ui + z {(fi/ai) + (fi+1/ai+1)(1 + z)} /2 Ui+1 —-1 + z + z2/2-. (13)

В частном случае при a(x) — const (13) совпадает с хорошо известной схемой второго порядка [23-25], а также со схемой работы [26].

z

4. Вычислительный эксперимент и его анализ

Результаты численного решения рассмотренной выше задачи Коши (3) с помощью специальной схемы (11) и ее аппроксимации (13) второго порядка показаны на рис. 1. Как и следует теоретически, специальная схема решает эту задачу точно. Ее аппроксимация (13) дает одинаковые со схемой [23-25] результаты, которые качественно и количественно превышают возможности неявной схемы Эйлера первого порядка. Рассмотрим тестовую задачу с переменным коэффициентом а(ж) [18]:

£и'(ж) + (1 + ж)и(ж) = (1 + ж), и(0) = 0, ж € [0, 2], (14)

и(ж) = 1 — ехр (—(2ж + ж2)/(2£)) . (15)

Ее аналитическое (15) и численное решение и (ж) при различных значениях £ показано на рис. 2 (при £ = 0.01, к = 0.05 расчеты приведены в области ж < 0.2). Хорошо видно, что малые значения £ приводят к образованию пограничного слоя. Для сравнения, согласно работе [18], здесь показаны результаты расчетов по неявной схеме

Эйлера и равномерной по s схеме первого порядка (коэффициенты a(x),/(x) в обеих схемах берутся с предыдущего i-го слоя [18]). Максимальная погрешность Д^ = max |u(x) — Uh(x)| , x G [0, 2] численного решения Uh задачи при различных значениях шага h и параметра s приведена в табл. 1.

Видно, что специальная схема точно (до машинного представления числа) воспроизводит решение независимо от h и s, что подтверждает теоретические выкладки. Погрешности схемы (13) и схемы [23-25] незначительно отличаются друг от друга. Изменение шага сетки демонстрирует их второй порядок. Как видно из рис. 2, эти схемы с приемлемой для практики точностью воспроизводят все особенности формирования (кривые 2, 3) пограничного слоя. В числе других преимуществ аппроксимации специальной схемы (13) — экономичность, так как используются только рациональные выражения.

а б

Рис. 2. Аналитическое и численное решение задачи (14) при различных значениях е: кривые 1-3 — точное решение при е = 1 (а); 0.1; 0.01 (б). Расчеты: Д — метод Эйлера [18]; V — схема [18]; * — схема второго порядка [23-25]; о, □ — специальная схема (11) и ее аппроксимация (13) второго порядка; е = 1, Л = 1; е = 0.1, Л = 0.25; е = 0.01, Л = 0.05 (при х < 0.2)

Рис. 3. Аналитическое и численное решение * — схема [23-25]; о, □ — специальная схема N = 5; У> = 100

задачи (16): расчеты: Д — схема Эйлера [25]; (11) и ее аппроксимация (13) второго порядка;

Таблица 1. Погрешность численного решения задачи (14) при различных значениях Н и е

Схема Эйлера [18] Схема [18] Схема [23-25 , второй порядок

Н е = 1 е = 0.1 е = 0.01 е = 0.001 е = 1 е = 0.1 е = 0.01 е = 0.001 е = 1 е = 0.1 е = 0.01 е = 0.001

1 2.8е-1 9.1е-2 1.0е-2 1.0е-3 1.4е-1 4.5е-5 0.0 0.0 2.7е-2 6.0е-3 6.6е-5 6.7е-7

0.1 3.7е-2 1.5е-1 9.1е-2 9.9е-3 2.0е-2 1.8е-2 1.8е-5 0.0 6.2е-4 3.1е-2 1.4е-2 1.7е-4

0.01 3.9е-3 2.0е-2 1.3е-1 9.1е-2 2.1е-3 2.7е-3 1.8е-3 2.2е-6 6.8е-6 5.4е-4 3.2е-2 1.6е-2

0.001 3.9е-4 2.1е-3 1.8е-2 1.3е-1 2.1е-4 2.8е-4 2.8е-4 1.8е-4 6.9е-8 5.8е-6 5.7е-4 3.2е-2

Окончание табл. 1

Специальная схема (11) Схема (13), второй порядок

Н е = 1 е = 0.1 е = 0.01 е = 0.001 е = 1 е = 0.1 е = 0.01 е = 0.001

1 0.0 0.0 0.0 0.0 5.3е-2 7.8е-3 8.8е-5 8.9е-7

0.1 1.1е-16 1.1е-16 0.0 0.0 1.2е-3 3.4е-2 1.5е-2 1.8е-4

0.01 4.4е-16 2.2е-16 1.1е-16 0.0 1.4е-5 6.3е-4 3.2е-2 1.6е-2

0.001 1.8е-15 1.9е-15 4.4е-16 1.1е-16 1.4е-7 6.7е-6 5.7е-4 3.2е-2

Таблица 2. Результаты численного решения задачи (16)

N Сд(*) V (*) иЛ(*), N = 5 иЛ(*), N = 10

узла точное Схема Схема Схема Схема Схема Схема Схема Схема

Эйлера [25] [23-25] (11) (13) Эйлера [25] [23-25] (11) (13)

0 0 0 0 0 0 0 0 0 0 0 0

1 0.352 8.927 0.229 0 0.238 0.224 0.228 0.217 0.230 0.226 0.228

2 0.726 25.838 3.016 2.393 3.037 2.872 2.928 2.898 3.017 2.976 2.991

3 1.376 30.689 10.487 7.652 10.391 10.052 10.065 9.145 10.467 10.382 10.377

4 3.077 15.079 16.078 9.893 15.458 16.055 15.309 12.687 15.934 16.106 15.863

5 1.6е16 9.3е-31 7.0е-15 5.030 2.9е-15 5.9е-15 5.9е-15 5.094 1.9е-15 3.8е-15 3.8е-15

Рассмотрим тестовый пример из [25] с более сложной зависимостью коэффициентов: du

— + CR(t)(u - V(t)) = 0, u(0) = 0, t G [0,n/2], (16)

CÄ(t)=tg(t), V(t) = Vot2 cos2(t), u(t) = cos(t) { Vo [2t sint - (t2 - 2) cos t] + uo - 2Vo} . (17)

Аналитическое (17) и численное решение задачи показано на рис. 3. Уравнение вида (16) используется, например, для описания движения частиц в потоке газа в условиях скоростной неравновесности [7, 23-25]. При этом u — скорость частиц, V(t) — скорость газа (заданная функция), CR(t) — коэффициент сопротивления (обратно пропорциональный диаметру частиц). Графики зависимостей u(t), V(t), CR(t) для данного примера приведены на рис. 3. Коэффициент Cr изменяется в диапазоне значений от 0 до то. Здесь же и в табл. 2 по данным [25] приведены результаты расчета задачи с шагом h = n/(2N) в пяти точках интервала [0,п/2]. Видно, что схема Эйлера (Cr берется на i-м слое, V — на (i + 1)-м слое [25]) сильно размазывает решение u(t) и имеет очень большую погрешность в конце интервала. Специальная схема (11), ее аппроксимация (13), равно как и схема [23-25] второго порядка, с хорошей для практики точностью воспроизводят все особенности решения при грубом шаге интегрирования.

Заключение

Предложена специальная одношаговая разностная схема для численного решения жесткого ОДУ первого порядка. Схема имеет второй порядок аппроксимации и обладает равномерной по малому параметру при производной сходимостью. В основе ее построения лежит применение асимптотического метода Лапласа для приближенной оценки интегрального решения задачи Коши. Предложена рациональная аппроксимация второго порядка специальной схемы. Показана связь с другими известными одношаговы-ми методами. Численными расчетами подтверждена хорошая сходимость схемы и ее аппроксимации к точным решениям в широком диапазоне шагов интегрирования.

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

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

[2] Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Наука, 1987.

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

[4] Хайрер Э., НЁрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений: Нежесткие задачи. М.: Мир, 1990.

[5] Каханер Д., Моулер К., Нэш С. Численные методы и математическое обеспечение: пер. с англ. М.: Мир, 1998.

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

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

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

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

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

[11] Боглаев Ю.П., СтАнисловский А.А. Обзор библиографии по численным методам решения сингулярно возмущенных задач / АН СССР. ИФТТ. Препр. Черноголовка, 1984. 47 с.

[12] Боглаев Ю.П. О численных методах решения сингулярно возмущенных задач // Дифференциальные уравнения, 1985. Т. 21, № 10. С. 1804-1806.

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

[14] Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М.: Мир. 1999.

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

[16] Miller J.J., O'Riordan E., Shishkin G.I. Solution of singularly perturbed problems with-uniform numerical methods — introduction to the theory of linear problems in one and two dimensions // World Scientific, 1995.

[17] ЛисЕйкин В.Д., Петренко В.Е. Адаптивно инвариантный метод численного решения задач с пограничными и внутренними слоями. Новосибирск: ВЦ СО АН СССР, 1989.

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

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

[20] Самарский А.А. Введение в теорию разностных схем. М.: Наука, 1971.

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

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

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

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

[25] Рычков А.Д. Об одной разностной схеме, используемой при расчетах неравновесных и двухфазных течений // Газовая динамика. Томск: Изд-во Том. ун-та. 1977. С. 89-92.

[26] Зверев В.Г. Разностная схема третьего порядка точности для решения жесткого обыкновенного дифференциального уравнения с линейными коэффициентами // Вычисл. технологии. 2006. Т. 11, № 2. С. 28-38.

[27] Зверев В.Г. Разностные схемы повышенного порядка точности для численного решения жесткого обыкновенного дифференциального уравнения с линейными коэффициентами // Математическое моделирование. 2007. Т. 19, № 9. С. 94-104.

[28] Найфэ А. Введение в методы возмущений. М.: Мир. 1984.

[29] Де Брейн Н.Г. Асимптотические методы в анализе. М: ИЛ. 1961.

[30] Федорюк М.В. Метод перевала. М.: Наука. 1977.

Поступила в редакцию 3 октября 2007 г., в переработанном виде — 14 января 2008 г.