§ 1.3. Методы Рунге

§ 1.3. Методы Рунге — Кутты

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

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

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

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

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

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

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

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

Попытаемся повысить порядок аппроксимации явной схемы Эйлера следующим образом. В исходном уравнении (E) заменим производную обычным конечно-разностным соотношением, а f(t, x) постараемся аппроксимировать "с более высоким порядком". Например, возьмем вместо f(t_i–1, x_i–1) "среднее направление" между векторами поля направлений уравнения (E) в точках (t_i–1, x_i–1) (вектор f₁ на рис. 3) и (t_i, x*_i), где (t_i, x*_i)— следующая вершина ломаной Эйлера (вектор f₂). Интуитивно ясно, что получившееся направление (вектор f₃) ближе к "истинно нужному" направлению (вектор f₄).

Рис. 3.

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

[F_τ(x)]_i = {
x_i – x_i–1
τ   –   1
2 (
f(t_i–1, x_i–1) + f[t_i, x_i–1 + τf(t_i–1,x_i–1)]

) ,  если i = 1, ..., n,
x_i –x⁰,  если i = 0;

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

(φ_τ)₀ = x⁰,

(φ_τ)_i = (φ_τ)_i–1 + τ
2 ( f[t_i–1, (φ_τ)_i–1] + f(t_i, (φ_τ)_i–1 + f[t_i–1, (φ_τ)_i–1]) ) .

Разностную схему Хойна (вернее рекуррентный переход от x_i–1 к x_i) часто записывают в виде двух полушагов:
x*_i= x_i–1 + τf(t_i–1, x_i–1), (1а)

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

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

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

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

φ(t_i) = φ(t_i–1 + τ) = φ(t_i–1) + τφ′(t_i–1) + τ²
2
φ′′(t_i–1) + O(τ³).

Ниже мы обозначаем ∂f(t, x)/∂t через f_t, (∂f(t, x)/∂x — через f_x, f(t, x) — через f, а значения этих функций в точке (t_i, φ(t_i)) — через f_{t, i}, f_{x, i} и f_i, соответственно. В силу уравнения (E)

φ′(t_i–1) = f_i–1,

φ′′(t_i–1) = f_{t, i–1} + f_{x, i–1}φ′(t_i–1) = f_{t, i–1} + f_{x, i–1}f_i–1.

Далее, поскольку f ∈ C²

f(t + τ, x + h) = f + f_tτ + f_xh + O(τ²) + O(|| h|| ²).

Поэтому

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

= f_i–1 + f_{t, i–1}τ+ f_{x, i–1}τf_i–1 + O(τ²) + O(|| τf_i–1|| ²) =

= f_i–1 + τφ′′(t_i–1) + O(τ²) + O(|| τf_i–1|| ²).

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

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

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

[F_τ(P_τφ)]_i = (P_τφ)_i – (P_τφ)_i–1
τ – 1
2 ( f[t_i–1, (P_τφ)_i–1]+

+ f(t_i, (P_τφ)_i–1 + τf[t_i–1, φ(t_i–1)]) =

=
φ(t_i–1) + τφ′(t_i–1) + τ²
2
φ′′(t_i–1) + O(τ³) –φ(t_i–1)

τ

–

– 1
2
[f_i–1 + f_i–1 + τφ′′(t_i–1) + O(τ²)] = φ′(t_i–1) +

τ
2 φ′′(t_i–1)+

+ O(τ³)
τ – φ′(t_i–1) – τ
2
φ′′(t_i–1) + O(τ²) = O(τ²).

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

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

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

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

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

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

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

(2)

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

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

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

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

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

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

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

Φ(t, x, τ) = α₁f(t, x) + α₂f[t + β₂τ, x + γ₂₁τf(t, x)],

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

Φ(t, x, τ) = (α₁ + α₂)f +τα₂(β₂f_t + γ₂₁f_xf) + O(τ²).

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

α₁ + α₂ = 1, α₂β₂ = 1
2 , α₂γ₂₁= 1
2 .

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

α₁ = 1 – α, α₂ = α, β₂ = γ₂₁ = 1
2α ,

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

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

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

x₀ = x⁰,

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

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

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

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

где
Φ(t, x, τ) = α₁k₁ + α₂k₂ = 2
∑
s = 1
α_sk_s,
(8)

k₁ = f(t, x), k₂ = f[t + β₂τ, x + γ₂₁τf(t, x)]. (9)

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

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

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

где
Φ(t, x, τ) = p
∑
s = 1
α_sk_s,
(11)

k₁ = f(t, x), k_s = f ( t + β_sτ, x +τ s–1
∑
r = 1
γ_srk_r ) , 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 (f_t + f_x f) + τ²
6 (f_tt + 2f_tx f + f_xx f ² + f_x f_t + f_x²f).

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

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

Φ(t, x, 0) = (α₁ + α₂ + α₃)f,

Φ′_τ(t, x, 0)= (α₂β₂ + α₃β₃)f_t + (α₂γ₂₁ + α₃γ₃₁ + α₃β₃γ₃₂)f_x f,

Φ′′_ττ(t, x, 0)= (α₂β₂²+ α₃β₃²)f_tt+ 2(α₂β₂γ₂₁ + α₃β₃γ₃₁ + α₃β₃γ₃₂)f_tx f +

+ [α₂γ²₂₁+ α₃(γ₃₁ + γ₃₂)²]f_xx f₂ + 2α₃β₂γ₃₂f_t f_x + 2α₃γ₂₁γ₃₂f_t²f.

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

α₁ + α₂ + α₃ = 1,
α₂β₂ + α₃β₃ = 1
2 ,

α₂β₂²+ α₃β₃²= 1
6 ,

α₃β₂γ₃₂ = 1
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))

x_i = x_i–1 + τΦ(t_i–1, x_i–1, τ), x₀ = x⁰,

где

Φ(t, x, τ) = p
∑
s = 1
α_sk_s,

k_s =f ( t + β_sτ, x +τ p
∑
r = 1
γ_srk_r ) , s = 1, ..., p.
(15)

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

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

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

Эта система при достаточно малых τ всегда однозначно разрешима. Доказать это утверждение можно так. Положим k = (k₁, ..., k_p) ∈ R^mp и F: R^mp→R^mp (= (R^m)^p) формулой

F(k) = (f₁(k), ..., f_p(k)),

где

f_s(k) = ( t + β_sτ, x +τ p
∑
r = 1
γ_srk_r ) , s = 1, ..., p.

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

Если теперь определить норму в R^mp, например, равенством ||k||_mp = max_1≤s≤p||k_s||, то относительно этой нормы, как легко видеть, оператор F удовлетворяет условию Липшица с константой l = τL·max_1≤s≤p∑^p_r=1|γ_sr|.

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

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

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

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

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

k₁ = f(t + β₁τ, x + τγ₁₁k₁)

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

k₂ = f(t + β₂τ, x + τγ₂₁k₁+ τγ₂₂k₂)

(с уже найденным k₁) относительно k₂, и т. д.

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

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

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