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

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

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

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

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

В этом и последующих параграфах в случаях, когда это необходимо, мы считаем уравнение (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) с разностным оператором

[Ft(x)]i = м п н п о xi – xi–1 t   –   1 2 ж и f(ti-1, xi-1) + f[ti, xi-1 + tf(ti-1, xi-1)] ц ш ,  если i = 1, ..., n, xi –x0,  если i = 0;

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

(jt)0 = x0,

(jt)i = (jt)i-1 +  t 2 ж и f[ti-1, (jt)i-1] +  f(ti, (jt)i-1 + f[ti-1, (jt)i-1]) ц ш .

Разностную схему Хойна (вернее рекуррентный переход от x_i-1 к x_i) часто записывают в виде двух полушагов:

x*i= xi-1 + tf(ti-1, xi-1),

(1а)

xi = xi-1 +  t 2 [f(ti-1, xi-1) + f(ti, x*i)].

(1б)

Поэтому она обычно называется схемой предиктор-корректор: на первом полушаге (1а) приближенное решение "предсказывается" (от англ. to predict — предсказывать) с первым порядком точности, а на втором (1б) — "корректируется" (от to correct — исправлять, корректировать) с целью повышения точности.

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

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

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

j(ti) = j(ti-1 + t) = j(ti-1) + tjў(ti-1) +  t2 2 jўў(ti-1) + O(t3).

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

jў(ti-1) = fi-1,

jўў(ti-1) = ft, i-1 + fx, i-1jў(ti-1) = ft, i-1 + fx, i-1fi-1.

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

f(t + t, x + h) = f + ftt + fxh + O(t2) + O(|| h|| 2).

Поэтому

f(ti, j(ti-1) + tf[ti-1, j(ti-1)]) =

= fi-1 + ft, i-1t+ fx, i-1tfi-1 + O(t2) + O(|| tfi-1|| 2) =

= fi-1 + tjўў(ti-1) + O(t2) + O(|| tfi-1|| 2).

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

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

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

[Ft(Ptj)]i =  (Ptj)i - (Ptj)i-1 t  –  1 2 ж и f[ti-1, (Ptj)i-1] +

+ f(ti, (Ptj)i-1 + tf[ti-1, j(ti-1)]) =

=   j(ti-1) + tjў(ti-1) +  t2 2 jўў(ti-1) + O(t3) -j(ti-1) t   –

– 1 2 [fi-1 + fi-1 + tjўў(ti-1) + O(t2)] = jў(ti-1) + t 2 jўў(ti-1) +

+  O(t3) t  - jў(ti-1) -  t 2 jўў(ti-1) + O(t2) = O(t2).

Таким образом, ||[F_t(P_tj)]_i|| Ј Ct² при малых t и при всех допустимых i.

Задача 1.3.2. Докажите, что эта оценка равеномерна по i О {0, ..., i} и, следовательно, ||F_t(P_tj)||_t = O(t²).

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

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

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

j(ti) = j(ti-1) + tjў(ti-1) +  t2 2 jўў(ti-1) + ...,

оборвем ряд, взяв первые k + 1 его членов:

j(ti) » j(ti-1) + tjў(ti-1) + ...+  tk k! j(k)(ti-1),

(2)

выразим все производные решения из уравнения (E) и подставим в правую часть (2):

j(ti) » j(ti-1) + td(ti-1, j(ti-1), t);

(3)

здесь через d обозначен результат такой подстановки. Например, для k = 2

d(t, x, t) = f +  t 2 (ft + fxf).

(4)

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

Разложение (3) порождает очевидным образом класс разностных схем

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

(5)

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

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

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

F(t, x, t) = a1f(t, x) + a2f[t + b2t, x + g21tf(t, x)],

(выбор обозначений для индексов у коэффициентов a, b, g станет ясен чуть ниже), подбирая коэффициенты a₁, a₂, b₂, g₂₁ так, чтобы главные члены в разложениях d и F по степеням t были одинаковы. Имеем

F(t, x, t) = (a1 + a2)f +ta2(b2ft + g21fxf) + O(t2).

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

a1 + a2 = 1,   a2b2 =  1 2 ,   a2g21=  1 2 .

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

a1 = 1 - a,   a2 = a,   b2 = g21 =  1 2a ,

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

Найденный набор решений порождает однопараметрическое семейство разностных схем

xi = xi-1 + t ж и (1 - a)f(ti-1, xi-1) + af й л ti-1 +  t 2a , xi-1 +  t 2a f(ti-1, xi-1) щ ы ц ш .

(6)

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

x0 = x0,

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

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

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

Схему (6) обычно записывают в виде

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

(7)

где

F(t, x, t) = a1k1 + a2k2 =  2 е s = 1 asks,

(8)

k1 = f(t, x), k2 = f[t + b2t, x + g21tf(t, x)].

(9)

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

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

Общая явная p-этапная схема Рунге — Кутты по определению имеет вид

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

(10)

где

F(t, x, t) =  p е s = 1 asks,

(11)

k1 = f(t, x),   ks = f ж и t + bst, x + t s–1 е r = 1 gsr kr ц ш ,   s = 2, ..., p.

(12)

Коэффициенты a_s, b_s и g_sr определяются (как и в предыдущем пункте) так, чтобы функция F наилучшим образом аппроксимировала функцию d в (5). Подробнее эта процедура выглядит так. Вычисляются частные производные функции F порядков 0, ..., p - 1 по t при t = 0 и приравниваются к производным точного решения. При этом для методов высокого порядка (p і 3) обычно предполагаются выполненными дополнительные условия вида

bs =  s–1 е r = 1 gsr,

(13)

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

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

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

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

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

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

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

F(t, x, 0) = (a1 + a2 + a3)f,

Fўt(t, x, 0)= (a2b2 + a3b3)ft + (a2g21 + a3g31 + a3b3g32)fx f,

Fўўtt(t, x, 0)= (a2b22+ a3b32)ftt+ 2(a2b2g21 + a3b3g31 + a3b3g32)ftx f +

+ [a2g221+ a3(g31 + g32)2]fxx f2 + 2a3b2g32ft fx + 2a3g21g32ft2f.

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

a1 + a2 + a3 = 1, a2b2 + a3b3 =  1 2 , a2b22+ a3b32=  1 6 , a3b2g32 =  1 3 , g21 = b2, g31 + g32 = b3.

(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 + tF(ti-1, xi-1, t),   x0 = x0,

где

F(t, x, t) =  p е s = 1 asks,

ks = f ж и t + bst, x + t p е r = 1 gsr kr ц ш ,   s = 1, ..., p.

(15)

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

xi = xi-1 + tf(ti, xi).

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

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

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

где

fs(k) =   ж и t + bst, x + t p е r = 1 gsr kr ц ш ,   s = 1, ..., p.

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

k = F(k).

(16)

Если теперь определить норму в Rmp, например, равенством ||k||mp = max1ЈsЈp||ks||, то относительно этой нормы, как легко видеть, оператор F удовлетворяет условию Липшица с константой l = tL·max1ЈsЈpеpr=1|gsr|.

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

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

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

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

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

k1 = f(t + b1t, x + tg11k1)

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

k2 = f(t + b2t, x + tg21k1+ tg22k2)

(с уже найденным 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.