научная статья по теме ПРОСТОЙ СПОСОБ ПОСТРОЕНИЯ ДВУХШАГОВЫХ МЕТОДОВ РУНГЕ–КУТТЫ Математика

Текст научной статьи на тему «ПРОСТОЙ СПОСОБ ПОСТРОЕНИЯ ДВУХШАГОВЫХ МЕТОДОВ РУНГЕ–КУТТЫ»

ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 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 рублей.

Показать целиком