close

Вход

Забыли?

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

?

Минимизация полной погрешности метода Эйлера для систем линейных дифференциальных уравнений с постоянными коэффициентами.

код для вставкиСкачать
Сер. 10. 2009. Вып. 4
ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА
УДК 519.622.2
Е. А. Калинина, О. Н. Самарина
МИНИМИЗАЦИЯ ПОЛНОЙ ПОГРЕШНОСТИ МЕТОДА ЭЙЛЕРА
ДЛЯ СИСТЕМ ЛИНЕЙНЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
С ПОСТОЯННЫМИ КОЭФФИЦИЕНТАМИ
1. Введение. На практике часто приходится иметь дело с системами линейных
дифференциальных уравнений с постоянными коэффициентами или с системами, в которых коэффициенты можно считать постоянными в течение довольно больших промежутков изменения переменной. В действительности при построении различных систем прогнозирования, АСУ, тренажеров и подобных систем разрабатываемые программные комплексы просто решают задачу интегрирования этой системы с постоянным шагом интегрирования. Иначе говоря, для решения задачи Коши применяется
метод Эйлера как наиболее простой и наименее трудоемкий. Как известно, при вычислении значения решения в точке с помощью метода Эйлера погрешность получается
довольно большой. Однако существуют только оценки погрешности метода (см., например, [1, 2]), нет достаточно хороших оценок полной погрешности – суммы погрешности метода и вычислительной погрешности. Но при уменьшении шага интегрирования,
что ведет к снижению погрешности метода, вычислительная погрешность становится
весьма существенной и не принимать ее во внимание нельзя. Можно было пренебрегать вычислительной погрешностью, когда уровень развития вычислительной техники
не позволял выполнять столь большое число операций, что влияние вычислительной
погрешности на результат становилось значительным. В современных условиях, когда компьютер выполняет вычисления в режиме реального времени, нужно учитывать
и ошибки округления. В связи с этим появился вероятностный подход к решению задачи оценки полной погрешности [3]. К сожалению, он не всегда применим, поскольку
ошибки округления нельзя считать независимыми случайными величинами [4]. Также
применялся метод поиска точного решения системы [5]. Другие подходы не использовались в связи с большой сложностью задачи. На это указывают многие авторы [4, 6–8].
В настоящей работе с помощью оценки верхней границы относительной погрешности
явно находится оптимальное значение шага интегрирования, обеспечивающее наименьшую полную погрешность. Показывается, что для некоторого круга систем линейных
дифференциальных уравнений интегрирование с помощью метода Эйлера допустимо,
оно обеспечивает нахождение значения решения в точке с достаточно малой погрешностью. Приводятся численные примеры.
Калинина Елизавета Александровна – кандидат физико-математических наук, доцент кафедры высшей математики факультета прикладной математики–процессов управления СанктПетербургского государственного университета. Количество опубликованных работ: 11. Научные направления: теория управления, устойчивость, компьютерная алгебра, символьные (аналитические) методы. E-mail: ekalinina69@gmail.com.
Самарина Ольга Николаевна – аспирант факультета прикладной математики–процессов управления. Научный руководитель: доктор физико-математических наук, проф. А. Ю. Утешов. Количество
опубликованных работ: 2. Научные направления: компьютерная алгебра, символьные (алгебраические)
методы. E-mail: samarina_o@mail.ru.
c Е. А. Калинина, О. Н. Самарина, 2009
95
2. Основной результат. Для решения задачи будем использовать следующую норму вектора:
⎞
⎛
x1
⎟
⎜
Xm×1 = ⎝ ... ⎠ : ||X|| = |x1 | + |x2 | + · · · + |xm |
xm
и соответствующую матричную норму
||A|| = max ||Aj ||,
1jm
где A = (A1 , A2 , . . . , Am ) – матрица порядка m.
Рассмотрим задачу Коши для системы линейных дифференциальных уравнений
Ẋ = AX,
X(t0 ) = X0
(X0 = 0),
(1)
здесь A = [aij ]m
i,j=1 – вещественная квадратная матрица порядка m с постоянными
элементами,
⎞
⎛
⎛ 0 ⎞
x1 (t)
x1
⎟
⎜
⎜
. ⎟
.
..
X(t) = ⎝
⎠ , X0 = ⎝ .. ⎠ .
xm (t)
x0m
Будем считать, что все элементы матрицы A и вектора X0 заданы точно.
Известно, что решение задачи (1) имеет вид X = eA(t−t0 ) X0 . Найдем полную относительную погрешность приближенного решения задачи Коши методом Эйлера в точке
t1 = t0 + 1.
При данном числе шагов n получим следующее приближенное значение решения:
n
1
X̃(t1 ) ≈ E + A X0 .
(2)
n
Пусть AJ – жорданова нормальная форма матрицы A, а T = [tij ]m
i,j=1 – матрица, со−1
ставленная из векторов канонического базиса, такая, что AJ = T AT .
Тогда равенство (2) можно записать в виде
n
1
T −1 X0 .
X̃(t1 ) ≈ T E + AJ
n
Следует учесть, что вектор X̃(t1 ) может быть нулевым только при достаточно больших значениях t1 , и в этом случае можно говорить о стремлении решения задачи Коши (1) к нулю с ростом t. Если же в нуль обращаются не все компоненты вектора
X̃(t1 ), то можно попытаться избежать этой ситуации изменением числа шагов n метода Эйлера, например увеличить или уменьшить n на 1. Обращение в нуль части
компонент вектора X̃(t1 ) возможно лишь в исключительных случаях, поскольку наименьшие модули чисел, отличных от нуля, имеют следующие значения: для вычислений
с одинарной точностью float ≈ 1.17 · 10−38, для двойной точности double ≈ 2.225 · 10−308,
для long double ≈ 3.36 · 10−4932 .
Далее рассмотрим общий случай: ни одна из компонент вектора X̃(t1 ) не обращается
в нуль.
Оптимальное (в смысле наименьшей полной погрешности) значение числа шагов
метода Эйлера можно найти, воспользовавшись следующей теоремой.
96
Теорема 1. Оптимальное число шагов n для нахождения значения решения системы дифференциальных уравнений (1) при t1 = t0 + 1 приближенно может быть
найдено по формуле
3
4
ẍm (t1 ) 4 ẍ1 (t1 ) ẍ2 (t1 ) 4
+
+
·
·
·
+
xm (t1 ) 5 x1 (t1 ) x2 (t1 ) n≈
.
(3)
2mε
Замечание 1. Здесь ε = 2υ, где υ («unit roundoff») – максимальная относительная
погрешность вычислений (так, для float (4 байта) ε ≈ 1.19 · 10−7, для double (8 байт) ε ≈
2.22 ·10−16 , для long double (10 или 12 байт в зависимости от системы) ε ≈ 1.08 ·10−19 )∗) .
Обозначим через Y (t) вектор вторых производных:
⎞ ⎛
⎞
⎛
ẍ1 (t)
y1 (t)
Y (t) = ⎝ . . . ⎠ = ⎝ . . . ⎠ = A2 X(t).
ym (t)
ẍm (t)
Тогда справедливо следующее утверждение.
Следствие 1. Оптимальное число шагов n для нахождения значения решения
системы дифференциальных уравнений (1) при t1 = t0 + 1 приближенно может быть
найдено по формуле
3
⎛
⎞
4
4
x1 (t1 )
4
⎜ x2 (t1 ) ⎟
signym (t1 )
signy1 (t1 ) signy2 (t1 )
4 1
⎟
,
,...,
n≈4
(4)
A2 ⎜
⎝ . . . ⎠.
5 2mε | x1 (t1 ) | | x2 (t1 ) |
| xm (t1 ) |
xm (t1 )
Здесь signy – знак числа y:
⎧
⎨
signy =
1, если y > 0,
0, если y = 0,
⎩
−1, если y < 0.
Замечание 2. При нахождении значения решения задачи Коши (1) в точке t1 =
t0 + τ методом Эйлера имеем
n
Aτ
X̃(t1 ) = E +
X0 .
n
Таким образом, этот случай сводится к рассмотренному случаю τ = 1 заменой матрицы A на матрицу Aτ . Следовательно, достаточно рассмотреть только случай τ = 1.
3. Доказательство теоремы 1. Рассмотрим сначала случай различных собственных чисел матрицы A.
Пусть A – квадратная матрица порядка m – имеет l пар комплексно-сопряженных
собственных чисел μ1 , μ̄1 , . . . , μl , μ̄l и m − 2l вещественных собственных числа
λ1 , λ2 , . . . , λm−2l . Для комплексных собственных чисел будем также использовать
тригонометрическую форму μj = ρj (cos ζj + i sin ζj ), μ̄j = ρj (cos(−ζj ) + i sin(−ζj )),
j = 1, 2, . . . , l, причем будем считать sin ζj > 0.
∗) Значения ε взяты из стандартного включаемого файла float.h для C-компилятора на архитектуре x86.
97
Тогда j-тая компонента вектора X(t1 ) – решения задачи Коши (1) будет иметь вид
(j)
(j)
(j)
(j)
(j)
(j)
(j)
c1 eλ1 + c2 eλ2 + · · · + cm−2l eλm−2l + d1 eμ1 + d¯1 eμ̄1 + · · · + dl eμl + d¯l eμ̄l .
Соответственно j-тую компоненту вектора X̃(t1 ) – приближенного решения задачи
Коши методом Эйлера можно записать так:
n
n
λ1
λm−2l
(j)
(j)
c1
+ · · · + cm−2l 1 +
+
1+
n
n
"
"
μ1 #n ¯(j) "
μ̄1 #n
μl #n ¯(j) "
μ̄l #n
(j)
(j)
1+
1+
+ d1 1 +
+ d1 1 +
+ · · · + dl
+ dl
.
(5)
n
n
n
n
Обозначим относительную погрешность вычисления j-й компоненты вектора X̃(t1 )
по формуле (5) через δ1 (n).
Лемма 1. δ1 (n) = nε.
Д о к а з а т е л ь с т в о. Действительно, поскольку элементы матрицы A заданы
точно, вычислительная погрешность складывается из относительных погрешностей одного деления, одного сложения и n умножений. Погрешность каждого действия равна
υ. Относительные погрешности всех слагаемых одинаковы.
Таким образом, нужно найти значение n, при котором достигается минимум функции
(j)
(j) μ1
(j) μ̄1
(j) μl
(j) μ̄l
m (j) λ1
λm−2l
¯
¯
+ d1 e +d1 e + . . . +dl e +dl e
c1 e + . . . +cm−2l e
−1
F (n)=
+ εmn,
m−2l
l
n
"
"
#
#
n
λk
μk n ¯(j)
μ̄k
j=1 (j)
(j)
1+
dk 1 +
ck
+
+ dk 1 +
n
n
n
k=1
k=1
или
m−2l
l
#
"
(j)
(j)
(j)
λ
ρ
cos
ζ
k
k
k
c
e
+
2r
e
cos
ρ
sin
ζ
+
χ
k
k
k
k
k
m k=1
k=1
F (n)=
−1
+εmn,
m−2l
n
l
n
#
"
2
2
(j)
ρk cos ζk
λk
ρk
j=1 (j)
(j)
ck
+
2rk
cos ωk n+χk
1+
1+ 2 +2
n
n
n
k=1
k=1
где
(j)
(j)
dk = rk
ζk
"
#
"
1 + ρk cos
μk #
(j)
(j)
n
cos χk + i sin χk , ωk = arg 1 +
= arccos 6
,
n
ρ2
ζk
1 + nk2 + 2 ρk cos
n
с учетом положительности sin ζk , k = 1, 2, . . . , l.
Для доказательства теоремы понадобятся следующие формулы:
n
1
λ2
λk
(j)
(j)
(j)
ck
= ck eλk − ck eλk k + o
1+
,
n
2n
n
"
#
"
μk #n ¯(j) "
μ̄k #n
(j)
(j)
(j)
(6)
+ dk 1 +
= 2rk eρk cos ζk cos ρk sin ζk + χk −
dk 1 +
n
n
#
"
1
1 (j)
(j)
− rk ρ2k eρk cos ζk cos ρk sin ζk + 2ζk + χk + o
.
n
n
98
Разложив F (n) в сходящийся ряд по степеням 1/n до членов второго порядка, используя равенства (6), находим
m−2l
l
#
"
(j) 2 λk (j) 2 ρk cos ζk
(j) ck λk e +
2rk ρk e
cos ρk sin ζk + 2ζk + χk m
1 k=1
k=1
F (n) =
+ εmn.
m−2l
l
#
"
2n j=1 (j)
(j)
(j)
λk
ρk cos ζk
ck e +
2rk e
cos ρk sin ζk + χk
k=1
k=1
Продифференцировав F (n) по n и приравняв это выражение к нулю, получим приближенное уравнение для нахождения n
3
#
"
4
(j) 2 ρk cos ζk
(j) m−2l (j) 2 λk
l
m 4
c
λ
e
+
2r
ρ
e
cos
ρ
sin
ζ
+
2ζ
+
χ
k
k
k
k=1
k=1
k
k
k
k
k
4 1
.
#
"
n≈5
(j) ρk cos ζk
(j)
m−2l (j) λk
l
2mε j=1 c
e
+
2r
e
cos
ρ
sin
ζ
+
χ
k
k
k
k
k
k=1
k=1
Осталось заметить, что
ẍj (t) xj (t) =
t=t0 +1
m−2l (j) 2 λk
k=1 ck λk e
+
m−2l (j) λk
k=1 ck e
"
#
(j)
(j)
2rk ρ2k eρk cos ζk cos ρk sin ζk + 2ζk + χk
#
"
.
(j) ρk cos ζk
(j)
l
cos ρk sin ζk + χk
k=1 2rk e
l
k=1
+
Тем самым справедливость формулы (3) установлена.
Теперь предположим, что у матрицы A имеются кратные собственные числа. В этом
случае понадобятся формулы (для кратного вещественного собственного числа λk кратности pk и для кратных комплексных собственных чисел μk и μ̄k кратностей pk )
(n−1)(n−2) . . . (n−pk +1)
pk !nk−1
n−pk
eλk
λk
λ2k +2pk λk +pk (pk −1)
1
=
1+
1−
+o
,
n
pk !
2n
n
" μ #n−pk
"
#n−pk
k
(j) (n−1)(n−2)...(n−pk +1)
¯(j) (n−1)(n−2)...(n−pk +1) 1+ μ̄k
+
d
1+
k
pk !nk−1
n
pk !npk −1
n
(j)
(j)
0
2r
r
(j)
(j)
= k eρk cos ζk cos(ρk sin ζk + χk ) − k eρk cos ζk ρ2k cos(ρk sin ζk + χk + 2ζk ) +
pk !
npk !
1
1
(j)
(j)
+ 2pk ρk cos(ρk sin ζk + χk + ζk ) + pk (pk − 1) cos(ρk sin ζk + χk ) + o
dk
=
n
в обозначениях, введенных ранее. Аналогично рассмотренному ранее случаю получаем,
что формула (3) остается справедливой и в случае наличия кратных собственных чисел
у матрицы A.
Теперь предположим, что элементы матрицы A заданы с относительными погрешностями, не превосходящими δA, а элементы вектора X0 – с относительными погрешностями, не превосходящими δX0 .
Лемма 2. В этом случае относительная погрешность вычислений для любой компоненты вектора X̃(t1 ) равна (ε + δA)n + δX0 .
Д о к а з а т е л ь с т в о. Действительно, относительная вычислительная погрешность для каждой компоненты складывается из суммы n погрешностей умножений, одной – сложения и одной – деления.
99
Следствие 2. Оптимальное число шагов n для нахождения значения решения
системы дифференциальных уравнений (1) при t1 = t0 + 1 приближенно может быть
найдено по формуле
3
⎞
⎛
4
4
x1 (t1 )
signy1 (t1 ) signy2 (t1 )
signym (t1 )
1
4
,
,...,
n≈5
A2 ⎝ . . . ⎠,
2m(ε + δA) | x1 (t1 ) | | x2 (t1 ) |
| xm (t1 ) |
xm (t1 )
где δA – максимальная относительная погрешность элементов матрицы A;
δX0 – максимальная относительная погрешность элементов вектора начальных данных X0 .
Д о к а з а т е л ь с т в о следствия 2 аналогично доказательству теоремы 1.
4. Применение полученного результата. Поскольку в правую часть формулы (4) входит значение приближенного решения в точке t1 , то возникает вопрос о возможности применения полученного результата на практике.
Для нахождения оптимального значения количества шагов интегрирования в методе Эйлера предлагается воспользоваться методом Банаха (простых итераций). Сначала определяется приближенное значение X1 (t1 ), которое получается интегрированием
с помощью метода Эйлера с одним шагом (n1 = 1). Затем по нему находится следующее
значение n2 числа шагов, и строится новое решение с числом шагов n2 и т. д. Каждое
следующее значение nk+1 числа шагов рассчитывается по формуле
⎛ (k)
⎞
3
x1 (t1 )
4
(k)
(k)
⎜ (k)
⎟
4 1
signy1 (t1 )
signym (t1 )
x2 (t1 ) ⎟
2 X (t ), X (t ) = ⎜
nk+1 = 5
,
.
.
.
,
A
⎜
⎟ . (7)
k 1
k 1
(k)
2mε | x(k) (t1 ) |
...
⎝
⎠
| xm (t1 ) |
1
(k)
xm (t1 )
Приведем простое достаточное условие сходимости последовательности (7).
Теорема 2. Для того чтобы последовательность (7) сходилась, достаточно,
чтобы
(1)
(1)
signym (t1 )
signy1 (t1 )
,...,
A2 X1 (t1 ) < 2mε.
(1)
(1)
| x1 (t1 ) |
| xm (t1 ) |
Д о к а з а т е л ь с т в о. При выполнении условия теоремы на первом шаге
получаем 0 < n2 < 1, следовательно, оптимальное число шагов метода Эйлера будет
равно 1.
В действительности справедливой является следующая теорема.
Теорема 3. Последовательность (7) всегда является сходящейся.
Д о к а з а т е л ь с т в о. Утверждение теоремы сразу же следует из непрерывности
решения системы линейных дифференциальных уравнений Ẋ = AX по начальным
данным, а также из непрерывной зависимости n от значения решения в точке t1 .
Оценим полную погрешность вычисления методом Эйлера для найденного значения n. Для этого рассмотрим, как соотносятся погрешность метода и вычислительная
погрешность для различного числа шагов метода Эйлера. При малых n погрешность
метода много больше вычислительной погрешности. Затем с увеличением числа шагов погрешности становятся сопоставимы, и тут и достигается наименьшее значение
полной погрешности. После этого при очень больших n вычислительная погрешность
становится много больше погрешности метода, которая стремится к нулю с ростом n.
100
Таким образом, при оптимальном n имеем
F (n) 2mnε.
Очевидно, что данная оценка является довольно грубой.
Пример. Определим, при каких значениях A, X0 и m – порядка матрицы A полная
погрешность будет не больше 1%, т. е. выполняется неравенство
2mnε 0.01.
Для точности float (ε = 1.19 · 10−7 ) и из условия теоремы 2
(1)
(1)
signym (t1 )
signy1 (t1 )
,...,
A2 X1 (t1 ) < 2mε
(1)
(1)
| x1 (t1 ) |
| xm (t1 ) |
получаем следующую оценку:
m < 106.
Таким образом, для порядка m < 106 полная погрешность метода составит в этом
случае не более 1%.
Как видим, даже при очень грубой оценке полной погрешности сверху существуют
матрицы, для которых вычисления с помощью метода Эйлера дают довольно точный
результат.
5. Численные примеры. Теперь рассмотрим несколько примеров.
Были проведены вычисления для матриц с различными собственными числами.
Результаты этих вычислений представлены в таблице.
Оптимальные значения числа шагов метода Эйлера
и соответствующие полные погрешности
№
Матрица,
собственные
числа
1
2
3
4
0 2
−1 3
1; 2
0 2
−1 3
1; 2
−4 3
−6 5
−1; 2
1
1
−2
−1
±i
X0
3
2
1
−3
2
1
1
2
Значение
решения
nопт
3527
4802
4293
2050
12.8199
10.1014
−29.9618
−40.8375
−6.28196
−13.6672
−1.9846
0.239192
Точное
значение
решения
12.82561976
10.10733793
−29.97713807
−40.85026538
−6.285417772
−13.67447388
−1.984110649
0.2391336272
δполн
0.00102941
0.000824155
0.00102941
0.000490666
101
Полная относительная погрешность решения
Приведем для второго примера график зависимости полной погрешности от числа
шагов метода Эйлера (рисунок).
1
0 2
.
Здесь матрица системы A =
, X0 =
−3
−1 3
Замечание 3. В этом случае при n = 32 000 полная относительная погрешность
вычислений равна 0.00402651, что превосходит найденное значение почти в 5 раз.
Заключение. В настоящей работе находится число шагов метода Эйлера, обеспечивающее наименьшую полную погрешность (сумму погрешности метода и вычислительной погрешности) значения решения задачи Коши в точке для системы линейных
дифференциальных уравнений с постоянными коэффициентами.
Как видно из приведенных примеров, решение задачи Коши для автономной системы линейных дифференциальных уравнений в точке находится довольно точно. Таким
образом, указанным в работе методом можно пользоваться в случаях, когда требуется найти решение с максимально возможной точностью. Алгоритм также может быть
применен в задачах, где требуется построить точное решение задачи Коши на каком-то
промежутке (возможно, с указанными характерными точками, где необходимо знать
точные значения решения), например при построении графиков с фиксированным временным интервалом, как они обычно и строятся.
Литература
1. Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. М.: Наука, 1987. 600 c.
2. Демидович Б. П., Марон И. А. Основы вычислительной математики. СПб.: Изд-во «Лань», 2006.
672 с.
3. Stuart A. M. Probabilistic and Deterministic Convergence Proofs for Software for Initial Value
Problems // Numerical Algorithms. 1997. Vol. 14, N 3. P. 227–260.
4. Björk A., Dahlquist G. Numerical mathematics and scientific computations: in 2 vol. Philadelphia:
SIAM, 2008. Vol. 1. XXVIII+717 p.
5. Tucker W. A Rigorous ODE Solver and Smale’s 14th Problem // Found. Comput. Math. 2002. N 2.
P. 53–117.
102
6. Higham N. J. Accuracy and stability of numerical algorithms. Philadelphia: SIAM, 1996. 718 p.
7. Golub G. H., Ortega J. M. Scientific Computing and Differential Equations. San Diego: Academic
Press, 1992. 335 p.
8. Ортега Дж., Пул У. Введение в численные методы решения дифференциальных уравнений / пер. с англ. Н. Б. Конюховой, под ред. А. А. Абрамова. М.: Наука, 1986. 288 с.
Статья рекомендована к печати проф. Л. А. Петросяном.
Статья принята к печати 28 мая 2009 г.
1/--страниц
Пожаловаться на содержимое документа