(i)
∂x
t=t
0
˙x
0
−
∂x
(i)
∂ ˙x
¨
x
t=t
0
.
Построим последовательность полиномов (2), приближающих ре-
шение (4):
x
(n)
(t, x
0
, ˙x
0
) = c
(n)
0
a
0
+ c
(n)
1
a
1
t + c
(n)
2
a
2
t
2
+ . . . + c
(n)
m
n
a
m
n
t
m
n
,
где c
(n)
m
n
– так называемые множители сходимости, a
i
– коэффициен-
ты Тейлора из (5).
Далее для численного решения поставленной задачи необхо-
димо согласовать начальные данные Коши и параметры ω и µ
уравнения (3). Уравнение Дюффинга, помимо периодических реше-
ний, допускает и решения апериодические.
Воспользуемся тем, что общий интеграл – интеграл энергии урав-
нения (3) выписывается в эллиптических функциях [5]
x = Csn (σt, k),
(6)
где C – амплитуда, σ – частота, k – модуль эллиптической функции,
0 ≤ k ≤ 1, t – время.
Параметры ω и µ в уравнении (3) могут быть определены через
C, σ и k следующим образом
σ
2
(1 + k
2
) = ω
2
,
(7)
2σ
2
k
2
C
2
= µ.
(8)
В [5] показано, что общий интеграл уравнения (3) имеет форму
(6), если только параметры C, σ и k удовлетворяют соотношениям
(7) и (8).
Характер колебаний определяется только начальной энергией.
Если начальная энергия меньше некоторого критического значения,
40
то процесс, описываемый уравнением (3), носит колебательный ха-
рактер. Значит достаточно подобрать такие начальные данные и па-
раметры ω, µ, чтобы решения были периодическим. Все необходимые
формулы для этого выписаны.
В качестве среды для вычисления использовалась система ком-
пьютерной алгебры MAPLE. Были проведены расчеты на компью-
тере AMD Athlon XP 1700+, 768MB of RAM. Коэффициенты схо-
димости были рассчитаны заранее и записаны в отдельных файлах
на компьютере. Вызвал серьезные трудности расчет коэффициентов
Тейлора. Хотелось получить их в символьном виде, но из-за пере-
грузки оперативной памяти компьютера и, в соответствии с этим,
очень большим временем счета программы пришлось от этого от-
казаться, и сохранять в памяти уже посчитанные коэффициенты
Тейлора. Начальные данные для уравнения Дюффинга были тоже
рассчитаны программным путем.
При t
0
= 0, x
0
= 0, ˙x
0
= 0, 1, ω = 1, 005024997, µ = 2, 010049994
в таблицах приведены рассчитанные оценки |x − x
(n)
| в зависимости
от n и ν.
Таблица 1. n = 10, ν = 1
Номер полинома
Кол-во членов
Оценка
1
1
0,25324998∗10
5
2
3
0,25303955∗10
5
3
7
0,25268883∗10
5
4
15
0,25246086∗10
5
5
31
0,25230724∗10
5
6
63
0,25219782∗10
5
7
127
0,25211622∗10
5
8
255
0,25205315∗10
5
9
511
0,25079821∗10
5
10
1023
0,24934176∗10
5
Таблица 2. n = 6, ν = 2
Номер полинома
Кол-во членов
Оценка
1
2
0,25324998∗10
5
2
8
0,25198739∗10
5
3
26
0,25175357∗10
5
4
80
0,25167175∗10
5
5
242
0,25163387∗10
5
6
728
0,25094359∗10
5
41
Литература
1. Маркушевич А.И. Теория аналитических функций. М.: Наука,
1968.
2. Голубев В.В. Лекции по аналитической теории дифференциаль-
ных уравнений. Л.: ГИТТЛ. 1950.
3. Сансоне Дж. Обыкновенные дифференциальные уравнения. М.,
1953.
4. Брумберг В.А. Ряды полиномов в задаче трех тел. Бюлл. ИТА, 9,
4 (107). 1963. С. 234–256.
5. Моисеев Н.Н. Асимптотические методы нелинейной механики. М.:
Наука, 1981.
6. Иванова О.А. Решение задачи Коши в звезде Миттаг–Леффлера
для аналитических дифференциальных уравнений // Процессы
управления и устойчивость: Труды 36-й научной конференции ас-
пирантов и студентов / Под ред. Н.В. Смирнова, В.Н. Старкова.
– СПб.: Изд-во СПбГУ, 2005. С. 53–59.
42
Калинина Е.А., Самарина О.Н.
Санкт-Петербургский государственный университет
Минимизация полной погрешности метода Эйлера
для уравнения ˙x = ax при расчетах с точностью
реальной процессорной арифметики
1. Постановка задачи. Рассмотрим задачу Коши для обыкно-
венного дифференциального уравнения
˙x = ax,
x(0) = 1.
(1)
Найдем значение решения этого уравнения при t = 1 с помо-
щью метода Эйлера. Разобьем отрезок [0,1] на n равных частей. При
данном числе шагов n получим следующее приближенное значение
решения (вычислительные формулы см. в [1]): x(1) = 1 +
a
n
n
.
Погрешность метода Эйлера стремится к нулю с ростом числа
шагов n. Однако, с увеличением n на полученный результат б´ольшее
влияние оказывает вычислительная погрешность (ошибки округле-
ния) [3, 4].
Задача. Найти число шагов n для задачи (1), при котором пол-
ная погрешность (погрешность метода в сумме с вычислительной
погрешностью) будет наименьшей.
Замечание 1. Везде в дальнейшем будем считать n натураль-
ным, x — вещественным.
2. Решение задачи. Пусть значение a задано с относительной
погрешностью δa. Найдем оптимальное значение n (в смысле наи-
меньшей относительной погрешности) для вычисления e
a
по форму-
ле f (n) = 1 +
a
n
n
.
Лемма. При n>|a|>0 справедливо неравенство
e
a
(1+
a
n
)
n
−1>0.
Доказательство. Найдем максимальное значение функции
g(x) = ln(1+x)−x при |x| < 1. Так как g (x) = −x/(x+1), то функция
g(x) достигает максимального значения при x = 0, g(0) = 0. Таким
образом, g(x) ≤ 0 при |x| < 1. При x = a/n (n > |a| > 0) получаем
неравенство ln(1 +
a
n
) <
a
n
=⇒ (1 +
a
n
)
n
< e
a
.
С учетом леммы полная относительная погрешность вычисления
e
a
по формуле f (n) = 1 +
a
n
n
находится следующим образом [2]:
43
δ(n) =
e
a
1 +
a
n
n
+ (ε + δa)n − 1,
где ε — точность вычислений (так, для f loat (4 байта) ε ≈ 1, 19·10
−7
,
для double (8 байт) ε ≈ 2, 22 · 10
−16
, для long double (10 или 12 байт
в зависимости от системы) ε ≈ 1, 08 · 10
−19
).
Дифференцируя функцию δ(x) по x, получаем необходимое усло-
вие экстремума
δ (x) = −
e
a
ln 1 +
a
x
−
a
x
(
1+
a
x
)
1 +
a
x
x
+ (ε + δa) = 0.
Утверждение. Функция δ (x) при x > |a| имеет единственный
вещественный корень, который является точкой минимума функ-
ции δ. Этот корень приближенно находится по формуле
x ≈
|a|
2(ε + δa)
(2)
.
Доказательство. Докажем сначала единственность корня функ-
ции δ (x).
При x > |a| рассмотрим функцию
h(x) = ln 1 +
a
x
−
a
x 1 +
a
x
.
Поскольку ее производная h (x) = −
a
2
x(x + a)
< 0 на рассматривае-
мом промежутке, то h строго убывает. Кроме того, lim
x→+∞
h(x) = 0.
Таким образом, h(x) > 0 при x > |a|.
Кроме того, при x > |a| для функции u(x) = −
e
a
1 +
a
x
x
имеем
u (x) =
e
a
1 +
a
x
x
ln 1 +
a
x
−
a
x 1 +
a
x
=
h(x)e
a
1 +
a
x
x
> 0.
Так как функция δ (x) − (ε + δa) является произведением поло-
жительной строго убывающей функции h и отрицательной строго
возрастающей функции u, то она отрицательная строго возрастаю-
щая функция. Кроме того,
44
lim
x→+∞
(δ (x) − (ε + δa)) = 0.
Значит, δ (x) имеет ровно один вещественный корень при x > |a|,
и при переходе через этот корень (при возрастании x) меняет знак с
отрицательного на положительный. Тем самым, корень производной
является точкой минимума функции δ.
Теперь покажем, что этот корень можно найти по формуле (2).
В самом деле, при x → +∞ имеем
ln 1 +
a
x
−
a
x 1 +
a
x
=
a
2
2x
2
−
2a
3
3x
3
+ o
a
3
x
3
,
и при x >> a также
e
a
(
1+
a
x
)
= 1 + o
a
x
.
Тогда δ (x) = − 1 + o
a
x
a
2
2x
2
+ o
a
2
x
2
+ (ε + δa).
Таким образом, для нахождения точки минимума функции δ(x)
получаем уравнение −
a
2
2x
2
+ o
a
2
x
2
+ (ε + δa) = 0, и при x >> |a|
имеем (2). Действительно, полученное значение x удовлетворяет
условию x >> |a|.
Таким образом, искомое значение n равно
n
1
=
|a|
2(ε + δa)
или n
2
=
|a|
2(ε + δa)
+ 1.
(3)
Здесь x обозначает целую часть x.
Замечание 2. Для минимизации полной относительной погреш-
ности вычислять оптимальное значение n по формуле (3) следует с
максимальной возможной точностью. В приведенных далее приме-
рах указана точность вычислений.
3. Численные примеры. Приведем несколько примеров вычис-
ления полной погрешности метода Эйлера для задачи Коши (1).
Пример 1. Найдем оптимальное n при a = 1 для δa = 0,
ε ≈ 1, 19 · 10
−7
. Получаем n
1
= 2049, n
2
= 2050 (точные значения).
45
Таблица 1. Значения относительных погрешностей δ(n) (float, C)
для различных значений n
n
δ(n)
n
δ(n)
n
δ(n)
2043
0,000250
2047
0,000243
2051
0,000243
2044
0,000247
2048
0,000247
2052
0,000247
2045
0,000247
2049
0,000243
2053
0,000250
2046
0,000243
2050
0,000243
2054
0,000250
Пример 2. Найдем оптимальное n при a = 40 для δa = ε ≈
≈ 1, 19 · 10
−7
. Получаем n
1
= 57926, n
2
= 57927 (вычислено с точно-
стью float).
Таблица 2. Значения относительных погрешностей δ(n) (float, Visual C++)
для различных значений n
n
δ(n)
n
δ(n)
n
δ(n)
57921
0,0286704
57925
0,0258728
57929
0,0300713
57922
0,0279703
57926
0,0251746
57930
0,0293703
57923
0,0272706
57927
0,0244769
57931
0,0286699
57924
0,0265715
57928
0,0307727
57932
0,0279698
Пример 3. Найдем оптимальное n при a = −50 для δa = ε ≈
≈ 1, 19 · 10
−7
. Получаем n
1
= 72407, n
2
= 72408 (вычислено с точно-
стью float).
Таблица 3. Значения относительных погрешностей δ(n) (float, Visual C++)
для различных значений n
n
δ(n)
n
δ(n)
n
δ(n)
72403
0,034704
72407
0,0331219
72411
0,0359336
72404
0,0354073
72408
0,0338241
72412
0,0366378
72405
0,0361112
72409
0,0345268
72413
0,0359553
72406
0,0368155
72410
0,03523
72414
0,0336573
Литература
1. Бахвалов Н.С. Численные методы. М.: Наука, 1973. 631 с.
2. Демидович Б.П., Марон И.А. Основы вычислительной математи-
ки. М.: Физматгиз, 1963. 660 с.
3. Ортега Дж., Пул У. Введение в численные методы решения диф-
ференциальных уравнений. М.: Наука, 1986. 288 c.
4. Higham N.J. Accuracy and stability of numerical algorithms. SIAM,
1996. 668 p.
46
Купцов С.Ю.
Санкт-Петербургский государственный университет
Об устойчивости однопараметрических семейств
систем линейных дифференциальных уравнений
Рассмотрим линейную стационарную систему
A(ε)x
(s)
(t) +
s−1
p=0
m
q=1
B
pq
(ε)x
(p)
t − h
q
(ε) = 0,
(1)
где x(t) – неизвестная вектор–функция из R
n
, A и B
pq
–
(n × n)-матрицы, h
q
– запаздывания системы, ε – скалярный па-
раметр. Будем считать, что ε меняется в пределах области E =
(−∞, +∞), и на этом множестве элементы всех матриц системы, а
также её запаздывания, задаются непрерывными функциями. Кроме
того, для всех ε ∈ E выполнено det A(ε) = 0, h
q
(ε) ≥ 0. Предположим
также, что известны оценки
A(ε + σ) − A(ε) ≤ a(σ, ε),
|h
q
(ε + σ) − h
q
(ε)| ≤ H(σ, ε),
B
pq
(ε + σ) − B
pq
(ε) ≤ b(σ, ε),
B
pq
(ε) ≤ b
0
(ε).
Здесь и далее в качестве нормы произвольной матрицы выбирает-
ся максимум из норм её столбцов, а норма столбца – эрмитова.
Причём, не умаляя общности, будем считать, что a(σ, ε), b(σ, ε) и
H(σ, ε) — непрерывные функции, заданные на множестве ε ∈ E,
σ ∈ (−∞, +∞), строго монотонно убывающие от +∞ до 0 при изме-
нении σ от −∞ до 0 и строго монотонно возрастающие от 0 до +∞
при изменении σ от 0 до +∞.
Пусть система (1) асимптотически устойчива по Ляпунову при
ε = ε
0
∈ E. Целью данной работы является построение последова-
тельностей ε
k
и ε
k
, обладающих свойствами:
1. ε
k
≤ ε
0
, и данная последовательность чисел монотонно убыва-
ет; ε
k
≥ ε
0
, и данная последовательность монотонно возраста-
ет;
2. Любая система, входящая в семейство (1) и отвечающая зна-
чениям ε ∈ [ε
k
, ε
k
], асимптотически устойчива по Ляпунову;
47
3. Либо последовательность ε
k
→ −∞, либо последовательность
ε
k
→ ε, и тогда система (1), отвечающая ε = ε, не является
асимптотически устойчивой по Ляпунову;
4. Либо последовательность ε
k
→ +∞, либо последовательность
ε
k
→ ε, и тогда система (1), отвечающая ε = ε, не является
асимптотически устойчивой по Ляпунову.
Характеристическая матрица системы (1) имеет вид
G(λ, ε) = A(ε)λ
s
+
s−1
p=0
m
q=1
B
pq
(ε)λ
p
e
−h
q
(ε)λ
.
Причем, асимптотическая устойчивость системы (1) по Ляпунову
равносильна (см., например, [2], стр. 21, 32) отсутствию корней урав-
нения g(λ, ε) = det G(λ, ε) = 0 на множестве Re λ ≥ 0.
Введём в рассмотрение функцию
Φ(A) =
n
| det A| + A
n
− A
и предложим следующий алгоритм построения ε
k
и ε
k
:
1. Полагаем k = 0 и фиксируем некоторое α ∈ (0, 1);
2. Находим R
k
= 1 + m max
σ∈[0,1]
b
0
(ε
k
± σ)
Φ(A(ε
k
± σ))
;
3. Определяем ∆
k
=
min
ω∈[−R
k
,R
k
]
Φ G(iω, ε
k
) ;
4. Для Θ
k
(σ) = a(σ, ε
k
)R
s
k
+ m
R
s
k
− 1
R
k
− 1
b(σ, ε
k
) + b
0
(ε
k
)R
k
H(σ, ε
k
)
находим σ
k
> 0 из условия Θ
k
(σ
k
) ∈ [α∆
k
, ∆
k
);
5. Полагаем ε
k+1
= ε
k
± min {σ
k
, 1} и переходим ко второму пунк-
ту алгоритма, увеличив значение k на 1.
Если при прохождении этого алгоритма во втором и пятом пунк-
тах всегда выбирать верхний (нижний) знак, то построенная после-
довательность ε
k
, k = 0, 1, . . . , фактически и будет последовательно-
стью ε
k
(соответственно, ε
k
).
Покажем это на примере последовательности ε
k
, имея в виду,
что для ε
k
доказательство будет аналогичным. Все рассуждения бу-
дут базироваться на утверждении из [2], согласно которому, для лю-
бых матриц A и B с комплексными элементами выполнение условия
B < Φ(A) влечёт за собой неособость A + B.
48
Прежде всего заметим, что выбор R
k
> 1 во втором пункте ал-
горитма гарантирует отсутствие корней уравнения g(λ, ε
k
+ σ) = 0
на множестве Re λ ≥ 0, |λ| ≥ R
k
для всех σ ∈ [0, 1]. Данное об-
стоятельство представляется очевидным после представления мат-
рицы G(λ, ε
k
+ σ) в виде λ
s
(A(ε
k
+ σ) + Y ) и оценки нормы Y на
указанном выше множестве. Вместе с этим, выбор σ
k
из условия
Θ
k
(σ
k
) < ∆
k
и свойство монотонности Θ
k
(σ) гарантируют отсут-
ствие при всех σ ∈ [0, σ
k
] корней уравнения g(iω, ε
k
+ σ) = 0 в пре-
делах отрезка ω ∈ [−R
k
, R
k
]. В этом легко убедиться, если предста-
вить матрицу G(iω, ε
k
+σ) в виде G(iω, ε
k
)+Z и оценить Z сверху,
воспользовавшись неравенством |e
iϕ
− 1| = 2| sin
1
2
ϕ| ≤ |ϕ|, справед-
ливом для всех действительных ϕ. В итоге, выполнение неравенств
0 ≤ σ ≤ min {σ
k
, 1} приводит к отсутствию чисто мнимых корней
уравнения g(λ, ε
k
+ σ) = 0, а это, в свою очередь, к совпадению
количества нулей на множестве Re λ ≥ 0 у функций g(λ, ε
k
+ σ) и
g(λ, ε
k
). Другими словами, для всех ε ∈ [ε
0
, ε
1
]∪[ε
1
, ε
2
]∪. . .∪[ε
k−1
, ε
k
]
количество нулей функции g(λ, ε) на множестве Re λ ≥ 0 остаётся
постоянным, и, следовательно, равно нулю, как у функции g(λ, ε
0
).
Наконец, если последовательность ε
k
сходится к ε, то очевидно,
что, во–первых, равенство min {σ
k
, 1} = 1 может реализовывать-
ся только конечное число раз, а значит σ
k
→ +0, и, во–вторых,
последовательность R
k
ограничена сверху некоторым числом R.
Последнее приводит к тому, что α∆
k
≤ Θ
k
(σ
k
) → +0. Однако,
G(iω, ε
k
) равномерно ограничена на множестве ω ∈ [−R, R], а то-
гда условие ∆
k
→ +0 равносильно существованию последовательно-
сти ω
k
∈ [−R, R], для которой выполнено det G(iω
k
, ε
k
) → 0. Выби-
рая из этой последовательности сходящуюся подпоследовательность,
приходим к наличию чисто мнимого корня у уравнения g(λ, ε) = 0,
что полностью завершает доказательство.
Литература
1. Зубов В.И. Лекции по теории управления. М.: Наука, 1975. 496 c.
2. Купцов С.Ю. Оценка радиуса устойчивости линейной стационар-
ной системы дифференциально–разностных уравнений // Про-
цессы управления и устойчивость: Труды 30-й научной конферен-
ции / Под ред. В.Н. Старкова. СПб.: НИИ Химии СПбГУ, 1999.
С. 103–107.
49
Купцова С.Е.
Санкт-Петербургский государственный университет
Асимптотически инвариантные множества
Рекомендовано к публикации профессором Жабко А.П.
Рассмотрим систему дифференциальных уравнений
˙x = f (t, x), Достарыңызбен бөлісу: |