Глава 1. Задача Коши

Назад § 1.3. Методы Рунге — Кутты Вперед

Признаюсь, я с некоторою торжественностью возвещаю об этом новом лице.

Ф.М. Достоевский. Село Степанчиково и его обитатели

В этом параграфе описываются и исследуются явные методы Рунге — Кутты, а также кратко - неявные методы Рунге — Кутты.

1.3.1. Задача повышения порядка сходимости.

В приложениях особенно важна задача повышения порядка сходимости разностных схем. Разностные схемы более высокого порядка позволяют уменьшить шаг сетки. Например, если при заданной точности ε схема первого порядка требует шага O(ε), то схема четвертого порядка — O(4ε). Практика показывает, что разностные схемы высокого порядка оказываются особенно эффективными при проведении прецизионных (с большой требуемой точностью) расчетов. Более того, схемами низких порядков такие расчеты часто вообще провести невозможно даже при сколь угодно малом (допустимом в ЭВМ) шаге.

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

В этом и последующих параграфах в случаях, когда это необходимо, мы считаем уравнение (E) скалярным, т. е. считаем, что m = 1.

1.3.2. Схема Хойна, или предиктор-корректор.

Попытаемся повысить порядок аппроксимации явной схемы Эйлера следующим образом. В исходном уравнении (E) заменим производную обычным конечно-разностным соотношением, а f(t, x) постараемся аппроксимировать "с более высоким порядком". Например, возьмем вместо f(ti–1, xi–1) "среднее направление" между векторами поля направлений уравнения (E) в точках (ti–1, xi–1) (вектор f1 на рис. 3) и (ti, x*i), где (ti, x*i) следующая вершина ломаной Эйлера (вектор f2). Интуитивно ясно, что получившееся направление (вектор f3) ближе к "истинно нужному" направлению (вектор f4).

К построению метода Хойна
Рис. 3.

Аналитически эта идея реализуется в виде разностной схемы (S) с разностным оператором

[Fτ(x)]i = {
xixi–1
τ
  –  1
2
(
f(ti–1, xi–1) + f[ti, xi–1 + τf(ti–1, xi–1)]

),  если i = 1, ..., n,
xix0,  если i = 0;

Эта схема (которую иногда называют разностной схемой Хойна) также является явной в том смысле, что ее решение вычисляется по рекуррентным формулам:

τ)0 = x0,

τ)i = (φτ)i–1 +  τ
2
( f[ti–1, (φτ)i–1] +  f(ti, (φτ)i–1 + f[ti–1, (φτ)i–1])) .

Разностную схему Хойна (вернее рекуррентный переход от xi–1 к xi) часто записывают в виде двух полушагов:
x*i= xi–1 + τf(ti–1, xi–1),(1а)
xi = xi–1 +  τ
2
[f(ti–1, xi–1) + f(ti, x*i)].
(1б)
Поэтому она обычно называется схемой предиктор-корректор: на первом полушаге (1а) приближенное решение "предсказывается" (от англ. to predict предсказывать) с первым порядком точности, а на втором (1б) "корректируется" (от to correct исправлять, корректировать) с целью повышения точности.

1.3.3. Аппроксимация схемы предиктор-корректор.

Покажем, что если f в (E) дважды непрерывно дифференцируема, то схема предиктор-корректор имеет второй порядок аппроксимации.

Решение φ трижды непрерывно дифференцируемо (см. утверждение о гладкости решений) и поэтому

φ(ti) = φ(ti–1 + τ) = φ(ti–1) + τφ′(ti–1) +  τ2
2

φ′′(ti–1) + O3).

Ниже мы обозначаем f(t, x)/∂t через ft, (∂f(t, x)/∂x через fx, f(t, x) — через f, а значения этих функций в точке (ti, φ(ti)) — через ft, i, fx, i и fi, соответственно. В силу уравнения (E)

φ′(ti–1) = fi–1,

φ′′(ti–1) = ft, i–1 + fx, i–1φ′(ti–1) = ft, i–1 + fx, i–1fi–1.

Далее, поскольку fC2

f(t + τ, x + h) = f + ftτ + fxh + O2) + O(|| h|| 2).

Поэтому

f(ti, φ(ti–1) + τf[ti–1, φ(ti–1)]) =

= fi–1 + ft, i–1τ+ fx, i–1τfi–1 + O2) + O(|| τfi–1|| 2) =

= fi–1 + τφ′′(ti–1) + O2) + O(|| τfi–1|| 2).

Задача 1.3.1. Докажите, что O(||τfi–1||2) = O2) равномерно по i ∈ {1, ..., i}.

Поэтому, фактически, правая часть последнего равенства равна fi–1 + τφ′′(ti–1) + O2).

Оценим теперь ||[Fτ(Pτφ)]i||. При i = 0 очевидно [Fτ(Pτφ)]i = 0. Если же i > 0, то

[Fτ(Pτφ)]i =  (Pτφ)i – (Pτφ)i–1
τ
 – 1
2
( f[ti–1, (Pτφ)i–1] +

+ f(ti, (Pτφ)i–1 + τf[ti–1, φ(ti–1)]) =



=  
φ(ti–1) + τφ′(ti–1) +  τ2
2

φ′′(ti–1) + O3) –φ(ti–1)


τ


  –

1
2

[fi–1 + fi–1 + τφ′′(ti–1) + O2)] = φ′(ti–1) +

τ
2
φ′′(ti–1) +

O3)
τ
 – φ′(ti–1) –  τ
2

φ′′(ti–1) + O2) = O2).

Таким образом, ||[Fτ(Pτφ)]i|| ≤ Cτ2 при малых τ и при всех допустимых i.

Задача 1.3.2. Докажите, что эта оценка равеномерна по i ∈ {0, ..., i} и, следовательно, ||Fτ(Pτφ)||τ = O2).

Таким образом, схема предиктор-корректор является схемой второго порядка аппроксимации, и поэтому, если бы мы доказали ее устойчивость, то в силу теоремы Лакса он была бы схемой второго порядка точности. К этому вопросу мы еще вернемся.

1.3.4. Схемы Рунге — Кутты. Общие соображения.

Схема предиктор-корректор, так же, как и явная схема Эйлера, является представителем обширного семейства явных схем Рунге — Кутты. Их построение основано на следующей идее. Разложим φ в ряд Тейлора

φ(ti) = φ(ti–1) + τφ′(ti–1) +  τ2
2
φ′′(ti–1) + ..., 

оборвем ряд, взяв первые k + 1 его членов:
φ(ti) ≈ φ(ti–1) + τφ′(ti–1) + ...+  τk
k!

φ(k)(ti–1),

(2)

выразим все производные решения из уравнения (E) и подставим в правую часть (2):
φ(ti) ≈ φ(ti–1) + τδ(ti–1, φ(ti–1), τ);(3)

здесь через δ обозначен результат такой подстановки. Например, для k = 2
δ(t, x, τ) = f τ
2
(ft + fxf). 
(4)

Задача 1.3.3. Найдите δ при k = 4.

Разложение (3) порождает очевидным образом класс разностных схем
xi = xi–1 + τδ(ti–1, xi–1, τ). (5)

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

1.3.5. Схемы Рунге — Кутты. Пример.

Попытаемся аппроксимировать δ при k = 2 выражением вида

Φ(t, x, τ) = α1f(t, x) + α2f[t + β2τ, x + γ21τf(t, x)],

(выбор обозначений для индексов у коэффициентов α, β, γ станет ясен чуть ниже), подбирая коэффициенты α1, α2, β2, γ21 так, чтобы главные члены в разложениях δ и Φ по степеням τ были одинаковы. Имеем

Φ(t, x, τ) = (α1 + α2)f +τα22ft + γ21fxf) + O2).

Сравнение с (4) приводит к системе уравнений на коэффициенты α1, α2, β2, γ21

α1 + α2 = 1,   α2β2 = 1
2
,   α2γ21= 1
2
.

Эта система имеет однопараметрическое семейство решений

α1 = 1 – α,   α2 = α,   β2 = γ21 =  1
,

Задача 1.3.4. Докажите!

Найденный набор решений порождает однопараметрическое семейство разностных схем
 xi = xi–1 + τ ( (1 – α)f(ti–1, xi–1) + αf [ ti–1 + τ
, xi–1 + τ
f(ti–1, xi–1) ]).
(6)

Схема (6), естественно, дополняется начальным условием

x0 = x0,

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

При α = 1/2 мы получаем рассмотренную выше схему предиктор-корректор.

Задача 1.3.5. Каков порядок аппроксимации схемы (6) при различных α?

Схему (6) обычно записывают в виде
xi = xi–1 + τΦ(ti–1, xi–1, τ), (7)

где
Φ(t, x, τ) = α1k1 + α2k2 =  2

s = 1
αsks, 
(8)

k1 = f(t, x), k2 = f[t + β2τ, x + γ21τf(t, x)].(9)

и называют явной двухэтапной (или двухстадийной) схемой Рунге — Кутты по числу слагаемых в представлении (8) для функции Φ (или, что то же, числу вычислений правой части уравнения).

1.3.6. Явные схемы Рунге — Кутты.

Общая явная p-этапная схема Рунге — Кутты по определению имеет вид
xi = xi–1 + τΦ(ti–1, xi–1, τ),   x0 = x0,(10)

где
Φ(t, x, τ) = p

s = 1
αsks, 
(11)
 k1 = f(t, x),   ks = f ( t + βsτ, x + τ s–1

r = 1
γsr kr) ,   s = 2, ..., p.
(12)

Коэффициенты αs, βs и γsr определяются (как и в предыдущем пункте) так, чтобы функция Φ наилучшим образом аппроксимировала функцию δ в (5). Подробнее эта процедура выглядит так. Вычисляются частные производные функции Φ порядков 0, ..., p – 1 по τ при τ = 0 и приравниваются к производным точного решения. При этом для методов высокого порядка (p ≥ 3) обычно предполагаются выполненными дополнительные условия вида
βs =  s–1

r = 1
γsr, 
(13)

которые сильно упрощают как решение, так и исследование системы уравнений на коэффициенты искомых схем.

1.3.7. Уравнения на коэффициенты явной трехэтапной схемы.

Для примера изложим в виде задач план вывода уравнений на коэффициенты явной трехэтапной схемы Рунге — Кутты.

Задача 1.3.6. Покажите, что при k = 3

δ(t, x, τ) = f τ
2
(ft + fx f) +  τ2
6
(ftt + 2ftx f + fxx f 2 + fx ft + fx2f).

(здесь и ниже индексы у f обозначают соответствующие производные, например, ftx = f 2(t, x)/∂tx.

Задача 1.3.7. Покажите, что при p = 3

Φ(t, x, 0) = (α1 + α2 + α3)f,

Φ′τ(t, x, 0)= (α2β2 + α3β3)ft + (α2γ21 + α3γ31 + α3β3γ32)fx f,

Φ′′ττ(t, x, 0)= (α2β22+ α3β32)ftt+ 2(α2β2γ21 + α3β3γ31 + α3β3γ32)ftx f +

+ [α2γ221+ α331 + γ32)2]fxx f2 + 2α3β2γ32ft fx + 2α3γ21γ32ft2f.

Задача 1.3.8. Приравнивая δ, δ′τ,δ′′ττ и Φ, Φ′τ,Φ′′ττ,соответственно, в точке (t, x, 0), а затем приравнивая коэффициенты при соотвествующих агрегатах из производных функции f (например, при ftx f) и учитывая соотношения (13), покажите, что коэффициенты этой явной трехэтапной схемы Рунге — Кутты удовлетворяют системе уравнений

α1 + α2 + α3 = 1,
α2β2 + α3β3 = 1
2
,
α2β22+ α3β32= 1
6
,
α3β2γ32 = 1
3
,
γ21 = β2,
γ31 + γ32 = β3.
(14)

Задача 1.3.9. Найдите хотя бы одно решение системы (14). Найдите все ее решения.

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

1.3.8. Порядок аппроксимации явных схем Рунге — Кутты.

Система уравнений для определения коэффициентов схемы явной p-этапной схемы Рунге — Кутты в общем случае имеет семейство решений. Интересен вопрос о том, какая из соответствующих схем имеет наивысший порядок аппроксимации. (Здесь мы говорим о порядке аппроксимации, хотя имеем в виду, как это будет показано в следующем параграфе, порядок сходимости.) Полный ответ на этот вопрос не известен. Известно точное минимальное число этапов p(k), необходимое для достижения порядка аппроксимации k явной схемы Рунге — Кутты для всех k ≤ 7:

p(k) =  { k при k = 1,2,3,4,

k + 1 при k = 5,6,

k + 2 при k = 7.

Максимальный достигнутый порядок аппроксимации для явных схем Рунге — Кутты равен 10. Такой порядок достигается на построенной в 1975 г. явной восемнадцатиэтапной схеме. Эта схема занесена в книгу рекордов Гиннесса. Позднее построена семнадцатиэтапная схема 10-го порядка.

1.3.9. Неявные методы Рунге — Кутты.

Прямым обобщением рассмотренных в этом параграфе методов являются так называемые неявные p-этапные методы Рунге — Кутты. Они определяются следующими формулами (ср. с (10)(12))

xi = xi–1 + τΦ(ti–1, xi–1, τ),   x0 = x0,

где

Φ(t, x, τ) = p

s = 1
αsks, 
ks = f ( t + βsτ, x + τ p

r = 1
γsr kr) ,   s = 1, ..., p.
(15)

Коэффициенты αs, βs и γsr находятся из тех же соображений, что и для явных методов. Простейшим представителем этого класса схем является неявный метод Эйлера:

xi = xi–1 + τf(ti, xi).

В отличие от явных методов Рунге — Кутты формулы (15) не позволяют вычислять (m-мерные) векторы ks "один за другим": они (формулы) представляют собой систему из p m-мерных уравнений для p m-мерных неизвестных k1, ..., kp, или, что то же, систему pm скалярных уравнений с pm неизвестными.

Эта система при достаточно малых τ всегда однозначно разрешима. Доказать это утверждение можно так. Положим k = (k1, ..., kp) ∈ Rmp и F: RmpRmp (= (Rm)p) формулой

F(k) = (f1(k), ..., fp(k)),

где

fs(k) =   ( t + βsτ, x + τ p

r = 1
γsr kr) ,   s = 1, ..., p.

В этих обозначениях система (15) записывается в виде
k = F(k).(16)

Если теперь определить норму в Rmp, например, равенством ||k||mp = max1≤sp||ks||, то относительно этой нормы, как легко видеть, оператор F удовлетворяет условию Липшица с константой l = τL·max1≤sppr=1sr|.

Задача 1.3.10. Докажите!

Поэтому при l < 1 при достаточно малых τ, т. е. при таких τ оператор F является сжимающим. Но тогда в силу принципа сжимающих отображений (эквивалентное (15)) уравнение (16) имеет единственное решение, которое может быть приближенно вычислено методом простой итерации.

Задача 1.3.11. Проведите подробное доказательство.

Необходимость решать на каждом шаге систему (15) резко увеличивает объем необходимой вычислительной работы.

Если в формулах (15) γsr = 0 при s > r и хотя бы одно γss ≠ 0, то метод называется полуявным p-этапным методом Рунге — Кутты. В случае полуявных методов система уравнений (15) распадается на систему p (m-мерных) уравнений, которые можно решать поочередно: сначала уравнение

k1 = f(t + β1τ, x + τγ11k1)

относительно k1, затем уравнение

k2 = f(t + β2τ, x + τγ21k1+ τγ22k2)

(с уже найденным k1) относительно k2, и  т. д.

Задача 1.3.12. Найдите все неявные двухэтапные методы Рунге — Кутты.

Возрастающий объем вычислительных затрат для неявных методов Рунге — Кутты частично компенсируется бóльшим в общем случае порядком сходимости. Например, доказано, что для любого pN+ существует неявный p-этапный метод Рунге — Кутты 2p-го порядка сходимости. Кроме того, неявные методы по сравнению с явными обладают лучшими свойствами устойчивости (об этом мы будем говорить позже).


File based on translation from TEX by TTH, version 3.05.
Created 23 May 2002, 20: 16.