close

Вход

Забыли?

вход по аккаунту

?

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

код для вставкиСкачать
Математические
структуры и моделирование
2009, вып. 19, с. 1927
УДК 517.91
МЕТОД ОБАТНОЙ ФУНКЦИИ ДЛЯ ЕШЕНИЯ
АВТОНОМНОО ОБЫКНОВЕННОО
ДИФФЕЕНЦИАЛЬНОО УАВНЕНИЯ:
ОДНОМЕНЫЙ СЛУЧАЙ
В.В. Коробицын, Ю.В. Фролова
A novel method for numerial solving autonomous ordinary dierential
equation dx/dt = f (x) is suggested. The solution is found as a funtion t(x)
for that the inrementation is independent from previous points of solution
trajetory. This fat is used for parallel realization of omputational method.
The numerial experiments show the appliability of suggested method for
area with rapid aeleration of right-hand funtion.
Численные методы решения обыкновенных диеренциальных уравнений
появились довольно давно. Еще Ньютон и Эйлер предложили методы поиска
приближенного решения. Классические работы унге, Кутты, Адамса посвящены построению одношаговых и многошаговых методов решения ОДУ. Казалось бы, всј уже разработано и внедрено, однако интерес к разработке новых
и модиикации существующих методов не ослабевает. С одной стороны, новые
задачи требуют новых методов. Так, в последнее время активно развиваются методы решения жестких, диеренциально-алгебраических и разрывных
задач [1? [4?, [8? [11?. А с другой стороны, развитие компьютерной техники
требует разработки эективных методов решения на многоядерных, многопроцессорных и распределенных компьютерных системах. В последнее время
стали появляться работы, связанные с параллельной реализацией методов решения ОДУ [5? [7?, [12?.
Задача данной статьи заключается в представлении нового метода решения
обыкновенных диеренциальных уравнений с высокой степенью распараллеливания вычислений. ешение вычисляется независимо на отдельных промежутках по пространственной переменной x, где определяется время ? , необходимое для преодоления промежутка.
Представленный метод применим для решения одного автономного нелинейного обыкновенного диеренциального уравнения первого порядка, а также
высказана идея о распространении этого метода на автономные системы ОДУ.
c
Copyright 2009 В.В. Коробицын, Ю.В. Фролова.
Омский государственный университет.
E-mail:
korobits@rambler.ru
20 В.В.
1.
Коробицын, Ю.В. Фролова.
Метод обратной ункции решения ОДУ
Аппроксимация решения автономного ОДУ
ассмотрим задачу
dx
= f (x), x(t0 ) = x0 , t > t0 .
(1)
dt
Необходимо найти решение x(t) : R ? R задачи Коши (1), учитывая, что ункция правой части f является непрерывной по переменной x и не зависит от t в
явном виде.
ешение x(t) можно разложить в ряд Тейлора вблизи точки t0
x(t0 + ? ) = x(t0 ) +
dx
d2 x
?2
(t0 ) · ? + 2 (t0 ) ·
+ ...
dt
dt
2!
(2)
Частичную сумму ряда (2) обозначим за Pk (t0 , ? ),
Pk (t0 , ? ) :=
тогда
Введем обозначения
dx
d2 x
?2
dk x
?k
(t0 ) · ? + 2 (t0 ) ·
+ · · · + k (t0 ) · ,
dt
dt
2!
dt
k!
x(t0 + ? ) ? x(t0 ) = Pk (t0 , ? ) + O(? k+1).
f0 = f (x0 ), f0? =
df
d2 f
(x0 ), f0?? = 2 (x0 )
dx
dx
и приведем выражения для Pk , k = 0, 1, 2, . . . :
P0 (t0 , ? ) = 0,
dx
(t0 ) · ? = f0 ?,
dt
d2 x
?2
?2
P2 (t0 , ? ) = P1 (t0 , ? ) + 2 (t0 ) ·
= f0 ? + f0? f0 ,
dt
2!
2
3
2
3
?3
?
?
dx
= f0 ? + f0? f0 + f0?? f02 + f0?2 f0
,···
P3 (t0 , ? ) = P2 (t0 , ? ) + 3 (t0 ) ·
dt
3!
2
6
В отличие от традиционного подхода, когда по заданному шагу ? находят
значение решения x(t0 + ? ), мы зададим малое приращение dx по переменной x
и будем искать ? такое, что
P1 (t0 , ? ) = P0 (t0 , ? ) +
x(t0 + ? ) ? x(t0 ) = dx .
В этом случае необходимо решить алгебраическое уравнение
Pk (t0 , ? ) = dx .
(3)
Функция Pk (t0 , ? ) является многочленом от ? степени k , поэтому уранение (3)
будет иметь k корней. Необходимо найти наименьший действительный положительный корень, который будет давать значение приращения по переменной t,
необходимое для достижения точки x0 + dx .
Математические структуры и моделирование. 2009. Вып. 19.
21
Отметим, что решение уравнения (1) перейдет из точки x0 в точку x0 + dx
при dx > 0, только когда f (x) > 0 для всех x ? [x0 ; x0 + dx ]. Однако если в
некоторой точке x1 этого интервала ункция f (x1 ) = 0, то это означает, что в
этой точке уравнение (1) имеет особую точку. Поэтому решение не достигнет
точки x0 +dx . Если же f (x0 ) < 0, то изучается интервал [x0 ; x0 ?dx ] аналогичным
образом.
2.
Формулы для вычисления
?
Далее будем считать, что f (x) > 0 для всех x ? [x0 ; x0 +dx ]. Введем обозначения
f1 = f (x0 + dx ), f?1 = f (x0 ? dx ).
Тогда значения производных от ункции f в точке x0 можно оценить ормулами
f1 ? f?1
f1 ? 2f0 + f?1
f0? ? f0? (dx ) :=
(4)
.
, f0?? ? f0?? (dx ) :=
2dx
d2x
Оба приближения имеют погрешность порядка O(d3x ).
Найдем аналитические выражения для корней уравнения (3) при разных k .
При k = 1 получаем уравнение f0 ? = dx . Корень этого линейного уравнения
определяется выражением ? = dfx0 . Погрешность решения, соответствующего
найденному ? , будет порядка O(d2x).
При k = 2 получаем квадратное уравнение относительно ?
f0? f0 2
? + f0 ? ? dx = 0,
2
корни которого определяются выражениями
D = f02 + 2f0 f0? dx ,
?1,2
s
?
1
?f0 ± D
f0?
= ? ?1 ± 1 + 2dx
.
=
f0 f0?
f0
f0
Если f0? > 0, то положительный действительный корень будет
s
1
f0?
?= ?
1 + 2dx ? 1 .
f0
f0
(5)
Если f0? < 0, то для существования действительного корня необходимо удовлетворение условия 1 + 2dx f0? /f0 ? 0, что будет выполнено при dx ? ?f0 /(2f0? ).
Учитывая выражение для f0? из (4), получаем dx ? ?f0 dx /(f1 ? f?1 ). Отсюда
?f0 /(f1 ? f?1 ) ? 1. Следовательно, dx нужно выбрать таким, чтобы было выполнено условие f1 ? f?1 ? ?f0 . Тогда минимальное ? для достижения x0 + dx
соответствует (5).
При k ? 3 получаем алгебраическое полиномиальное уравнение относительно ? степени k . Это уравнение можно решать численно. При этом необходимо
взять наименьший положительный корень ? , соответствующий первому пересечению кривой полинома с точкой x + dx .
22 В.В.
Коробицын, Ю.В. Фролова.
Метод обратной ункции решения ОДУ
Однако мы предлагаем повысить точность вычисления ? , не повышая степени полинома. Для этого используем разложение ункции решения в точке
x1 = x(t1 ), t1 = t0 +? , вычисляя значение ункции решения в точке x0 = x1 ?dx ,
получаем
dx
d2 x
?2
?k
dk x
(t1 ) · ? + 2 (t1 ) ·
+ · · · + (?1)k k (t1 ) ·
+ ...
dt
dt
2!
dt
k!
x(t0 ) = x(t1 ) ?
Тогда x(t0 ) ? x(t1 ) = Pk (t1 , ?? ) + O(? k+1 ). Учитывая выражение (3), получаем два соотношения x1 ? x0 = Pk (t0 , ? ) и x0 ? x1 = Pk (t1 , ?? ). Вычитая эти
соотношения, получаем
2dx = Pk (t0 , ? ) ? Pk (t1 , ?? ).
(6)
Введем обозначения
f2 = f (x0 + 2dx ),
f1? =
df
f2 ? f0
(x1 ) ?
,
dx
2dx
f1?? =
d2 f
f2 ? 2f1 + f0
. (7)
(x1 ) ?
2
dx
d2x
Вычислим выражение для ? как корень уравнения (6) для разных значений
параметра k .
При k = 1 получаем линейное уравнение 2dx = f0 ? ? f1 · (?? ). ешение
этого уравнения ? = 2dx /(f0 + f1 ). Погрешность определения величины ? будет
порядка O(d2x).
При k = 2 получаем квадратное уравнение
f0? f0 ? f1? f1
?2
+ (f0 + f1 )? ? 2dx = 0.
2
(8)
При k ? 3 получаем полиномиальное уравнение относительно ? , наименьший положительный вещественный корень которого можно вычислять численно.
Далее мы найдем аналитическое выражение для наименьшего положительного вещественного корня уравнения (8) и оценим погрешность его вычисления
при аппроксимации производных ункции f , входящих в уравнение, конечными разностями вида (7).
Выражение для коэициента уравнения при квадрате можно вычислить
f0? f0 ? f1? f1 =
f2 ? f0
?f?1 f0 + 2f0 f1 ? f1 f2
f1 ? f?1
f0 ?
f1 =
.
2dx
2dx
2dx
Найдем корни уравнения (8)
D = f0 + f1
?1,2
2
2
? 4 f0? f0 ? f1? f1 ?2dx = f0 + f1 + 4 ?f?1 f0 + 2f0 f1 ? f1 f2 ,
?
f0 + f1 ?
? f0 + f1 ± D
= dx
=
?
?
2(f0 f0 ? f1 f1 )
q
2
? 4 f?1 f0 ? 2f0 f1 + f1 f2
.
f?1 f0 ? 2f0 f1 + f1 f2
f0 + f1
Математические структуры и моделирование. 2009. Вып. 19.
23
Введем обозначение
F =
f0 + f1
f?1 f0 ? 2f0 f1 + f1 f2
и перепишем выражение для корней уравнения
s
?1,2 = dx F 1 ±
4
.
1?
f0 + f1 F
Для существования вещественных корней уравнения необходимо выполнение
4 условия
? 1, это будет выполнено в одном из случаев (f0 + f1 )F < 0
f0 +f1 F
или (f0 + f1 )F ? 4.
Наименьший положительный вещественный корень всегда определяется выражением
s
4
.
? = dx F 1 ? 1 ?
(9)
f0 + f1 F
Причем если произведение dx F >
0, то должно быть f0 + f1 F ? 4, а если
dx F < 0, то должно быть f0 + f1 F < 0.
Каким должен быть dx ? Знак величины dx определяется знаком ункции f .
Если f0 > 0, то должно быть dx > 0. Длина шага dx будет зависеть от оценки
погрешности.
Оценим погрешность найденного корня (9). При вычислении корня производные ункции f были заменены конечными разностями (7), погрешность вычисления которых порядка O(d2x ). Следовательно, погрешность вычисления F
будет порядка O(d3x ), а погрешность определения ? порядка O(d4x).
3.
Практическая оценка погрешности
Для практической оценки погрешности вычисления значения ? проведем вычисления с двумя шагами длиной dx и одним шагом длиной 2dx , тогда
t(x0 + 2dx ) ? t0 = ?1 + ?2 + C · (dx )p ,
t(x0 + 2dx ) ? t0 = ?3 + C · (2dx )p ,
где ?1 , ?2 два значения, вычисленные при шаге dx , ?3 одно значение, соответствующее шагу 2dx . Вычитая первое выражение из второго, получим
?1 + ?2 ? ?3 = C · (dx )p · (2p ? 1).
Тогда погрешность вычисления ? с шагом dx будет
Errdx = C · (dx )p =
?1 + ?2 ? ?3
.
2p ? 1
24 В.В.
4.
Коробицын, Ю.В. Фролова.
Метод обратной ункции решения ОДУ
Численные эксперименты
Точность предложенного метода нахождения численного решения автономного
ОДУ оценим на двух примерах.
Первое уравнение линейное:
dx
= x, x(0) = x0 , t > 0.
dt
(10)
Аналитическое решение этого уравнения определяется ормулой x(t) = x0 et .
Откуда можно выразить ункцию t(x) = ln(x/x0 ).
Для численного эксперимента использовалась ормула (9). Были заданы
параметры: x0 = 1, dx = 0.19. раик решения уравнения (рис. 1-а) показывает
экспоненциальный рост x(t). Порядок относительной погрешности численного
решения уменьшается от -1.833 до -2.845 (рис. 1-б). Судя по данному граику,
можно сделать вывод, что предложенная схема вычисления решения лучше
работает в условиях больших значений ункции правой части уравнения.
а)
б)
ис. 1. а) раик решения системы (10), б) Порядок относительной погрешности
численного решения
Второе рассматриваемое уравнение модель Верхюльста-Пјрла:
dx
= k(a ? x)x, x(0) = x0 , t > 0.
dt
(11)
Аналитическим решением уравнения является логистическая кривая, определяемая ормулами:
x(t) =
aeka(t?C1 )
при 0 < x < a,
1 + eka(t?C1 )
x(t) =
aeka(t?C2 )
при x > a.
1 ? eka(t?C2 )
Значения констант определяются исходя из начальных данных
C1 =
x0
? ln a?x
0
ka
, C2 =
0
? ln x0x?a
ka
.
25
Математические структуры и моделирование. 2009. Вып. 19.
Численное решение построено при заданных параметрах: a = 3, k = 2, dx = 0.04.
Начальные данные для построения двух траекторий решения: x10 = 1, x20 = 5.
Первая траектория (рис. 2-а) помечена квадратиками, вторая кружками.
а)
б)
ис. 2. а) раик решения системы (11), б) Порядок относительной погрешности
численного решения; линии с квадратиками для
0 < x < a,
линии с кружками для
x > a.
Уравнение (11) имеет особую точку x = a, при приближении к которой
рост времени t приводит к малым изменениям переменной x. Поэтому в нашем
случае, когда шаг dx задан постоянным и вычисляется шаг по времени ? , при
приближении к особой точке происходит интенсивный рост величины ? , что
приводит к повышению погрешности решения. Это можно увидеть на обеих
кривых порядка относительной погрешности решения (рис. 2-б).
Полученные результаты вычислений показывают рациональность использования предложенного метода для поиска решения в области интенсивного роста
значений ункции правой части уравнения. Поэтому разумно предположить
эективность использования этого метода для решения жестких уравнений.
лавным преимуществом представленного метода является его параллельность на каждом i-м интервале вычисление шагаP?i можно проводить параллельно, после чего решение можно составить x(t + ni=1 ?i ) = x0 + ndx . А также
для нахождения одного шага требуется только одно новое вычисление ункции
правой части, что по вычислительным затратам сравнимо с методом Эйлера,
но точность предлагаемого метода значительно выше.
5.
Идея об определении траектории решения автономной
системы ОДУ
ассмотрим автономную систему обыкновенных диеренциальных уравнений
первого порядка размерности 2:
? dx
1
= f (x1 , x2 ),
?
dt
?
? dx
2
= g(x1 , x2 ),
dt
(12)
x1 (t0 ) = x01 ,
?
?
?
x2 (t0 ) = x02 .
26 В.В.
Коробицын, Ю.В. Фролова.
Метод обратной ункции решения ОДУ
Применяя метод разделения переменных к каждому уравнению в отдельности,
получим
Z
Z
Z
dx1
dx1
= dt ? t1 (x1 , x2 ) =
+ C1 ,
f (x1 , x2 )
f (x1 , x2 )
из второго уравнения
Z
dx2
+ C2 .
t2 (x1 , x2 ) =
g(x1 , x2 )
Введенные здесь ункции t1 (x1 , x2 ) и t2 (x1 , x2 ) описывают поведение интегралов
от обратных ункций правых частей системы (12). Константы C1 , C2 определяются исходя из начальных условий системы.
Предполагаем, что решение системы (12) можно определить как параметрическую кривую ? в пространстве R2 вида
? = {(x1 , x2 ) ? R2 | t1 (x1 , x2 ) ? t2 (x1 , x2 ) = 0}.
Параметр s кривой ? определяется s = t1 (x1 , x2 ) = t2 (x1 , x2 ), (x1 , x2 ) ? ?.
Для нахождения численного решения системы (12) необходимо разбить азовое пространство системы на квадратные участки Sij , (i ? Z, j ? Z) со стороной dx :
Sij = {(x1 , x2 ) ? R2 | x01 + idx ? x1 ? x01 + (i + 1)dx , x02 + jdx ? x2 ? x02 + (j + 1)dx }.
Для каждого квадрата Sij вычисляются значения ункций правых частей системы в углах этого квадрата. Используя билинейную интерполяцию, вычисляются интегралы для определения t1 , t2 . Зная точку входа траектории решения
в квадрат, определяются константы C1 и C2 , после чего находится точка выхода
траектории решения. Эта точка, в свою очередь, является точкой входа в соседний квадрат, где процедура повторяется. Таким образом, траектория решения
определяется точками, на границах квадратов.
Отметим, что вычисление интегралов в каждом квадрате Sij является независимым от значений ункций в других квадратах, что дает возможность производить вычисления в параллельном режиме. Определение конкретной траектории решения, соответствующей заданной начальной точке, будет основано на
определении точек выхода для каждой клетки, для чего потребуется находить
только константы C1 и C2 по заранее вычисленным интегралам. Причем эти
интегралы не нужно вычислять повторно, если необходимо найти другую траекторию, достаточно вычислить константы. Что дает ускорение в определении
азавого портрета системы.
Для увеличения точности определения точки пересечения границы необходимо заменить билинейную интерполяцию на биквадратичную, бикубическую
и т.д.
Предложенную идею определения траектории решения систему ОДУ можно распространить на n-мерный случай. Траектория будет определяться как
пересечение n ? 1 гиперповерхностей размерности n ? 1.
27
Литература
1. Деккер, К. Устойчивость методов унге-Кутты для жестких нелинейных диеренциальных уравнений / К. Деккер, Я. Вервер. М.: Мир, 1988.
2. Коробицын, В.В. Алгоритм численного решения диеренциальных уравнений с
разрывной правой частью / В.В. Коробицын, Ю.В. Фролова // Математические
структуры и моделирование. 2005. Вып. 15. С.4654.
3. Коробицын, В.В. Исследование поведения явных методов унге-Кутты при решении систем обыкновенных диеренциальных уравнений с разрывной правой частью / В.В. Коробицын, В.Б. Маренич, Ю.В. Фролова // Математические структуры и моделирование. 2007. Вып. 17. С.1925.
4. Коробицын, В.В., Алгоритм численного решения кусочно-сшитых систем /
В.В. Коробицын, Ю.В. Фролова, В.Б. Маренич // Вычисл. технологии. 2008.
Т.13, N.2. С.7081.
5. Новиков, Е.А. еализация полуявного (4,2)-метода решения жестких систем на
параллельной ЭВМ / Е.А. Новиков, Л.П. Каменщиков // Вычисл. технологии. 2004. Т.9. Вестник КазНУ им. аль-Фараби. Сер. Математика, механика, инорматика. N.3 (42). (Совм. выпуск. Ч.1). С.235241.
6. Фельдман, Л.П. Сходимость и оценка погрешности параллельных одношаговых
блочных методов моделирования динамических систем с сосредоточенными параметрами / Л.П. Фельдман // Научные труды ДонНТУ. Серия: Инорматика,
моделирование и вычислительные методы, (ИКВТ-2000), выпуск 15. Донецк:
ДонНТУ, 2000. С.3439.
7. Фельдман, Л.П. Параллельные алгоритмы численного решения задачи Коши для
систем обыкновенных диеренциальных уравнений / Л.П. Фельдман, И.А. Назарова // Матем. моделирование. 2006. Т.18, N.9. С.1731.
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.201222.
11. Gear, C.W. Solving ordinary dierential equations with disontinuities / C.W. Gear,
O. Шsterby // ACM Trans. Math. Software. 1984. V.10. P.2344.
12. Jakson, K. The potential for parallelism in Runge-Kutta methods. Part I: RK formulas
in standard form / K. Jakson, S. Norsett // SIAM J. Numer. Anal. 1995. V.32. P.4982.
Документ
Категория
Без категории
Просмотров
8
Размер файла
218 Кб
Теги
обыкновенное, автономное, решение, уравнения, дифференциальной, обратное, метод, функции, одномерных, случай
1/--страниц
Пожаловаться на содержимое документа