ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2009, том 49, № 11, с. 1920-1930
УДК 519.622
ПРОСТОЙ СПОСОБ ПОСТРОЕНИЯ ДВУХШАГОВЫХ МЕТОДОВ РУНГЕ-КУТТЫ
© 2009 г. Л. М. Скворцов
(105005Москва, ул. 2-я Бауманская, 5, МГТУим. Н.Э. Баумана) e-mail: lm_skvo@rambler.ru Поступила в редакцию 12.12.2008 г. Переработанный вариант 04.05.2009 г.
Предложен способ построения двухшаговых методов Рунге—Кутты на основе одношаговых методов. Рассматриваются явные и диагонально неявные двухшаговые методы, имеющие второй или третий стадийный порядок. На тестовых задачах показано преимущество предложенных методов в сравнении с обычными одношаговыми методами. Библ. 15. Табл. 5.
Ключевые слова: двухшаговые методы Рунге—Кутты, стадийный порядок, явные методы, диагонально-неявные методы, жесткие системы уравнений.
1. ВВЕДЕНИЕ
Двухшаговые методы Рунге—Кутты (TSRK, см. [1]—[7]) обобщают обычные одношаговые методы, используя на очередном шаге интегрирования информацию, полученную не только на текущем, но и на предыдущем шаге. Благодаря этому двухшаговые методы могут иметь более высокий стадийный порядок, что позволяет повысить их точность, особенно при решении жестких и дифференциально-алгебраических задач. Общий класс двухшаговых методов Рунге—Кутты предложен в [1]. Мы рассмотрим один частный подкласс этого класса, отличающийся простотой построения таких методов и простотой их программной реализации.
Будем решать задачу Коши для системы обыкновенных дифференциальных уравнений (ОДУ)
У' = f(t, У), У(to) = У0, (1.1)
где y — вектор переменных, f — векторная функция, t — независимая переменная. Рассмотрим двухшаговые методы Рунге—Кутты, задаваемые формулами
уП = y„, Fl = f(t„, У), (1.2а)
i +1
Yl = Уп + hn£ j_ i + bF), Fl = f(tn + cthm Y), i = 2, 3, i, (1.2б)
j = i
i +1
Уi +1 = Уп + hn£ bs + iX, Fl +1 = f(tn + hn, У +1), Уп +1 = У +1, (1.2в)
j = 1
где hn — размер очередного шага, а Y'n и Fn , i = 1, 2, ..., s + 1, — стадийные значения и их производные на этом шаге. Будем называть стадии (1.2б) внутренними, а стадию (1.2в) заключительной. Отметим, что стадия (1.2а) не требует вычислений, поскольку результат ее выполнения совпадает с результатом выполнения заключительной стадии предыдущего шага. Поэтому говорят, что такие методы обладают свойством FSAL (First Same As Last). Представим коэффициенты метода (1.2) в виде двух матриц и вектора
A = [ ajj], B = [ bij], c = [ ci], i, j = 1, 2, ..., i + 1,
1920
при этом примем с1 = 0, с5 + 1 = 1, аь = а5 +1; I = Ьц = 0, I = 1, 2, ..., 5 + 1.
В предлагаемом способе построения двухшаговых методов вида (1.2) за основу берется обычный одношаговый метод Рунге-Кутты, задаваемый матрицей В и вектором с. Назовем его исходным методом. Матрицу А примем в виде
А = dgт, (1.3)
где векторы d и g определяются исходя из условий повышения стадийного порядка исходного метода и обеспечения необходимых свойств устойчивости. Такой подход позволяет повысить стадийный порядок исходного метода на 1 или на 2. В общем случае все коэффициенты двухша-говых методов зависят от соотношения размеров шагов w = к„/к„ _ 1. В наших методах от w зависят только коэффициенты вектора g, что упрощает их реализацию с переменным шагом. Первый шаг выполняется исходным одношаговым методом, поэтому предложенные методы не нуждаются в специальной стартовой процедуре.
С учетом (1.3) формулы (1.2) можно записать в виде, более удобном для их реализации:
* +1
V = У„, = уП), и = /Г-1,
1 = 1
С * +1 Л
VI = Уп + К
и
\и + XИ , К = ЦК + фп, VI), I = 2, 3, ..., *, (1.4)
4 1= 1 у
* +1
у*+1 = У + нТь , г г*+1 = ха + К У*+1) У , = У*+1
1 = 1
2. УСЛОВИЯ ПОРЯДКА
Вывод условий порядка основывается на сравнении разложений в ряд Тейлора точного и численного решений. Часто используют предположения, позволяющие упростить условия порядка. При этом важным понятием является стадийный порядок, определяемый как наименьший порядок на всех стадиях. Условия порядка двухшаговых методов Рунге-Кутты рассматривались в [1]_[3]. Основываясь на этих работах, приведем некоторые условия порядка для методов вида (1.2) и (1.4).
Примем следующие обозначения (5 + 1)-мерных векторов:
е = [ 1.....1 ]т, е* +1 = [0, ..., 0, 1 ]т, Ь = [Ь* +1, ь Ь* +1,* +1 ]т.
Предположим, что стадийный порядок исходного метода равен д , т.е. выполняются условия
кВск-1 = с, к = 1, 2, ..., д, (д + 1)Всд Ф сд +1 (2.1)
(здесь и далее предполагается покомпонентное выполнение операции возведения вектора в степень). Примем также обычное для методов Рунге-Кутты предположение, что д > 1, тогда с = Ве и метод однозначно задается матрицей В. Пусть р - порядок исходного метода, причем будем предполагать, что р > д .
Рассмотрим двухшаговый метод (1.2), и пусть он имеет стадийный порядок Тогда должны выполняться равенства
к
А|с—е)к -1 + Вск -1 к
= ск, к = 1, 2, ..., д. (2.2)
Предполагая, что q > д, и подставляя (1.3) в (2.2), при к = д + 1 получаем
(д + 1)dg
т* с - е
№
= сд +1 - (д + 1)Всд.
Это равенство будет выполняться, если принять
Из (2.1), (2.2) имеем
а из (2.4), (2.5) получим
й = сд +1 - (д + 1)Всд,
g (с - е) =
д + 1
gтc = о, к = о, 1.....д-1,
(2.3)
(2.4)
(2.5)
т д №
g с =
д + 1
(2.6)
Равенства (2.3), (2.5), (2.6) будем использовать как условия, обеспечивающие стадийный порядок q = д + 1 метода (1.4). Для обеспечения стадийного порядка q = д + 2 дополнительно должны выполняться условия
сд +2 - (д + 2 )Всд +1 = а( сд +1 - (д + 1) Всд),
д+1
д +1
д +1
т д +1 № д
g с = а =—- + № , д + 2
(2.7)
(2.8)
где а — некоторая константа.
Если q — стадийный порядок метода (1.2), то его порядок (не ниже) р = q. Метод (1.2) имеет порядок р = q + 1, если
(д + 1) Ъсд = 1.
Метод (1.2) имеет порядок р = q + 2, если выполняются условия (2.9) и
(д + 2) Ъсд +1 = 1,
( д + 2)( д + 1) Ът
с - е
№
+ Всд
= 1.
(2.9)
(2.10) (2.11)
Если р > д + 1, то (2.11) можно упростить. В этом случае из условий порядка одношагового метода имеем
Ъй = ъ т( сд+1 - (д + 1) Всд) = о.
Подставляя это выражение и (1.3) в (2.11), получаем
(д + 2)(д + 1) Ът Всд = 1. (2.12)
Таким образом, метод (1.4) имеет порядокр = q + 2, если р > д + 1 и выполняются условия (2.9), (2.10), (2.12).
д
д
3. УСТОЙЧИВОСТЬ
Устойчивость методов решения ОДУ принято исследовать на модельном уравнении Далкви-ста у' = Ху, у(?0) = у0. Используя формулы (1.2) для решения этого уравнения, получаем
V = (ее*т+1 + гА)У _ 1 + гВУп, г = НЛХ. (3.1)
Разрешая (3.1) относительно Уи и учитывая (1.3), имеем
Уп = Н (г) V -1, Н(г) = (I - г В)-1 (ее*т+1 + ), (3.2)
где I - единичная матрица. Устойчивость разностного уравнения (3.2) определяется спектром матрицы Н(г), т.е. корнями характеристического многочлена Р(п, г) = |п1 - Н(г)|.
В общем случае Р(п, г), как многочлен от п, имеет 5 + 1 различных корней, зависящих от г, но при выборе матрицы А в виде (1.3) ранг матрицы Н(г) равен 2, а в этом случае
Р(п, г) = п*- 1[п2 -Р1 (г)п + Ро(г)]. (3.3)
Таким образом, устойчивость предлагаемых методов определяется двумя нулями
П1,2 (г) = 1 (Р1 (г) ±Л/р2( г) - 4ро (г)) (3.4)
многочлена (3.3). Согласно определению У.1.1 из [8], область устойчивости метода задается неравенствами
(г )|< 1, | п 2 (г )|< 1, (3.5)
при этом оба нуля не должны одновременно равняться 1 или -1.
Чтобы получить выражения дляр0(г), Р1(г), представим матрицу Н(г) в виде
Н (г) = ХУт, X = (I - гВ)-1[ е г d ], У = [ е * +1 g ].
Воспользовавшись тем, что матрицы Н(г) и У(г) = Ут X имеют одинаковые ненулевые собственные значения, получим
Р1 (г) = V!!(г) + V22(г), Ро(г) = Vll(г^(г) - Vl2(г)v2l(г),
^(г) = е*т+1 (I - гВ )-1е, Vl2(г) = е*т+1 (I - гВ )-1^, (3.6)
V2l (г) = ^ (I - г В )-1е, V22 (г) = g т( I - гВ)-1 г d.
Заметим, что V11(z) является функцией устойчивости исходного одношагового метода. Соотношения (3.4)-(3.6) использовались нами для построения областей устойчивости двухшаговых методов.
4. ЯВНЫЕ МЕТОДЫ
Явные методы имеют Ьу = 0 при] > I. В этом случае исходный метод имеет первый стадийный порядок д = 1 и может быть представлен в виде таблицы Бутчера
0
С2 Ь21
с* Ь*1 . ■ Ь*, * - 1
Ь* + 1,1 • ■ Ь* + 1, * - 1 Ь* + 1, *
где ж — число стадий, равное числу вычислений правой части на одном шаге интегрирования. В соответствии с (2.3) принимаем
й = с2 - 2Вс, (4.1)
а из (2.5), (2.6) получаем
gтe = 0, gтc = w /2.
(4.2)
Двухшаговые методы 2-го стадийного порядка можно построить на основе любого одношаго-вого метода, порядок которого не ниже 2-го. Для этого достаточно задать d по формуле (4.1), а g определить из (4.2). При решении уравнений (4.2) используются только два коэффициента вектора g, остальные можно задать нулевыми либо выбрать из других соображений. Например, приравняв нулю коэффициенты g2, ..., gs, получим
g = |[-1, 0, ..., 0, 1 ]'.
(4.3)
Формулы (4.1), (4.3) задают наиболее простой способ построения явных двухшаговых методов.
Приведем конкретный метод. При ж = 2 оптимальный по точности одношаговый метод 2-го порядка получаем, задав с2 = 2/3. Построенный на его основе двухшаговый метод задается коэффициентами
В=
0 0 0 2/3 0 0 , с = 0 2/3 , й = 0 4/9 w , g = 22 -1 0
1/4 3/4 0 1 _ 0_ _1_
и имеет порядок р = 3 и стадийный порядок # = 2. Длина его области устойчивости вдоль вещественной оси (при постоянном шаге) I = 2.0. Обозначим этот метод через Т8ЯК23 (первая цифра ж, вторая р).
Двухшаговые методы 3-го стадийного порядка можно построить на основе одношаговых методов, коэффициенты которых удовлетворяют соотношению (2.7) при а = с2. В этом случае из (2.8) получаем
g с = 1 + —^. (4.4)
Вектор d определяем согласно (4.1), а g находим из (4.2), (4.4). При ж > 2 число коэффициентов вектора g больше, чем необходимо для решения уравнений (4.2), (4.4). Поэтому мы задавали ненулевыми только три коэффициента этого вектора: g1, gk и gs +1, а к выбирали так, чтобы область устойчивости была как можно больше. Решение уравнений (4.2), (4.4) относительно выбранных ненулевых коэффициентов имеет вид
(3 + 2 мс2 Л (3 + 2 ) (3 + 2 Л
81 = 71--3!, 8к = т—т-77", 8 +1 = 71 —.-+ 3!.
6 4 ск ; 6 Ск( Ск-1) 6 V 1 - Ск ;
Примем ж = 3 и пос
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.