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

УДК 517.91

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

В.В. Коробицын, Ю.В. Фролова

A novel method for numerical solving autonomous ordinary differential equation dx/dt = f (ж) is suggested. The solution is found as a function t(x) for that the incrementation is independent from previous points of solution trajectory. This fact is used for parallel realization of computational method.

The numerical experiments show the applicability of suggested method for area with rapid acceleration of right-hand function.

Численные методы решения обыкновенных дифференциальных уравнений появились довольно давно. Еще Ньютон и Эйлер предложили методы поиска приближенного решения. Классические работы Рунге, Куп ы. Адамса посвящены построению одношаговых и многошаговых методов решения ОДУ. Казалось бы, всё уже разработано и внедрено, однако интерес к разработке новых и модификации существующих методов не ослабевает. С одной стороны, новые задачи требуют новых методов. Так, в последнее время активно развиваются методы решения жестких, дифференциально-алгебраических и разрывных задач [1]- [4], [8]- [11]. А с другой стороны, развитие компьютерной техники требует разработки эффективных методов решения на многоядерных, многопроцессорных и распределенных компьютерных системах. В последнее время стали появляться работы, связанные с параллельной реализацией методов решения ОДУ [5]- [7], [12].

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

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

Copyright © 2009 В.В. Коробицын, Ю.В. Фролова.

Омский государственный университет.

E-mail: [email protected]

1. Аппроксимация решения автономного ОДУ

Рассмотрим задачу

йх ,

— = ;{х), х(Ь0) = х0, Ь> Ь0. (1)

Необходимо найти решение х(£) : К ^ К задачи Коши (1), учитывая, что функ-

ция правой части / является непрерывной по переменной х и те зависит от £ в явном виде.

Решение х(£) можно разложить в ряд Тейлора вблизи точки ¿о

йх й2х т2

х(г0 + т) = ж(*о) + ^(¿о) • Г + ^-(¿о) ■ • (2)

Частичную сумму ряда (2) обозначим за Рк(£о,т),

йх й2 х т2 йк х т к

у . у . у ч I у ч I

(*0’т) := д(ад ■т + да ^ ■ 2Т + ■ ■ ■+ ' ¥■

тогда

х(£о + т) - х(£о) = Рк(¿о, т) + 0(тк+1).

Введем обозначения

/о = /Ы, /' = ^Ы, /о = ^гЫ

и приведем выражения для Рк, к = 0,1, 2,...:

Ро(£о, т) = 0, йх

Р1(г0,т) = Р0(г0,т) + ^-(¿о) • т = ¡от,

й2х т2 т2

Р2^о,т) = А(£о,т) + (¿о) • = /о'Г + /о/оу,

й3х т3 т2 т3

Рз&о,т) = Р2^0,т) + ^3 (¿о) • = /о'Г + /о/о — + (/о/о + /о2/о)“^г> • • •

т

значение решения х(£о + т ), мы зададим малое прир ащение йх то перемен ной х т

х(£о + т) - х(£о) = йх.

В этом случае необходимо решить алгебраическое уравнение

Рк (¿о,т )= йх. (3)

Функция Рк(¿о,т) является многочленом от т степени к, поэтому уранение (3) будет иметь к корней. Необходимо найти наименьший действительный положительный корень, который будет давать значение приращения по переменной необходимое для достижения точки хо + йх.

Отметим, что решение уравнения (1) перейдет из точки ж0 в точку ж0 + dx при dx > 0 только когда f (ж) > 0 для всех ж £ [ж0; ж0 + dx]. Однако если в некоторой точке ж1 этого интервала функцня f (ж1) = 0, то это означает, что в этой точке уравнение (1) имеет особую точку. Поэтому решение не достигнет точки ж0+dx. Если же f (ж0) < 0, то изучается интервал [ж0; ж0 — dx] аналогичным образом,

2. Формулы для вычисления т

Далее будем считать, что f (ж) > 0 для всех ж £ [ж0; ж0 + dx]. Введем обозначения

fi = f(ж0 + dx), f-i = f(ж0 — dx).

f ж0

лами

Й « /¿(4) := /о ”/о№) ■= /l~2^ + /-1. (4)

Оба приближения имеют погрешность порядка O(dx).

Найдем аналитические выражения для корней уравнения (3) при разных k. При k = 1 получаем уравнение f0T = dx. Корень этого линейного уравнения определяется выражением т = Погрешность решения, соответствующего найденному т, будет порядка O(dx).

При k = 2 получаем квадратное уравнение относительно т

f0 f0 2 I г 1 n

-^-T + for -cix = 0,

корни которого определяются выражениями

D = /о2 + 2/„/Я. = 4 (-1 ± Jl+2d.,f).

f0f0 f0 V V f0/

Если f0 > 0, то положительный действительный корень будет

I-1)- (5)

Если f0 < 0, то для существования действительного корня необходимо удовлетворение условия 1 + 2dxfg /f0 > 0 что будет выполнен о при dx < — f0 /(2f^), Учитывая выражение для f из (4), получаем dx < —f0dx/(f1 — f-1). Отсюда

— f0 / ( f1 — f- 1 ) > 1 dx

полнено условие f1 — f-1 > — f^. Тогда минимальное т для ^^стижепия ж0 + dx

соответствует (5),

При k > 3 получаем алгебраическое полиномиальное уравнение относительно т степени k. Это уравнение можно решать численно. При этом необходимо

т

сечению кривой полинома с точкой ж + dx.

Однако мы предлагаем повысить точность вычисления т, не повышая степени полинома. Для этого используем разложение функции решения в точке х1 = £(¿1), ¿1 = ¿0 + т, ВЫЧИСЛЯЯ значение функции решения В точке х0 = х1 — 4, получаем

т2 ^к X тк

,((0) = *(« - Ж((1). г + ¥ + • • • + (-1)‘^((1) . ¥ + . . .

Тогда х(£0) — х(^) = Рк(£ь —т) + 0(тк+1), Учитывая выражение (3), получаем два соотношения х1 — х0 = Рк(£0,т) и х0 — х1 = Рк(¿1, — т), Вычитая эти

2dx = Pk(to, т) - Pk (ti, -т). (6)

Введем обозначения

î ^ и df / ^,/2-/0 ,// ^,/2-2/1 +/о

/2 = /(жо + 24), Л = у№ ~ —, Л = № ~-----«-----• (7)

dx 2dx dx2 dx

т

параметра k.

При k = 1 получаем линейное уравнение 2dx = /0т — fi • (—т), Решение этого уравнения т = 2dx/(/0 + /1 ), Погрешность определения величины т будет порядка O(dX).

При к = 2 получаем квадратное уравнение

(/о/о — /i/i) -у + (/о + /i)t — 2dx = 0. (8)

При к > 3 получаем полиномиальное уравнение относительно т, наименьший положительный вещественный корень которого можно вычислять численно.

Далее мы найдем аналитическое выражение для наименьшего положительного вещественного корня уравнения (8) и оценим погрешность его вычисления

/

ми разностями вида, (7),

Выражение для коэффициента уравнения при квадрате можно вычислить

t t, t _ f\~ /-1 t /2 — /0 t _ —f-ifo + 2/o/i — /1/2 J0J0 - J1J1 - wi J0 wi J1 - wi •

2dx 2dx 2dx

Найдем корни уравнения (8)

D = (f0 + fi) — 4(/0/0 — f1 fi) ( —2dx) = (f0 + fi) + 4( — f-i/o + 2/0/1 — /1/2) , _-(/„ + /1)±v/D_J (Л + /1) T \/(/о + Л)2 - 4(/-i/o - 2/0/1 + /1/2)

rl,2 — —Г7777----— — dx

2(/0/0 — /i/i) x (/-i/0 — 2/0/1 + /1 /2 )

Введем обозначение

/с + /і

/-1 /с — 2/С/і + /і/2

и перепишем выражение для корней уравнения

Для существования вещественных корней уравнения необходимо выполнение условия -7--4 ч < 1, это будет выполнено в одном из случаев (/о + /0^ < о

Наименьший положительный вещественный корень всегда определяется выражением

Причем если произведение > 0, то должно быть /0 + /0 ^ > 4, а если < 0, ТО ДОЛЖНО быть /0 + /0 ^ < 0,

Каким должен быть 4? Знак величины ¿х определяется этаком функции /, Если /0 > 0, то должно быть ¿х > 0, Длина шага ¿х будет зависеть от оценки погрешности.

Оценим погрешность найденного корня (9), При вычислении корня произ-/

чпсленпя которых порядка О(^Х). Следовательно, погрешность вычисления ^ будет порядка О(^Х), а погрешность определения т — порядка О(^Х).

3. Практическая оценка погрешности

т

числения с двумя шагами длиной ¿х и одним шагом длиной 2^х, тогда

где 71, т2 — два значения, вычисленные при шаге 4, т3 — одно значение, соответствующее шагу 2^х, Вычитая первое выражение из второго, получим

УиШБИЙ 7-----ч- ^

(л+л) ^ или (/с + /1)^ > 4,

(9)

¿(хс + 2^х) - ¿с = Ті + Т2 + С • (4)р

¿(хс + ЭД — ¿0 = Т3 + С • (24)Р

Т1 + Т2 — Т3 = С • (4Г • (2? — 1).

Тогда погрешность вычисления т с шагом ¿х будет

Егг^ = С • (4)р

Т~1 + т2 - Тз 2Р - 1

4. Численные эксперименты

Точность предложенного метода нахождения численного решения автономного ОДУ оценим па двух примерах.

Первое уравнение — линейное:

&х

— = х, х(0) = хо, I > 0. (10)

Аналитическое решение этого уравнения определяется формулой х(£) = х0е*, Откуда можно выразить функцию ¿(х) = 1п(х/х0).

Дня численного эксперимента использовалась формула (9), Были заданы параметры: х0 = 1, йх = 0.19, График решения уравнения (рис, 1-а) показывает экспоненциальный рост х(£). Порядок относительной погрешности численного решения уменьшается от -1,833 до -2,845 (рис, 1-6), Судя по данному графику, можно сделать вывод, что предложенная схема вычисления решения лучше работает в условиях больших значений функции правой части уравнения.

\gBerr

а) б)

Рис. 1. а) График решения системы (10). б) Порядок относительной погрешности

численного решения

Второе рассматриваемое уравнение — модель Верхюльста-Пёрла:

(^х

— = к (а — х)х, ж(0) = жо, I > 0. (11)

Аналитическим решением уравнения является логистическая кривая, определяемая формулами:

= I + екаа-Сг) ЦРИ 0 < * < = I _ еМ*-С2) 1фИ ^ > а-

Значения констант определяются исходя из начальных данных

— 1п — 1п ■

/~1 __ _______а—хо _ ______хр—а

' ка ' 2 ка

Численное решение построено при заданных параметрах: а = 3, к = 2, ^х = 0.04, Начальные данные для построения двух траекторий решения: х° = 1, х° = 5, Первая траектория (рис, 2-а) помечена квадратиками, вторая — кружками,

X

а) б)

Рис. 2. а) График решения системы (11). б) Порядок относительной погрешности численного решения; линии с квадратиками для 0 < х < а, линии с кружками для х > а.

Уравнение (11) имеет особую точку х = а, при приближении к которой рост времени £ приводит к малым изменениям переменной х. Поэтому в нашем случае, когда шаг йх задан постоянным и вычисляется шаг по времени т, при приближении к особой точке происходит интенсивный рост величины т, что приводит к повышению погрешности решения. Это можно увидеть па обеих кривых порядка относительной погрешности решения (рис, 2-6),

Полученные результаты вычислений показывают рациональность использования предложенного метода дня поиска решения в области интенсивного роста значений функции правой части уравнения. Поэтому разумно предположить эффективность использования этого метода дня решения жестких уравнений. Главным преимуществом представленного метода является его параллельность — на каждом г-м интервале вычисление шага т можно проводить параллельно, после чего решение можно составить х(£ + ^™=1 т{) = хо + пАх:. А также дня нахождения одного шага требуется только одно повое вычисление функции правой части, что но вычислительным затратам сравнимо с методом Эйлера, по точность предлагаемого метода значительно выше,

5. Идея об определении траектории решения автономной системы ОДУ

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

^ = /(*1.*2),

*%- = ФъХг), „.-л

х?(£°) = х?,

Х2(£о) = х°.

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

I ТГ~^~~\= I М => х2) = [ г^Х1 \ + съ

3 / (х?,х2) 3 ] / (х?,х2)

из второго уравнения

г2(х1,х2)= [ 2 + С2.

о 9 (х1, х2)

Введенные здесь функции £?(х?, х2) и £2(х?, х2) описывают поведение интегралов от обратных функций правых частей системы (12), Константы С? С2 определяются исходя из начальных условий системы.

Предполагаем, что решение системы (12) можно определить как параметрическую кривую Г в пространстве М2 вида

Г = {(х?,х2) е М2| ¿1 (х?,х2) - ¿2(х1,х2) = 0}.

Параметр в кривой Г определяется в = ¿1(х1,х2) = ¿2(х1,х2), (х1,х2) е Г,

Для нахождения численного решения системы (12) необходимо разбить фазовое пространство системы на квадратные участки Бу, (г е ^, 3 е Z) со стороной 4:

Б^ = {(х1,х2) е М2| х° + г^х < х1 < х1 + (г +1)^х,х2 + 34 < х2 < х2 + (3 + 1)4}.

Для каждого квадрата Б^ вычисляются значения функций правых частей системы в углах этого квадрата. Используя билинейную интерполяцию, вычисляются интегралы для определения ¿1; ¿2. Зная точку входа траектории решения

С1 С2

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

Отметим, что вычисление интегралов в каждом квадрате Б^ является независимым от значений функций в других квадратах, что дает возможность производить вычисления в параллельном режиме. Определение конкретной траектории решения, соответствующей заданной начальной точке, будет основано на определении точек выхода для каждой клетки, для чего потребуется находить С1 С2

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

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

И Т.Д.

Предложенную идею определения траектории решения систему ОДУ можно распространить на п-мерный случай. Траектория будет определяться как пересечение п — 1 гиперповерхностей размерноети п — 1,

Литература

1. Деккер, К. Устойчивость методов Рунге-Кутты для жестких нелинейных дифференциальных уравнений / К. Деккер, Я. Вервер. - М.: Мир, 1988.

2. Коробицын, В.В. Алгоритм численного решения дифференциальных уравнений с разрывной правой частью / В.В. Коробицын, Ю.В. Фролова // Математические структуры и моделирование. - 2005. - Вып. 15. - С.46-54.

3. Коробицын, В.В. Исследование поведения явных методов Рунге-Кутты при решении систем обыкновенных дифференциальных уравнений с разрывной правой частью / В.В. Коробицын, В.Б. Маренич, Ю.В. Фролова // Математические структуры и моделирование. - 2007. - Вып. 17. - С. 19-25.

4. Коробицын, В.В., Алгоритм численного решения кусочно-сшитых систем / В.В. Коробицын, Ю.В. Фролова, В.Б. Маренич // Вычисл. технологии. - 2008.

- Т.13, N.2. - С.70-81.

5. Новиков, Е.А. Реализация полуявного (4,2)-метода решения жестких систем на параллельной ЭВМ / Е.А. Новиков, Л.П. Каменщиков // Вычисл. технологии. -2004. - Т.9. Вестник КазНУ им. аль-Фараби. Сер. Математика, механика, информатика. N.3 (42). (Совм. выпуск. 4.1). - (’.235 2 11.

6. Фельдман, Л.П. Сходимость и оценка погрешности параллельных одношаговых блочных методов моделирования динамических систем с сосредоточенными параметрами / Л.П. Фельдман // Научные труды ДонНТУ. Серия: Информатика, моделирование и вычислительные методы, (ИКВТ-2000), выпуск 15. - Донецк: ДонНТУ, 2000. - С.34-39.

7. Фельдман, Л.П. Параллельные алгоритмы численного решения задачи Коши для систем обыкновенных дифференциальных уравнений / Л.П. Фельдман, И.А. Назарова // Матем. моделирование. - 2006. - Т.18, N.9. - С.17-31.

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

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

10. Cash, J.R. A variable order Runge-Kutta method for initial value problems with rapidly varying right-hand sides / J.R. Cash, A.H. Karp // ACM Trans. Math. Software. - 1990. - V.16, N.3. - P.201-222.

11. Gear, C.W. Solving ordinary differential equations with discontinuities / C.W. Gear, O. 0sterbv // ACM Trans. Math. Software. - 1984. - V.10. P.23-44.

12. Jackson, K. The potential for parallelism in Runge-Kutta methods. Part I: RK formulas in standard form / K. Jackson, S. Norsett // SIAM J. Numer. Anal. - 1995. - V.32. -P.49-82.