§ 1.5. Анализ погрешностей

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

1.5.1. Пример.

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

xў = lx

(1)

с l = -10³. Упомянутая оценка содержит быстро растущий с ростом T множитель » (e^|l|T - 1)/|l| (поскольку константа Липшица для правой части уравнения (1) равна, очевидно, |l| и, следовательно, L = |l|). Поэтому для вычисления решения на больших промежутках в соответствии с заключением теоремы 1.4.6 требуется сильное уменьшение шага: t = O(e[M(e^|l|T - 1)/|l]^-1) » O(ee^-lT), где e — требуемая точность. Например, при e » 1 и T » 1 шаг должен быть порядка e^-1000 » e^-400.

С другой стороны, в силу того, что уравнение (1) экспоненциально устойчиво, явный метод Эйлера обладает свойством уменьшения погрешности на каждом шаге. Для того чтобы понять это явление, рассмотрим простейшую ситуацию. Допустим мы ищем нулевое решение уравнения (1) и обозначим через e₀ погрешность при выборе начального условия: e₀ = x₀ - 0 = x₀. Проследим, как меняется эта погрешность в процессе счета (предполагая, для простоты, что мы ведем вычисления без ошибок округления). Имеем

ei = xi - 0 = xi-1 + ltxi-1 = (1 + lt)xi-1 = (1 + lt)(xi-1 - 0) = (1 + lt)ei-1,

(2)

откуда получаем

ei = (1 + lt)ie0.

Так как l < 0, при 0 < t Ј -2/l множитель (1 + lt)ⁱ по модулю не превосходит единицы и, следовательно, погрешность не превосходит e₀ независимо от длины T промежутка, на котором ведутся вычисления (более того, погрешность e стремится к нулю при i ® Ґ).

1.5.2. В чем причина?

Суть этого явления, огрубляя ситуацию, можно описать так. При выводе оценки (2) мы считали, что погрешность с шага на шаг переносится, подчиняясь уравнению (1), а при выводе оценки (1.4.16) — подчиняясь уравнению

eў = Le,

где L — константа Липшица функции F (рис. 5).

Рис. 5.

Таким образом, метод, использованный при доказательстве теоремы 1.4.6, не учитывает более тонких свойств уравнения (E), а именно, свойств устойчивости. Поэтому мы сейчас сначала оценим погрешность на каждом шаге, а затем, учитывая свойства уравнения (E), оценим общую погрешность метода.

1.5.3. Локальная и глобальная погрешности.

Пусть (t, x) О R×R^m — произвольная точка, j — решение уравнения (E), проходящее через эту точку, а j_t — приближенное решение, начинающееся в точке (t, x), определяемое явным одношаговым методом

xi = xi-1 + tF(ti-1, xi-1, t).

(3)

Локальной погрешностью метода (3) (в точке (t, x)) называется величина

e(t) = e(t, x, t) = j(t + t) - (jt)1

(рис. 6). Очевидно,

Рис. 6.

e(t, x, t) = x + tr(t, x, t) - x - tF(t, x, t) = t[r(t, x, t) - F(t, x, t)]

(напомним, что r(t, x, t) = [j(t + t) - j(t)]/t, если t > 0). Условие (1.4.15) в теореме 1.4.6 в точности означает, что

e(t, x, t) Ј Mtk+1.

(4)

Как следует из примера п. 1.4.7, локальная погрешность характеризует свойства аппроксимации метода.

Глобальной погрешностью метода (3) называют величину

En(t) = j(tn) - (jt)n = j(tn) - jt(tn).

Несколько огрубляя ситуацию, можно говорить, что в теореме 1.4.6 утверждается, что

||En(t)|| Ј  eLT - 1) L M||e(t)||.

В следующем пункте мы выведем более тонкую оценку глобальной погрешности через локальные.

1.5.4. Перенос погрешностей.

Пусть j — как обычно, решение задачи (E) – (C), а j_t — решение явной одношаговой схемы (3). Пусть, кроме того, для простоты, уравнение (E) скалярное. Обозначим через jⁱ (i = 1, ..., n) решение

Рис. 7.

уравнения (E), удовлетворяющее начальному условию x(t_i) = (j_t)_i (рис. 7). На рисунке видно, что глобальная погрешность E_n(t) получается в результате суммирования "перенесенных" вдоль решений уравнения (E) локальных погрешностей на каждом шаге:

En(t) = j(tn) - (jt)n =  n е i = 0 [ji-1(tn) - ji(tn)].

(5)

Оценим каждое слагаемое. Поскольку jⁱ — решения уравнения (E),

(ji-1)ў(t) - (ji)ў(t) є f[t, ji-1(t)] - f[t, ji(t)].

(6)

Обозначим j^i-1 - jⁱ через yⁱ, а (f[t, j^i-1(t)] - f[t, jⁱ(t)])/[j^i-1(t) - jⁱ(t)] через aⁱ(t).

Задача 1.5.1. Почему если jⁱ(t_i) - j^i-1(t_i) № 0, то jⁱ(t) - j^i-1(t) № 0 при всех t?

В этих обозначениях равенство (6) перепишется в виде

(yi)ў(t) = ai(t)yi(t),

откуда

yi(tn) = yi(ti)·exp ж и т tn ti ai(s) ds ц ш ,

или

ji-1(tn) - ji(tn) = e(ti-1, xi-1, t)·exp ж и т tn ti ai(s) ds ц ш ,

Поэтому, продолжая (5), получим

En(t) =  n е i = 1 e(ti-1, xi-1, t)·exp ж и т tn ti ai(s) ds ц ш .

Именно второй множитель в каждом слагаемом содержит информацию о свойствах дифференциального уравнения, позволяющую оценивать как переносятся погрешности вдоль его решений. В самом деле, по теореме Лагранжа

ai(s) = fx[s, xi(s)],

где xⁱ(s) О (j^i-1(s), jⁱ(s)) и, таким образом,

En(t) =  n е i = 1 e(ti-1, xi-1, t)·exp ж и т tn ti fx[s, xi(s)]ds ц ш .

(7)

Например, если о функции f неизвестно ничего, кроме непрерывности и условия Липшица, то (см. п. 1.1.1) |f_x(t, x)| Ј L при всех (t, x). Поэтому, в случае выполнения оценки (4) (т. е. условия (1.4.15) в теореме 1.4.6), оценка (7) может быть продолжена так:

En(t) Ј  n е i = 1 Mtk+1·exp  ж и т tn ti–1 L·ds ц ш .

Задача 1.5.2. Докажите, что из последней оценки вытекает оценка E_n(t) Ј e^LT·Ct^k с некоторой константой C.

Задача 1.5.3. Как можно получить аналогичные оценки в случае системы, т. е. при m і 2?

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

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

Пусть правая часть (скалярного) уравнения (E) дифференцируема по x и для некоторого L > 0 и всех (t, x) выполнено неравенство

fx(t, x) Ј -L.

(8)

Пусть, кроме того, выполнена оценка (4) для локальной погрешности. Тогда

En(t) Ј  e-Lt 1 - e-Lt Mtk+1.

(9)

Д о к а з а т е л ь с т в о.  В силу (8)

т tn ti fx[s, xi(s)]ds Ј  т tn ti (-L )ds = -L (tn - ti-1) = -L (n - i + 1)t.

Поэтому, продолжая (7), с учетом (4), получаем

En(t) Ј  n е i = 1 Mtk+1e-L(n-i+1)t = Mtk+1 n е i = 1 e-Lit <

< Mtk+1 Ґ е i = 1 (e-Lt)i =  e-L t 1 - e-Lt Mtk+1,

что и требовалось.

1.5.6. Замечания.

а) Оценка (9) в отличие от оценки (1.4.16) не ухудшается с ростом промежутка, на котором ищется решение.

б) (k+1)-й порядок малости по t в правой части оценки (9) лишь кажущийся; на самом деле, поскольку lim_t®0[t/(1 - e^-Lt)] = -1/L, правая часть этой оценки имеет порядок O(t^k).

в) Условие (8) для скалярного уравнения в точности означает, что все решения уравнения (E) являются экспоненциально устойчивыми.

г) Аналогом условия (8) в многомерном случае (m > 1) может служить условие

Re l Ј -L

при всех l О s[f_x(t, x)], где s(A) — спектр оператора A. Другими словами, все собственные значения оператора f_x(t, x) при всех (t, x) лежат в полуплоскости {l О C: Re l Ј L}.

1.5.7. Главный член локальной погрешности.

Условие (4) означает (при достаточной гладкости входящих в наши рассмотрения функций), что

¶le(t, x, t) ¶tl к к t=0  = 0 при 0 Ј l Ј k.

Поэтому, раскладывая e по степеням t по формуле Тейлора до (k+1)-го порядка в точке t = 0, имеем

e(t, x, t) = h(t, x)tk+1 + O(tk+2).

(10)

Величина h(t, x)t^k+1 называется главным членом локальной погрешности. Она может быть вычислена через решение и его производные, а следовательно, через производные правой части.

Задача 1.5.4. Покажите, что для явного метода Эйлера

h(t, x) =  1 2 (ft + fx f).

(11)

Найдите h(t, x) для метода предиктор-корректор

1.5.8. Задача контроля локальной погрешности.

При проведении вычислений весьма желательна информация о величине погрешности, вносимой на каждом шаге. С одной стороны, мы должны следить, чтобы погрешность не была слишком большой, чтобы в результате не получить решение с неудовлетворяющей нас точностью. С другой стороны, погрешность, следуя принципу "лучшее — враг хорошего", не должна быть слишком маленькой; иначе будет проведен большой объем излишней вычислительной работы.

Теорема 1.4.6, а для устойчивых уравнений — теорема 1.5.5, позволяет оценивать глобальную погрешность через локальную. Поэтому один из способов контроля точности вычислений заключается в контроле локальной погрешности. Если локальная погрешность излишне мала, шаг можно увеличить, если велика — уменьшить. Естественно считать, что за величину локальной погрешности отвечает ее главный член. Однако использовать для его вычисления формулы типа формулы (11) крайне невыгодно, поскольку они требуют знания высших производных функции f (методы же Рунге — Кутты придумывались именно для того, чтобы избавиться от вычисления этих производных). Имеется ряд приемов, позволяющих обойтись без формул для вычисления главного члена локальной погрешности. Один такой прием мы описываем ниже.

1.5.9. Оценка локальной погрешности с помощью экстраполяции Ричардсона.

Нашей задачей является оценка выражения h(t, x)t^k+1 в произвольной точке (t, x). С помощью метода (3) сделаем сначала один шаг длины t из точки (t, x), обозначив полученное значение решения через (j_t)₁, а затем из той же точки (t, x) сделаем два шага длины t/2, обозначив результат через (j_t/2)₂. Обе полученные величины являются приближениями решения j в точке t + t. Погрешность вычисления (j_t)₁ есть локальная погрешность в точке (t, x) и, в соответствии с (10),

j(t + t) = (jt)1 + h(t, x)tk+1 + O(tk+2).

Обозначим для краткости h(t, x) через H:

j(t + t) = (jt)1 + Htk+1 + O(tk+2).

(12)

Погрешность же вычисления (j_t/2)₂ складывается (см. (7)) из "перенесенной" локальной погрешности e(t, x, t/2) и локальной погрешности метода в точке (t + t/2, (j_t/2)₁):

j(t + t) = (jt/2)2 + E2 ж и t 2 ц ш  = (jt/2)2 +

+ e ж и t, x,  t 2 ц ш ·exp ж и т t+t t+t/2 fx[s, x1(s)] ds ц ш  + e й л t +  t 2 , (jt/2)1,  t 2 щ ы .

Заметим теперь, что в этой формуле (мы используем обозначение f_x = f_x(t, x))

e ж и t, x,  t 2 ц ш  = H· ж и t 2 ц ш k+1  + O(tk+2),

(13)

exp ж и т t+t t+t/2 fx[s, x1(s)] ds ц ш  = exp ж и т t+t t+t/2 [fx + O(t)] ds ц ш  = = exp ж и t 2 fx + O(t2) ц ш  = 1 +  t 2 fx + O(t2),

(14)

e ж и t +  t 2 , (jt/2)1,  t 2 ц ш  =h ж и t +  t 2 , (jt/2)1 ц ш ж и t 2 ц ш k+1 + + O(tk+2) = [H + O(t)] ж и t 2 ц ш k+1  + O(tk+2).

(15)

Задача 1.5.5. Докажите формулы (13) – (15).

Тогда, учитывая (13) – (15),

j(t + t) = (jt/2)2 + E2 ж и t 2 ц ш  = (jt/2)2 +  й л H· ж и t 2 ц ш k+1  + O(tk+2) щ ы  × ×  й л 1 +  t 2 fx + O(t2) щ ы  + [H + O(t)] ж и t 2 ц ш k+1  + O(tk+2) =  = (jt/2)2 + 2H ж и t 2 ц ш k+1  + O(tk+2).

(16)

Отбрасывая в (12) и (16) члены (k+2)-го порядка малости, можно, во-первых, "найти" значение H и, во-вторых, "найти" значение j(t + t). Кавычки здесь означают, что мы можем найти эти значения с точностью O(t^k+2). Для этого нужно решить систему

j~= (jt)1 + h~tk+1,

j~= (jt/2)2 + h~tk+1/2k

относительно неизвестных j~и h~:

h~tk+1=  2k 2k - 1 [(jt/2)2 - (jt)1],

j~= (jt/2)2 +  (jt/2)2 - (jt)1 2k - 1 .

Теперь величина h~tk+1может использоваться (с известной осторожностью) взамен главного члена погрешности h(t, x)tk+1 и применяться для контроля погрешности на шаге, а величина j~приближает значение j(t + t) с увеличенным на единицу порядком точности.

File based on translation from T_EX by T_TH, version 3.05.
28 May 2002, 11: 10.