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

2010

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА Математика и механика

№ 4(12)

МАТЕМАТИКА

УДК 519.62

В.Г. Зверев

СПЕЦИАЛЬНЫЕ РАЗНОСТНЫЕ СХЕМЫ ДЛЯ РЕШЕНИЯ ОБЫКНОВЕННОГО ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ СО ЗНАКОПЕРЕМЕННЫМ КОЭФФИЦИЕНТОМ

Предложены специальные разностные схемы I и II порядка для сквозного расчета задачи Коши для обыкновенного дифференциального уравнения со знакопеременным коэффициентом при искомой функции. В основе их построения лежит использование точного интегрального решения уравнения и аппроксимаций сеточных функций, инвариантных к знаку коэффициента уравнения. Рассмотрены свойства, асимптотика и рациональная аппроксимация схем. Тестовыми расчетами показана возможность применения схем при грубых шагах интегрирования и сходимость численных результатов к точным решениям.

Ключевые слова: обыкновенное дифференциальное уравнение, первый порядок, знакопеременный коэффициент, специальные разностные схемы, сквозной расчет.

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

Как известно [1-6], для дифференциального уравнения вида u' = f(t, и) устойчивому по начальным данным решению отвечают отрицательные значения якобиана J = fu(t, u), неустойчивому - положительные значения. Несмотря на то, что предметом численного анализа традиционно являются устойчивые задачи, подобные уравнения, связанные с локальным ростом решения на отдельных участках области интегрирования, могут возникать при описании процессов в механике, химической кинетике, биологии и других приложениях [1, 4].

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

Цель работы - в развитие исследований [7, 8] разработать специальные разностные схемы I и II порядка точности для сквозного расчета задачи Коши для линейного неоднородного ОДУ первого порядка со знакопеременным коэффициентом при искомой функции.

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

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

ёи

Ьи(х) = е---+ а(х)и(х) = / (х), и(0) = и0, х > 0, х е [0, X], (1)

ёх

где е >0 - параметр из (0, 1]; а(х), _Дх) - гладкие функции; а(х) может менять знак на (0, X), причем области знакопостоянства а(х) известны. При а(х) > 0 имеем устойчивую, при а(х)< 0 - неустойчивую по начальным данным ветвь решения [2, 4].

Для иллюстрации рассмотрим пример [4], имеющий точное решение (рис. 1): и'+10(х- 1)и = 0, и(0) = Сехр(-5), х е [0,2], (2)

и(х) = Сехр(-5(х -1)2). (3)

Рис. 1. Точное решение уравнения (2) с областями устойчивости (1< х <2) и умеренной по начальным данным неустойчивости (0< х <1). Кр. 1-5 - С = 1;

2; 4; 6; 10

Кривые 1-5 этого семейства отвечают значениям С = 1; 2; 4; 6; 10 соответственно. На интервале хе(0, 1) они сильно расходятся. Здесь якобиан уравнения 3 = -10(х - 1) > 0, что соответствует области умеренной по начальным данным неустойчивости решения. При х > 1, 3 < 0, кривые сближаются, и имеет место область устойчивого решения.

Рассмотрим численное решение уравнения (1). На отрезке [0, X] определим произвольную сетку так, чтобы точки смены знака а(х) совпадали с узлами сетки {х,}, I = 0,1,...,N; х0 = 0, хы = X; = (хм -х{) > 0.

Выделим класс наиболее простых в реализации и применении одношаговых разностных схем. Для интегрирования (1) при а(х) < 0 и любом шаге И, можно ис-

пользовать явную схему Эйлера I порядка [1-6]:

и+1 = и(1 - є_1аЛ) + , (4)

где и, и+1 - известное и искомое значение и(х) на і и і+1 слоях по аргументу х. При смене знака а(х) > 0 и некотором значении і„ может измениться знак множителя (1 -е1аі Иі ), что ведет к осцилляциям решения [1]. Для его сохранения необходимо использовать неявную схему Эйлера I порядка:

и+1 = (иі + ) /(1 + 6-1 «г +Л ) , (5)

которая будет обладать устойчивостью для любого Ии так как (1 + еГ1аі+\И) > 1 [1]. Однако она не годится в случае а(х) < 0. Таким образом, каждая из схем напрямую не подходит для сквозного расчета.

На рис. 2 показано численное решение задачи при С = 1 по явной (4) и неявной (5) схемам Эйлера при постоянном шаге интегрирования И = 0,5 (а), И = 0,2 (Ь), И = 0,1 (с). Видно, что обе схемы (кривые 2 и 3) при грубом шаге теряют качество и подвержены осцилляциям (рис. 2, а, Ь). Кроме того, для неявной схемы на восходящей ветви решения имеет место вычислительная неустойчивость (кривая 3, рис. 2, с, и(1) = 18,6). Для обеих схем в силу I порядка характерна большая погрешность, поэтому требуется применение более мелкого шага, причем предпочтительнее вначале (0 < х < 1) выглядит явная схема (кривая 2 рис. 2, с[).

Рис. 2. Численное решение уравнения (2) при С = 1: а - И = 0,5; Ь - И = 0,2; с - И = 0,1; d -И = 0,025. Кр. 1 - точное решение; 2 (А) - явная; 3 (□) - неявная схема Эйлера; 4 - специальная схема (12) I порядка; 5 (о), 6 - специальная схема (15) и ее аппроксимация (27) II порядка

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

Специальная разностная схема I порядка

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

X X X

и(х) = и(х)ехр[-А(х)] +е- | /(%)ехр[-е-11а(п)ёг{\ё%, А(х) = е-11 а(%)^%. (6)

X % X

Оно не накладывает каких-либо ограничений на знак коэффициента а(х). Функция ехр[-А(х)] > 0 для Ух, поэтому ее аппроксимация не должна допускать смены знака при любых значениях сеточных параметров. Эти особенности должны прослеживаться и в разностном аналоге (6).

В схеме Эйлера используются кусочно-постоянные коэффициенты: а = аг,

/ = / для явной (4) и а = аг+1, / = /+1 - для неявной схемы (5). При этих условиях (6) легко интегрируется и позволяет получить распределение и(х) на сеточном интервале:

/

и(х) = и(хг)ехр[-е-1а(х-хг)] + — {1 -ехр[-е-1а(х-хг)]} . (7)

а 1 ’

Полагая в (7) х = х1+1, получим разностное выражение для вычисления искомого и1+1:

ui+1 = ui ехр(-2) + е-1Лг /{—ех^^ Г)}, г = е-1ййг. (8)

Здесь 2 - сеточный параметр. Видно, что (8) справедливо для любого знака г (определяется знаком а ). Следует отметить, что для устойчивой задачи (г > 0) решение (8) в [10, 11] положено в основу равномерной по параметру е разностной схемы I порядка.

Явная схема Эйлера (4) следует из (8) при традиционном разложении экспоненты:

а = а, / = /г, ехр(-) ~1 - , =е-1 аЛ-. (9)

При 2^ < 0 (а,- < 0) имеем ехр(—2) и 1 + 21. При больших 2^ > 1 (а,- > 0, жесткий случай, малые е) данная аппроксимация непригодна из-за смены знака выражения (9).

Неявная схема Эйлера (5) вытекает из (8) при следующем приближении экспоненты:

а = аг/ = ехр(-2, +1) ~ 1/(1 + 2+1X 2+1 =е~Ч+1Л . (10)

Это дает асимптотику редуцированной части задачи (1) иг+1 = /+1/аг+1 при 2 ^ да, однако не подходит при 2 <0 (а(х)<0). Таким образом, традиционный вид записи аппроксимации экспоненты (9), (10) препятствует проведению сквозного расчета и требуется его модификация.

Выделим для этой цели знак сеточного параметра г = х) X и запишем (8)

в виде

Мг+1 = и [ехр(| х| ХТ^М +в~1И/ ~ [^|7(Х) } • (11)

Используем разложение экспоненты в ряд Тейлора с I порядком

ехр(|Г|) и 1 + IX , получим

и+,= и[1+1 х1Г- х > +.-*/11 ~ ^ ^ }• "2>

Коэффициенты а, 7 определим по левому (для явной схемы) или по правому

(для неявной схемы) концу сеточного интервала, например:

[(а, 7 X (а ^0, а+1 ^0),

(а, 7Х 1 (а+1, 7+1), (а ^0, а+1 ^0).

Из (12) при г < 0 следует явная (4), при г > 0 - неявная схема Эйлера (4). Случай Х| ^ 0, согласно (8), (12), в пределе дает щ+х = и1 + .~1к{7, поэтому определим его и при г = 0.

Результаты сквозного расчета по схеме (12) при различных шагах интегрирования показаны кривыми 4 на рис. 2. Видно, что отсутствуют осцилляции, кривые качественно соответствуют точному решению (кривая 1), но в силу I порядка имеют большую погрешность. Для количественного соответствия требуется сильное уменьшение шага, рис. 2, ё. Погрешность схемы (12) связана с приближенным вычислением экспоненты и использованием кусочно-постоянных коэффициентов. Чем больше сеточный параметр |х| > 1, тем хуже воспроизводится функция е(х) = ехр(—х), рис. 3. Поэтому для повышения точности численных результатов при грубых шагах интегрирования требуется применение других аппроксимаций.

Рис. 3. Сеточная функция е(х) = ехр(-х) и ее аппроксимации I и II порядков. ег = 1-х; е1 = 1/(1+х); е2 = [1+И+И2/2]-явп(г)

Специальная разностная схема II порядка

Рассмотрим интегральное решение (6) уравнения (1) и выделим случай а(х) Ф 0 . Используем подход [8], основанный на применении асимптотического метода Лапласа [12-14] для вычисления интегралов с экспоненциальной функцией. Применяя к (6) интегрирование по частям, получим

и(х) = и(х) ехр[-А(х)] + /(х) - /(х) ехр[-А(х)] -

-1

Ш

.а©.

ехр

а( х) а( хі)

X

-є-11 а(п)й П

(13)

Следует заметить, что при кусочно-постоянных а, / интеграл в (13) равняется нулю и сразу следуют выражения (7), (8). Возьмем в подынтегральном выражении сомножитель [Д|)/а(|)]7а(|) в средней точке 0.5(х, + х) и точно вычислим оставшуюся часть интеграла с экспонентой. В результате получим приближенную зависимость и(х) на сеточном интервале:

и(х) = и(х) ехр[-А(х)] + /(х) - /(Х) ехр[-А(х)] -

а( х) а (х)

/ (х) - /(х)

а( х) а( хі)

(1 - ехр[- А( х)]) 'I

(14)

є 'а(х-х)

Г

Положим в (14) х = х,+1 и вычислим А(х ,+1) по формуле трапеций, которая точна для линейной функции а(х):

х+1

А(х+!) = е-1 | а(п)ёП ~ е1ак,, а = (ai + а+) / 2 .

В результате получим специальную одношаговую схему:

/

и+1 = иів 2 + -!^1[1 -Р(г)] + ^~[Р(2)-е 2П, Р(2) =

/

(1 - е-2 )

2 = є Іакі. (15)

Здесь Р(2) - сеточная функция [7]. Схема (15) для устойчивой жесткой задачи (2^-ю, є^-0) рассматривалась в [8]. Она также не имеет ограничений на знак коэффициента а . Вводя сеточные функции ^(2), п(2) [7, 8], (15) можно переписать в другом виде:

(16)

е(2)=Є2, §(2)=

1 -Р(2) [е~2 + 2 -1]

П( 2) =

Р(2) -е-2 [1 - е- (2 + 1)]

2 2 2 2"

Вид функций е(2), Р(2), |(2), п(2) и их аппроксимаций показан на рис. 3-5. Случай а(х) = 0 и переход а(х) через ноль требуют отдельного рассмотрения. Запишем решение (6) для и м на сеточном шаге в окрестности нуля а(х) в следующем виде:

“г+1

"і+і

= иі ехр[-А(х+1)] +є-1 ехр[-А(х+1)] | /(|)ехр[АОй\ . (17)

х

Рис. 4. Сеточная функция Р(2) = [1 - ехр(-2)]/2 и ее аппроксимация II порядка.

Р2 = [1 - б2(2)]/2

Рис. 5. Сеточные функции ^(2), п(2) и их аппроксимации (^2, п2) II порядка. ^2 = [62(2) + 2 - 1]^2; П2 = Р2 - ^2.

При а, = 0, а1+1 = 0 из (17) со II порядком следует

и,+1 = и, +е-1/мпк1, /,+1/2 = [/+, + /]/2 . (18)

Рассмотрим случай а, = 0, а1+1 Ф 0. Коэффициент а(х) можно хорошо приблизить линейной функцией

а(х) = к(х - х,), к = а1+1 / А,, х е [х,, х1+1 ].

Подставляя а(х) в (17), получим

и,+, = и,

Л1 +1

ехр(-2) +е-1ехр(-2) | /(5)ехр[2(§-х)2/А,2]ё5,

2 =

е 'а+Л

(19)

Выделим знак z = sign(z)|z|, возьмем f(Q в средней точке, преобразуем (19) к виду

S 1 f+1/-Д. exp[-sign(z) | z |] f 2

uM = ut exp[-sign(z) | z |] +-г+^2-г—-j=-------------- f exp[sign(z)y ]dy . (20)

VI z | 0

Отсюда, при z >0 (a г =0, a г+1 > 0) следует

Иг+1 = Ыг exp(-z) +8-1f+1/2hJ(z), J (z) = £XP(JL) f ey2 dy , (21)

<z 0

X

где J(z) выражается через интеграл Досона D(х) = exp(-x2)fexp(t2)dt [15]. Для

0

растущей ветви решения при z < 0 (аг = 0, аг+1 < 0) (20) принимает вид

Ыг+1 = Ыг exp(| z |) + s^f+1/2hG(z),

exp(| z |)^ -y2 exp(| z D-Jn in; (22)

G(z) = A J e y dy = A -TerfW z).

Vkl 0 vkl 2

2 х

Здесь erf( x) =-=f exp(-t 2)dt - интеграл вероятностей [15]. y/n 0

Рассмотрим другой случай приближения коэффициента а(х) к нулю, когда аг Ф 0, аг+1 =0. Запишем а(х) на сеточном интервале в виде

а(х) = аг (1 - (х - хг) / h), х е [хг, хм ].

Подставим а(х) в (17), возьмемf(^) в средней точке, выделим знак z = sign(z)|z|. В результате получим

е—1 f h ''i\z\

Ыг+1 = U exp[-sign(z) | z |] +S—г!!!-2 г J exp[-sign(z)y2 ]dy, z = s-^. /2 .

v| z | 0

Отсюда, при z >0 (а г > 0, а /+1 = 0) следует разностное выражение

1 / rf ( / )

Ыг+1 = Ыг exp[-z] +s^lf+l/2hlK(z), К(z) = -= f exp[-y2]dy = —. (23)

Vz 0 2 Vz

В другом случае при z <0 (соответствует аг <0, аг+1 =0) разностное выражение принимает вид

1 2

Ыг+1 = Ыг exp[| z |] +s~lf+1/2hL(zX L(z) = -/= f ey2dy . (24)

v| z | 0

Таким образом, расчетная формула (15) для определения и г+1 при а(х) Ф 0 дополняется выражениями (18), (21)-(24), описывающими различные случаи перехода а(х) через ноль.

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

Рассмотрим свойства специальной схемы (15), (16). Легко показать, что при а(х)=сош^ Ах) - линейной функции, а также другом важном случае - f.х)/а(х) = = const, а(х) - линейная функция - схема дает точное решение задачи. Разностные выражения (21) - (24) при переходе коэффициента а(х) через ноль будут точны

для линейной функции а(х) и А(х) = const. В общем случае для произвольных а(х), Ах) она будет иметь II порядок.

При z—да (а > 0) функции |(z) ~ 1/z, n(z) ~ 1/z2 и решение (16) асимптотически выходит на и г+1 ~ А+1/а г+1. Таким образом, для жесткой неоднородной задачи при а(х) >0 преобладающим является весовой множитель у источника на г+1 слое. При отрицательных z — -да функции |(z) ~exp(|z|)/z2, n(z) ~ exp(|z|)/|z| и решение (16) имеет экспоненциальный рост иг+1 ~ (иг -f /а;)exp(|z|). Из этого следует, что в неоднородной задаче доминирует множитель при источнике на г-м слое.

При |z| ——0 (а Ф 0, h/s—0), |(z) — 1/2, n(z) — 1/2 и схема (16) принимает вид, аналогичный формуле трапеций:

Ыг +1 = Ы г + ^

■ 1_ г+1 г J

В точке перехода коэффициента а(х) через ноль при |г| ^0 (А/е^-О) выражения (21) - (24) стремятся к

иг +1 = иг + ^/г+1/2Ьг .

Представляет практический интерес рациональная аппроксимация сеточных функций в (15), (16), независимая от знака аргумента г (коэффициента а,+1/2). Запишем их с учетом г = Б1§п(г)|г| и разложения ехр(|г|) и 1+Х|+Х|2/2 со II порядком:

Г . . . ,9 г) 1 — 6 (г) е (г) + г — 1

= |1 + И + И /2] , Р2ОО =-2--, §2(г) = “-------2-----. (25)

Л 2 г2

Непосредственной проверкой приходим к выводу, что в схеме (15), (16), в конечном счете, имеют место следующие разложения сеточных функций (рис. 3-5):

z > 0: ^2 (z) =----—2—, Р2 (z) = {1 + ~}^2 (z),

1 + z + z2 / 2 I 2)

^2(z) = e2(zX П2(z) = 2e2(z); (26)

z < 0: e2(z) = 1 + |z| + |z|2/2, P2(z) = (1 + |z|/2),

1 (1 + 1 z|)

^2( z ) = ^ П2 (z) =------—.

Подставляя (25) в (15) или (16), получим следующий вид разностного выраже-

ния:

Ыг + 2 {(f+1 / аг+1 )(1 + z) + (f / аг )}

z > 0: иг+1 =-----------2------------------------------------------- -; (27)

г+: (1 + z + z 2 /2)

.12 , \z\ J f+1 . f

z ^ 0: иг+1 = иг (1 + l A +lzl /2)+ V Ir^ + Or (1 + l zl)f .

2 иаг+11 |аг| J

В случае а = const при z >0 (27) совпадает со схемой, приведенной в работах [7, 16, 17].

Рассмотрим аппроксимацию схемы (21) - (24) в окрестности перехода коэффициента а(х) через ноль. В случае (а г =0, а г+1 > 0) из (21) со II порядком следует

2 > 0:

Щ+1 (1 + 2 + 2 2/2)

+ 6 /г+1/2 Иг32( 2 ),

-1

32 (2) =-

1 + 2/3

2 =

Б- аг +А

(28)

1 + 2 + 22/2 2

Для возрастающей ветви решения (а г =0, а/+1 < 0) для (22) можно использовать следующий вариант знакопостоянной аппроксимации:

2<0: и+1 = и(1+1 А+1212/2)+е1/+1/2л^2(2х С2(2) =1+|2^+12/,/2

г 1+р|/з

В случае приближения а(х) к нулю, для а г > 0, а /+ = 0 из (23) следует

2 > 0: и+1 = -

- + В-1/+1/2 ^2(2 ), К 2 (2 ) =-

1

(29)

(30)

(1 + 2 + 22/2) ■""1'2 ' ^ ^ 1 + 2/Г 2

Соответственно для (аг < 0, аг+1 = 0) из (24) со II порядком получим

2 < 0: и+1 = и(1+ | 2 | +1 212 /2) +6-1/+1/2Л^(2), 2) = (1+ | 2 |/3) . (31)

Аппроксимации (28)—(31) неоднородного слагаемого в окрестности смены знака коэффициента а(х) показаны на рис. 6. Нетрудно видеть, что для устойчивой ветви решения 32(2), К2(2) хорошо подходят до 2 ~ 5 (рис. 6, а). Для экспоненциально растущей ветви решения уже трудно ожидать хорошей точности, 02(2), Ь2(2) годятся до 2 ~ 1, рис. 6, Ь.

Рис. 6. Аппроксимации (28) - (31) неоднородного слагаемого в окрестности смены знака коэффициента а(х). Кр. 1 - 4 - точные значения 3(2), К(2), 0(2), Ь^); кр. 1’- 4’ - их аппроксимация II порядка

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

Результаты численного решения задачи (2) на основе специальной схемы (15) и ее рациональной аппроксимации (27) показаны на рис. 2 кривыми 5 (кружочки) и 6 (пунктир). Так как для линейного коэффициента а(х) аппроксимация Л(х1+1) является точной, то точными получаются и результаты по специальной схеме. Аппроксимация экспоненты, несмотря на II порядок, вносит свою погрешность (кривые 6), которая становится все заметнее с ростом шага И. Тем не менее на всех рисунках кривая 6 выглядит гораздо лучше, чем кривая 4 для сквозной схемы Эйлера I порядка.

Рассмотрим пример из работы [10]:

ем'(х) + (1 + х)и(х) = (1 + х), и(0) = 0, х е [0,2], (32)

и (х) = 1 - ехр(-(2х + х2) /(2е)),

В аналитическом решении положим е = -1 (рис. 7). В табл. 1 приведены абсолютная и относительная ошибки экспоненциально убывающего численного решения иИ при различных И по схеме (8), причем а = а, / = ^ [10], по специальной схеме (15) и ее аппроксимации (27). Для сравнения по данным [10] приведены результаты решения задачи по неявной схеме Эйлера и специальной схеме Титова, Шишкина (постоянный подгоночный параметр), имеющей равномерную по е сходимость О(И).

Рис. 7. Ошибка численного решения задачи (32) при є = -1. Н = 0.1. 2 - схема Эйлера [10]; 3 - схема Титова, Шишкина [10]; 4 - схема (8) [11]; 5, 6 - специальная схема (15) и ее аппроксимация (27) II порядка

Погрешность численного решения задачи (32) при є = -1

Вид схемы шах| и( х) - иН (х)| 0< х<2 шах 0< х<^ и(х) - иН (х)

и( х)

є = -1 Н = 1 Н = 0,1 Н = 0,01 Н = 1 Н = 0,1 Н = 0,01

Схема Эйлера [10] - 26,0 1,88 - 0,486 3,50-10-2

Схема Титова, Шишкина [10] 64,89 8,89 0,74 1,21 0,166 1,38-10-2

Схема (8) [11] 34,51 5,2 0,543 0,644 9,69-10-2 1,01-10-2

Аппроксимация (27) схемы (15) 30,58 1,5 1,79-10-2 0,571 2,8-10-2 3,33-Ю-4

Схема (15) 4,44-10-16 2,84-10-14 3,55-10-14 1,28-10-16 6,27-10-16 1,08е-14

Из рисунка и таблицы видно, что схема (15), как и следует теоретически, решает эту задачу с точностью до машинного представления числа. Вполне приемлемой с практической точки зрения выглядит и ее аппроксимация (27) II порядка. Здесь максимальная относительная ошибка при И = 0,1 составляет 2,8 %, в то время как по схеме Эйлера - 48,6 %, схеме [10] - 16,6 %, схеме (8) [11] - 9,7 %. Следует заметить, что при грубом шаге И = 1 неявная схема Эйлера вообще не имеет решения, а схема [10] с постоянным подгоночным параметром допускает частичную потерю качества.

Рассмотрим задачу с более сложной зависимостью коэффициентов уравнения: и'(х) + лсо8(лх)и(х) = {лсо8(лх)-2(х-2}}ехр(-(х-2)2), хе[0,4], (33)

и(0) = 1 + ехр(-4).

Аналитическим решением (33) является функция

и(х) = ехр(- Бт(лх)) + ехр(-(х - 2)2).

В точках х = 0,5; 1,5; 2,5; 3,5 происходит смена знака коэффициента при искомой функции в левой части уравнения. В специальной схеме (15) в этих узлах дополнительно используются выражения (21-24), а для ее аппроксимации (27) - выражения (28-31).

Численные расчеты при грубом шаге к = 0,25 показаны на рис. 8. Видно, что при этих условиях результаты по явной схеме Эйлера (кривая 2) лишь отдаленно напоминают точное решение (кривая 1), неявная схема (кривая 3) имеет сильную неустойчивость. В противоположность этому гораздо лучше выглядит сквозная схема Эйлера (светлые ромбики 4). Однако в силу I порядка и накопления ошибки имеет место сглаживание хода кривой. Результаты по специальной схеме II порядка (светлые кружочки 5) и ее аппроксимации (символ *) в количественном отношении хорошо воспроизводят все особенности хода точного решения.

Рис. 8. Численное решение уравнения (33): 1 - точное решение; 2 (А) - явная схема Эйлера; 3 (□) - неявная схема Эйлера;

4 (◊) - специальная схема (12) I порядка; 5 (о), 6(*) - специальная схема (15) и ее аппроксимация (27) II порядка

Заключение

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

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

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

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

ЛИТЕРАТУРА

1. Самарский А.А. Введение в численные методы. М.: Наука, 1982. 272 с.

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

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

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

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

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

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

8. Зверев В.Г., Гольдин В.Д. Об одной специальной разностной схеме для решения жесткого обыкновенного дифференциального уравнения // Вычислительные технологии. 2008. Т. 13. № 3. С. 54-64.

9. Камке Э. Справочник по обыкновенным дифференциальным уравнениям. М.: Наука, 1976. 576 с.

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

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

12. Найфэ А. Введение в методы возмущений. М.: Мир, 1984. 535 с.

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

14. ФедорюкМ.В. Метод перевала. - М.: Наука, 1977. 369 с.

15. Абрамовиц М., Стиган И. Справочник по специальным функциям. М.: Наука, 1979. 832 с.

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

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

Статья принята в печать 28.09.2010 г.

Zverev V.G. THE SPECIAL DIFFERENCE SCHEMES FOR NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATION WITH THE ALTERNATING-SIGN COEFFICIENT. The special difference first and second order schemes for the through numerical calculation of the Cauchy problem for an ordinary differential equation with the alternating-sign coefficient at the unknown function are proposed. The construction of schemes is based on the use of the exact integral solution of the equation and approximations of the grid functions invariant to the sign of the equation’s coefficient. Properties, asymptotics, and the rational approximation of the schemes are considered. An opportunity of application of the schemes at rough steps of integration and convergence of numerical results to exact solutions is shown on test examples.

Key words: ordinary differential equation, first order, alternating-sign coefficient, special difference schemes, through numerical calculation

ZVEREV Valentin Georgievich (Tomsk State University)

E-mail: [email protected]