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

АНАЛИЗЪ м. греч. разборъ, раздробка, разрЂшенiе, разложенiе цЂлаго на составныя части его... Анализировать что, разлагать, разбирать цЂлое на начала, основы, стихiи...

Владимiр Даль.. Толковый словарь живаго великорусскаго языка

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

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

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

ε_i = x_i – 0 = x_i–1 + λτx_i–1 = (1 + λτ)x_i–1 = (1 + λτ)(x_i–1 – 0) = (1 + λτ)ε_i–1,

(2)

Так как λ < 0, при 0 < τ ≤ –2/λ множитель (1 + λτ)ⁱ по модулю не превосходит единицы и, следовательно, погрешность не превосходит ε₀ независимо от длины T промежутка, на котором ведутся вычисления (более того, погрешность ε стремится к нулю при i → ∞).

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

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

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

x_i = x_i–1 + τΦ(t_i–1, x_i–1, τ).

(3)

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

ε(t, x, τ) ≤ Mτ^k+1.

(4)

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

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

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

E_n(τ) = φ(t_n) – (φ_τ)_n =

n
∑
i = 0

[φ^i–1(t_n) – φⁱ(t_n)].

(5)

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

(φ^i–1)′(t) – (φⁱ)′(t) ≡ f[t, φ^i–1(t)] – f[t, φⁱ(t)].

(6)

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

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

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

φ^i–1(t_n) – φⁱ(t_n) = ε(t_i–1, x_i–1, τ)·exp

(

∫

t_n

t_i

aⁱ(s)ds

)

E_n(τ) =

n
∑
i = 1

ε(t_i–1, x_i–1,τ)·exp

(

∫

t_n

t_i

aⁱ(s)ds

)

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

где ξⁱ(s) ∈ (φ^i–1(s), φⁱ(s)) и, таким образом,

E_n(τ) =

n
∑
i = 1

ε(t_i–1, x_i–1,τ)·exp

(

∫

t_n

t_i

f_x[s, ξⁱ(s)]ds

)

(7)

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

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

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

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

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

f_x(t, x) ≤ –Λ.

(8)

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

E_n(τ) ≤

e^–Λτ

1 – e^–Λτ

Mτ^k+1.

(9)

∫

t_n

t_i

f_x[s, ξⁱ(s)]ds ≤

∫

t_n

t_i

(–Λ )ds = –Λ (t_n – t_i–1) = –Λ (n – i + 1)τ.

E_n(τ) ≤

n
∑
i = 1

Mτ^k+1e^{–Λ(n–i+1)τ} =Mτ^k+1

n
∑
i = 1

e^–Λiτ <

<Mτ^k+1

∞
∑
i = 1

(e^–Λτ)ⁱ =

e^{–Λ τ}

1 – e^–Λτ

Mτ^k+1,

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

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

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

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

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

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

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

ε(t, x, τ) = h(t, x)τ^k+1 + O(τ^k+2).

(10)

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

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

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

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

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

φ(t + τ) = (φ_τ)₁ + Hτ^k+1 + O(τ^k+2).

(12)

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

+ ε

(

t, x,

)

·exp

(

∫

t+τ

t+τ/2

f_x[s, ξ¹(s)] ds

)

+ ε

[

t +

, (φ_τ/2)₁,

]

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

(

t, x,

)

= H·

(

)

^k+1

+ O(τ^k+2),

(13)

exp

(

∫

t+τ

t+τ/2

f_x[s, ξ¹(s)] ds

)

= exp

(

∫

t+τ

t+τ/2

[f_x + O(τ)] ds

)

= exp

(

f_x + O(τ²)

)

= 1 +

f_x + O(τ²),

(14)

(

t +

, (φ_τ/2)₁,

)

(

t +

, (φ_τ/2)₁

)

(

)

^k+1

+ O(τ^k+2) =[H + O(τ)]

(

)

^k+1

+ O(τ^k+2).

(15)

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

φ(t + τ) = (φ_τ/2)₂ + E₂

(

)

= (φ_τ/2)₂ +

[

H·

(

)

^k+1

+ O(τ^k+2)

]

[

1 +

f_x + O(τ²)

]

[H + O(τ)]

(

)

^k+1

+ O(τ^k+2) =

= (φ_τ/2)₂ + 2H

(

)

^k+1

+ O(τ^k+2).

(16)

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

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