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

Сер. 10. 2009. Вып. 4

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

УДК 519.622.2

Е. А. Калинина, О. Н. Самарина

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

1. Введение. На практике часто приходится иметь дело с системами линейных дифференциальных уравнений с постоянными коэффициентами или с системами, в которых коэффициенты можно считать постоянными в течение довольно больших промежутков изменения переменной. В действительности при построении различных систем прогнозирования, АСУ, тренажеров и подобных систем разрабатываемые программные комплексы просто решают задачу интегрирования этой системы с постоянным шагом интегрирования. Иначе говоря, для решения задачи Коши применяется метод Эйлера как наиболее простой и наименее трудоемкий. Как известно, при вычислении значения решения в точке с помощью метода Эйлера погрешность получается довольно большой. Однако существуют только оценки погрешности метода (см., например, [1, 2]), нет достаточно хороших оценок полной погрешности - суммы погрешности метода и вычислительной погрешности. Но при уменьшении шага интегрирования, что ведет к снижению погрешности метода, вычислительная погрешность становится весьма существенной и не принимать ее во внимание нельзя. Можно было пренебрегать вычислительной погрешностью, когда уровень развития вычислительной техники не позволял выполнять столь большое число операций, что влияние вычислительной погрешности на результат становилось значительным. В современных условиях, когда компьютер выполняет вычисления в режиме реального времени, нужно учитывать и ошибки округления. В связи с этим появился вероятностный подход к решению задачи оценки полной погрешности [3]. К сожалению, он не всегда применим, поскольку ошибки округления нельзя считать независимыми случайными величинами [4]. Также применялся метод поиска точного решения системы [5]. Другие подходы не использовались в связи с большой сложностью задачи. На это указывают многие авторы [4, 6-8]. В настоящей работе с помощью оценки верхней границы относительной погрешности явно находится оптимальное значение шага интегрирования, обеспечивающее наименьшую полную погрешность. Показывается, что для некоторого круга систем линейных дифференциальных уравнений интегрирование с помощью метода Эйлера допустимо, оно обеспечивает нахождение значения решения в точке с достаточно малой погрешностью. Приводятся численные примеры.

Калинина Елизавета Александровна — кандидат физико-математических наук, доцент кафедры высшей математики факультета прикладной математики—процессов управления Санкт-Петербургского государственного университета. Количество опубликованных работ: 11. Научные направления: теория управления, устойчивость, компьютерная алгебра, символьные (аналитические) методы. E-mail: [email protected].

Самарина Ольга Николаевна — аспирант факультета прикладной математики-процессов управления. Научный руководитель: доктор физико-математических наук, проф. А. Ю. Утешов. Количество опубликованных работ: 2. Научные направления: компьютерная алгебра, символьные (алгебраические) методы. E-mail: [email protected].

© Е. А. Калинина, О. Н. Самарина, 2009

2. Основной результат. Для решения задачи будем использовать следующую норму вектора:

І Х1

X,

m

xl

: \\Х || = \х! \ + 1x21 +-+ \хт |

\ хп

и соответствующую матричную норму

РН = 1т1.ах ^3 ^

где А = (А1, А2,..., Ат) - матрица порядка т.

Рассмотрим задачу Коши для системы линейных дифференциальных уравнений

Х = АХ, X (10)= Х0 (Хо = 0), (1)

здесь А = [а^]т3=1 - вещественная квадратная матрица порядка т с постоянными элементами,

/ х1(Ь) \ / х1

X (t) =

xo

m

Хо =

\ хт(Ь)

Будем считать, что все элементы матрицы А и вектора Хо заданы точно.

Известно, что решение задачи (1) имеет вид Х = ел(г-Ьо)Хо. Найдем полную относительную погрешность приближенного решения задачи Коши методом Эйлера в точке tl = tо + 1.

При данном числе шагов п получим следующее приближенное значение решения:

Х(пХо. (2)

Пусть AJ - жорданова нормальная форма матрицы А, а Т = [Ь^]"3=1 - матрица, составленная из векторов канонического базиса, такая, что AJ = Т-1АТ.

Тогда равенство (2) можно записать в виде

X(t1)^T[E+-Aj) Т-1Хо.

1

n

Следует учесть, что вектор X(ti) может быть нулевым только при достаточно больших значениях ti, и в этом случае можно говорить о стремлении решения задачи Коши (1) к нулю с ростом t. Если же в нуль обращаются не все компоненты вектора XT (ti), то можно попытаться избежать этой ситуации изменением числа шагов n метода Эйлера, например увеличить или уменьшить n на 1. Обращение в нуль части компонент вектора X(ti) возможно лишь в исключительных случаях, поскольку наименьшие модули чисел, отличных от нуля, имеют следующие значения: для вычислений с одинарной точностью float « 1.17 • 10~38, для двойной точности double « 2.225 • 10~308, для long double « 3.36 • 10~4932.

Далее рассмотрим общий случай: ни одна из компонент вектора X(ti) не обращается в нуль.

Оптимальное (в смысле наименьшей полной погрешности) значение числа шагов метода Эйлера можно найти, воспользовавшись следующей теоремой.

Теорема 1. Оптимальное число шагов п для нахождения значения решения системы дифференциальных уравнений (1) при ¿і = і0 + 1 приближенно может быть найдено по формуле

Xl(íl) Xl(íl) + x2(ti) x2(íl) + ••• + (^1) Xm (^1)

2тє

(3)

Замечание 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) вектор вторых производных:

Y (t) =

Xi(t)

Xm (t)

= A2X (t).

Тогда справедливо следующее утверждение.

Следствие 1. Оптимальное число шагов п для нахождения значения решения системы дифференциальных уравнений (1) при ^ + 1 приближенно может быть найдено по формуле

1 (signyi(ti) signy2(ti) signym(ti)

A2

( xi(ti) \ X2 (ti)

У xm (ti) J

(4)

2тє \ I xi(ti) I ’ I X2(ti) I I xm(ti) \

Здесь signy - знак числа y:

{1, если y > 0,

0, если y = 0,

— 1, если y < 0.

Замечание 2. При нахождении значения решения задачи Коши (1) в точке ti = to + т методом Эйлера имеем

вд = (е+^у*0.

Таким образом, этот случай сводится к рассмотренному случаю т =1 заменой матрицы A на матрицу Ат. Следовательно, достаточно рассмотреть только случай т =1.

3. Доказательство теоремы 1. Рассмотрим сначала случай различных собственных чисел матрицы A.

Пусть А - квадратная матрица порядка m - имеет l пар комплексно-сопряженных собственных чисел ці, li,..., ці,h и m — 2l вещественных собственных числа Xi,X2,...,Xm-2¡. Для комплексных собственных чисел будем также использовать тригонометрическую форму ц = pj(cos Q + i sin Q), lj = pj (cos(— Q) + i sin(—Q)), j = 1, 2,...,l, причем будем считать sinZj > 0.

*) Значения £ взяты из стандартного включаемого файла АоаЬ.Н для С-компилятора на архитектуре x86.

n

Тогда ^’-тая компонента вектора X (¿1) - решения задачи Коши (1) будет иметь вид с[з)вХ1 + е()еХ2 +-+ ¿¿—21 еХт-2‘ + еМ1 + еД1 +------+ ет ^ещ .

Соответственно ^’-тую компоненту вектора X (¿1) - приближенного решения задачи Коши методом Эйлера можно записать так:

Si)

1 + —) +--- + с(£_21

i +

n

+

+ d

(i)

fi\n

i

+ d

(i)

i + ^У + Ф (i + Щ

n J \ n J

(5)

Обозначим относительную погрешность вычисления ^’-й компоненты вектора ХГ(¿1) по формуле (5) через 51(и).

Лемма 1. 5]_(и) = пе.

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

Таким образом, нужно найти значение п, при котором достигается минимум функции

F (n)=]T

j=1

cj)eXl + .. . +сШ—:lexm-21 + djeP1 +d11j)ep1 + . .. +d(j)ePl +d(j)ePi

или

F (n)=E

j = 1

ш — 21 i ckj) f 1

k =1 V

ш — 21

E« k=1

-1

k=1

f^k

n

+ £mn,

^ ckj)eXfc + ^2rj)e№ cosCfc cos (pfc sin Cfc + x|j))

k = 1

ш—2l

k=1

n l

--1

1+ii +S2r<» i+a+2

k=1

+emn,

где

j(j) _ Jj) Ai) 1 „• ^ ,,0'Л .. „„„ Л , Mfc\ 1 +

<4 = rk [cosXk + *sinXfc ) , ^fc = arg ( Ц-------------------) = arccos

£fc_cosCfc

^ _L £js. \ 2 Pfc CQS Cfc

П2 П

с учетом положительности sin (k, & = 1, 2,...,l.

Для доказательства теоремы понадобятся следующие формулы:

4j)(i + -Y + 4J} (1 + —) = cos Cfc eos (Рк sin а+xf) -

&)

n

fk \

(i)

n

(i)

(6)

^rí3)PkePkCOS<:k cos (pfcSÍnCfc + 2Cfc +X(fcj)) +0 .

n

Лш — 2l

1

1

n

n

n

n

Разложив F (n) в сходящийся ряд по степеням 1/n до членов второго порядка, используя равенства (6), находим

j = 1

г-21

Е ckj)XleXk ^E2rkj)PÍePk COS Cfc cos (pk sin Zk + 2Zk + xïj))

k=1

k = 1

г-21

^2 cj)eXk + ^2rj)ePk COs Zk cos p sin Zk + xkP)

k=1

k=1

+ emn.

Продифференцировав F(п) по п и приравняв это выражение к нулю, получим приближенное уравнение для нахождения п

1

2me

Е

ЕГ=12г ckj)xîeXk + EU 2rk )p\ePk COs Zk cos (pksin Zk + 2<Ck + xkj))

J=1

Осталось заметить, что

ЕГ=12г c(kj)eAfc + Ek=i 2rj)epk COs Zk cos (pk sin Zk + X(kj)

xj(t)

i(t)

t = tü + 1

ЕГ=12г ck)xleXk + EU 2rk]plePk cosCk cos (рл sinCfc + 2ü + x|j))

EU^ 4j)eAfc + EU Zr^ePb cos & cos ipfc sin Cfc + x|j)

Тем самым справедливость формулы (3) установлена.

Теперь предположим, что у матрицы А имеются кратные собственные числа. В этом случае понадобятся формулы (для кратного вещественного собственного числа Хь кратности рь и для кратных комплексных собственных чисел и кратностей рь)

(n-1)(n-2) ...(n-pk+1) Л , А

Pk!nk 1

n-pk

1+— ) = —г ( i-

Pk !

eXk (-, X\+2pkXk+Pk(Pk-1)

2n

+ .|i

n

d(

du

(j) (n—l)(n—2)...(n—pfc+l) / /ifc\"-№ j(j) (n—l)(n—2)...(n—pfc+l)

2r

(j)

pk ! nk-1

1+^J

(j)

pk ! npk

1

1+^

n

Ve№ œsCfc cos(/°fc sin Cfc + xU - ^-ге№ cosCfc pi cos(pfc sin Cfc + xïP + Kk) +

r

Pk! k npk !

(j)

(j)

+ 2pkpk cos(pk sin Zk + xkj) + Zk) + Pk(Pk - 1) cos(pk sin Zk + Xj))

+ o | -n

в обозначениях, введенных ранее. Аналогично рассмотренному ранее случаю получаем, что формула (3) остается справедливой и в случае наличия кратных собственных чисел у матрицы А.

Теперь предположим, что элементы матрицы А заданы с относительными погрешностями, не превосходящими 5А, а элементы вектора Хо - с относительными погрешностями, не превосходящими ЗХо.

Лемма 2. В этом случае относительная погрешность вычислений для любой компоненты вектора Х^) равна (е + 5А)п + 5Х0.

Доказательство. Действительно, относительная вычислительная погрешность для каждой компоненты складывается из суммы п погрешностей умножений, одной - сложения и одной - деления.

n

x

n-pk

1

Следствие 2. Оптимальное число шагов п для нахождения значения решения системы дифференциальных уравнений (1) при ¿1 = ¿0 + 1 приближенно может быть найдено по формуле

(signyl(íl) signy2(íl) signym(íl)

2ш(є + ¿А) \ | хі(іі) | ’ | Ж2(іі) |

| хт(іі) \

А2

Хі(іі)

хт(іі)

где 5А - максимальная относительная погрешность элементов матрицы А; 5Х0 - максимальная относительная погрешность элементов вектора начальных данных Х0.

Доказательство следствия 2 аналогично доказательству теоремы 1.

4. Применение полученного результата. Поскольку в правую часть формулы (4) входит значение приближенного решения в точке ¿1, то возникает вопрос о возможности применения полученного результата на практике.

Для нахождения оптимального значения количества шагов интегрирования в методе Эйлера предлагается воспользоваться методом Банаха (простых итераций). Сначала определяется приближенное значение Х1(^), которое получается интегрированием с помощью метода Эйлера с одним шагом (п1 = 1). Затем по нему находится следующее значение п2 числа шагов, и строится новое решение с числом шагов п2 и т. д. Каждое следующее значение пь+1 числа шагов рассчитывается по формуле

пк + і —

/ (к) I signy1

(¿і)

2шє V | х[к) (¿і) |

А2Хк(іі), Хк(¿і) —

( х[к)(іі) N х2к) (іі)

(к) \ хт

(7)

(іі) /

Приведем простое достаточное условие сходимости последовательности (7). Теорема 2. Для того чтобы последовательность (7) сходилась, достаточно, чтобы

„(і)

(іі)

„(і)

(іі)

А2Хі(іі) < 2шє.

Доказательство. При выполнении условия теоремы на первом шаге получаем 0 < п2 < 1, следовательно, оптимальное число шагов метода Эйлера будет равно 1.

В действительности справедливой является следующая теорема.

Теорема 3. Последовательность (7) всегда является сходящейся.

Доказательство. Утверждение теоремы сразу же следует из непрерывности решения системы линейных дифференциальных уравнений X = АХ по начальным данным, а также из непрерывной зависимости п от значения решения в точке ¿1.

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

1

п

1

Таким образом, при оптимальном n имеем

F(n) ^ 2mne.

Очевидно, что данная оценка является довольно грубой.

Пример. Определим, при каких значениях A, Хо и m - порядка матрицы A полная погрешность будет не больше 1%, т. е. выполняется неравенство

2mne ^ 0.01.

Для точности float (е = 1.19 • 10~7) и из условия теоремы 2

Vi xi1)(ii) і і хт1 (ti) і

41>(i0 i ^(i)

получаем следующую оценку:

т < 106.

Таким образом, для порядка т < 106 полная погрешность метода составит в этом случае не более 1%.

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

5. Численные примеры. Теперь рассмотрим несколько примеров.

Были проведены вычисления для матриц с различными собственными числами. Результаты этих вычислений представлены в таблице.

Оптимальные значения числа шагов метода Эйлера и соответствующие полные погрешности

№ Матрица, собственные числа Хо и-опт Значение решения Точное значение решения ¿полн

1 (-! D 1;2 ю со 3527 / 12.8199 \ \ 10.1014 ) / 12.82561976 \ \ 10.10733793 ) 0.00102941

2 (-! D 1;2 (-S) 4802 ( -29.9618 \ \ -40.8375 ) ( -29.97713807 \ \ -40.85026538 ) 0.000824155

3 (- 5) -1;2 (?) 4293 ( -6.28196 \ \ -13.6672 ) ( -6.285417772 \ \ -13.67447388 ) 0.00102941

4 (::?) ±г (5) 2050 ( -1.9846 \ \ 0.239192 ) ( -1.984110649 \ \ 0.2391336272 ) 0.000490666

7-10 3

6-10 3

510-3

410-3

310“3

210-3 10 3

0 104

3-104 5-104 7-104

Полная относительная погрешность решения

Приведем для второго примера график зависимости полной погрешности от числа шагов метода Эйлера (рисунок).

Замечание 3. В этом случае при n = 32 000 полная относительная погрешность вычислений равна 0.00402651, что превосходит найденное значение почти в 5 раз.

Заключение. В настоящей работе находится число шагов метода Эйлера, обеспечивающее наименьшую полную погрешность (сумму погрешности метода и вычислительной погрешности) значения решения задачи Коши в точке для системы линейных дифференциальных уравнений с постоянными коэффициентами.

Как видно из приведенных примеров, решение задачи Коши для автономной системы линейных дифференциальных уравнений в точке находится довольно точно. Таким образом, указанным в работе методом можно пользоваться в случаях, когда требуется найти решение с максимально возможной точностью. Алгоритм также может быть применен в задачах, где требуется построить точное решение задачи Коши на каком-то промежутке (возможно, с указанными характерными точками, где необходимо знать точные значения решения), например при построении графиков с фиксированным временным интервалом, как они обычно и строятся.

Литература

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

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. Bjork 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.

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 г.