Явные методы второго порядка с согласованными областями устойчивости Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 6, № 4, 2001

ЯВНЫЕ МЕТОДЫ ВТОРОГО ПОРЯДКА С СОГЛАСОВАННЫМИ ОБЛАСТЯМИ УСТОЙЧИВОСТИ

Е.А. Новиков, Л.Н. Контарева Институт вычислительного моделирования СО РАН

Красноярск, Россия e-mail: [email protected]

For arbitrary m the coefficients of explicit m — stage methods of Runge — Kutta type are obtained. For these methods the stability domains of intermediate numerical schemes are conformable with the stability domain of the basic scheme.

Введение

Для численного решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений

у' = f (t,y), y(to) = Уо, ta < t < tk (1)

в [1] предлагается применять явные методы типа Рунге — Кутты

m / i—1 \

Vn+1 = Vn + Pmiki, ki = hf I tn + aih, yn + ^ fiijkj I , (2)

i=i V j=i J

где y и f — гладкие вещественные N-мерные вектор-функции; t — независимая переменная; h — шаг интегрирования; ki (1 < i < m) — стадии метода; pmi, ai, fîj (1 < i < m, 1 < j < i — 1) — коэффициенты, определяющие свойства точности и устойчивости схемы (2).

В [1] для произвольного m получены коэффициенты методов с первого по третий порядок, а в [2] — коэффициенты методов четвертого порядка с максимальным интервалом устойчивости. Там же численно показано значительное повышение эффективности алгоритмов интегрирования за счет комбинирования таких численных формул в процессе вычислений на основании критерия устойчивости. Отметим, что в [1, 2] не рассмотрен вопрос о выборе коэффициентов , которые влияют на устойчивость промежуточных схем и в конечном счете на эффективность алгоритма интегрирования.

Здесь для произвольного m получены коэффициенты ai, pmi, fîj (1 < i < m, 1 < j < i — 1 ) явных методов типа Рунге — Кутты второго порядка точности, при которых области устойчивости промежуточных численных схем согласованы с основной.

© Е.А. Новиков, Л.Н. Контарева, 2001.

1. Численные схемы

Для упрощения выкладок рассмотрим задачу Коши для автономной системы обыкновенных дифференциальных уравнений:

У = /(у), у(*о) = Уо, и < г < 4, (3)

для решения которой применим методы типа Рунге — Кутты, записанные в виде

Уп,г = Уп + вг+1,3к, 1 < г < т - 1, Уп+1 = Уп + Ртгкг 3=1 г=1

(4)

где кг = к/(Уп,г_1), 1 < г < т, Уп,0 = Уп. Все полученные ниже результаты можно использовать для неавтономных систем, если в (2) положить

г_1

аг = , 2 < г < т, а1 = 0.

3=1

(5)

В дальнейшем нам потребуется верхняя треугольная матрица Вт с элементами Ьу вида [1]

Ь1г = 1, 1 < г < т,

г_1

Ькг = ^^ Ьк_1,з, 2 < к < т, к < г < т, 3 = к_ 1

Ькг = 0, 2 < к < т, 1 < г < к - 1,

(6)

где вгз — коэффициенты схемы (2) или (4).

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

У' = Л У, у(0) = Уо, г > о с комплексным Л, КеЛ < 0. Применяя вторую формулу (4) к (7), получим

Уп+1 = Ят(г )Уп, г = Лк. Здесь Ят(г) есть многочлен степени т относительно г вида

(7)

Ят(г) = 1 + ^

стгг ,

г=1

где

У^ Ьг3Рт3 , 1 < г < т.

(8)

(9)

3=г

Для упрощения этой записи введем обозначения

рт (рт1, • • • , ртт) , Ст (ст1, • • • , Стт) •

Тогда коэффициенты схемы (4) и многочлена устойчивости (8) связаны соотношением

Втрт cm, (10)

т

С

где элементы матрицы Бт определены формулами (6). Для промежуточных численных схем (4) имеем

Уп,к = Як(г)уп, Як (г) = 1 + ^ с^/, (11)

i=1

где

Cki = ^ ь^вк+1,з, 1 < к < т - 1. (12)

Используя обозначения

Рк = (вк+1,1) • • • ) Рк+1,к)Т, ск = (ск1, • • • , скк)Т

получим, что коэффициенты в^ промежуточных схем (4) и соответствующих многочленов устойчивости связаны соотношениями

Бквк = ск, 1 < к < т - 1 (13)

Заметим, что из сравнения (6) и (12) следует, что Ь^ = ci-1)k-1, т.е. элементы (к + 1)-го столбца матрицы Бт совпадают с коэффициентами многочлена устойчивости Як (г). Отсюда следует, что если заданы коэффициенты многочленов устойчивости основной и промежуточных численных схем, то все коэффициенты методов (4) однозначно определяются из линейных систем (10) и (13). Ниже этот факт будет использоваться при определении коэффициентов формулы (4) второго порядка точности при условии, что области устойчивости промежуточных схем согласованы с основной.

Теперь запишем соотношения, обеспечивающие второй порядок точности схемы (4). Разлагая точное и приближенное решения в ряды Тейлора по степеням Л, можно записать

Л2 Л3

у(<п+1) = У(<п) + Л/ + у/'/ + - [/"/2 + Г2/] + 0(Л4),

т

Уп+1 = Уп + ( ^ЬцРтЛ Л/п + ( ^Ь2]Ртз\ /п+

кЗ=1 ) \3=2

т

+ ^ЗЬЬззРт^ //п + 0-5 ^3=2Ь2зРт^ // + ), (14)

где элементарные дифференциалы вычислены на точном у(^п) и приближенном уп решениях соответственно. Отсюда следует, что (4) будет иметь второй порядок точности, если

тт

^Т Ь13РтЗ = 1, У]Ь23РтЗ = 0А (15)

3=1 3=2

Это требование будет выполнено, если в (10) положить ст1 = 1 и ст2 = 0^5. Остальные коэффициенты с^ (3 < г < т) обычно используются для задания области устойчивости необходимых формы и размера. В [1] подробно изложен алгоритм вычисления нужных значений с^ (3 < г < т), и поэтому ниже будем считать их определенными. Теперь при

т

заданных коэффициентах (2 < г < т, 1 < ] < г — 1) решение линейной алгебраической системы (10) дает коэффициенты ртг (1 < г < т) т-стадийных методов (4) второго порядка точности с заданной областью устойчивости. Неравенство для контроля точности данных методов построим по аналогии с тем, как это сделано в [1]. Для этого необходимо, чтобы выражение для локальной ошибки 8пт т-стадийной численной формулы (4) имело вид

5п,т = + 0(Л4), (16)

где й — некоторая вычисленная постоянная. Из представления (14) следует, что требование (16) будет выполнено, если одновременно с (10) при ст1 = 1 и ст2 = 0.5 будет иметь место соотношение

т1

Ртз = х. (17)

3=2 3

Таким образом, задача о построении m-стадийных методов (4) второго порядка точности с заданной областью устойчивости и неравенствами для контроля точности вычислений сводится к совместному решению системы (10) и уравнения (17). Так как второе уравнение (10) совпадает со вторым уравнением (15), система (10), (17) будет совместна, если совместны вторые уравнения из (15) и (17), которые учитывая, что Ь2? = а? (2 < ] < т), можно переписать в виде

т т 1

X азРтз = X] аЗРтЗ = 3, (18)

3=2 3=2

где аг (2 < г < т) определены формулами (5).

2. Определение коэффициентов численных методов

Будем предполагать, что заданы коэффициенты многочленов устойчивости

к

Qk (г) = 1 + X скггг, 2 < к < т, (19)

г=1

причем ст1 = 1 и ст2 = 0.5. Коэффициент с11 полинома ^1(г) = 1 + с11 г выберем из условия совместности уравнений (18).

При выборе коэффициентов многочленов (19) можно руководствоваться различными соображениями. Если, например, основная численная схема (4) применяется с максимальным интервалом устойчивости, то имеет смысл выбрать полиномы устойчивости промежуточных схем с аналогичным свойством. Если же область устойчивости основной схемы "раздута" по мнимой оси, то области устойчивости промежуточных схем должны быть устроены подобным образом. Это позволяет избежать неустойчивости промежуточных численных формул, когда шаг интегрирования по устойчивости основной схемы близок к максимальному. Коэффициенты нужных многочленов можно вычислить предложенным в [1] алгоритмом. В [3] приведены многочлены Qk(г) (2 < к < 13) с максимальным интервалом устойчивости.

Для каждого к (2 < к < т) обозначим через 7к длину такого максимального интервала [7к, 0], что для всех г € [7к, 0] имеет место неравенство ^к(г)| < 1. Учитывая, что г = АЛ,

в (19) для каждого Qk(2 < к < т) проведем замену к на 7кк/7т. В результате вместо (19) рассмотрим набор многочленов

Q'k (*) = 1 + X 4*г, 4 = (7*/(2 < к < т).

(20)

г=1

Далее данные полиномы будут использоваться в качестве многочленов устойчивости методов (4). Замена к на 7*к/7т означает, что приближенное решение по промежуточным схемам (4) вместо точек (¿п + с*1к) (2 < к < т — 1) будет вычисляться в точках Ьп + ск1к (2 < к < т — 1). В этом случае максимальный шаг из условия устойчивости основной схемы будет максимальным и для промежуточных численных формул.

Теперь определение коэффициентов методов (4) будем осуществлять по следующему алгоритму. С использованием [1] вычислим коэффициенты полиномов (19), удовлетворяющие некоторым заданным свойствам. Вычислим коэффициенты многочленов (20) с применением соответствующей замены переменных. Учитывая, что элементы (к + 1)-го столбца матрицы Вт совпадают с коэффициентами многочлена устойчивости Q'k (¿), сформируем матрицу Вт, которая имеет вид

/ 1 1 1 ■ ■1 \

0 с'11 с21 ■ ■ с т- 1, 1

Вт = 0 0 с22 ■ ■ Ст—1,2

V 0 0 0 ■ ■ с т- 1 ,т- 1 )

где с'11 есть пока неопределенный коэффициент полинома QZ = 1 + с'11г. Из линейной системы уравнений (10) с верхней треугольной матрицей Вт последовательно определим Рт,т, Рт,т- 1, • • • , Рт3. Учитывая, что аг = Сг-11 (2 < г < т), перепишем уравнения (18) в виде

с11рт2

с121Рт2

1/2 — £

3=3 т

1/3— £

3=3

Разделив второе уравнение на первое, получим

1/3 — Е С3-1,1рт3

7 = 3=3_

"11 т '

1/2 — Е с3-1,1Рт3

3=3

т

Данное значение с'11 обеспечивает совместность уравнений (18). Значения Рт2 и Рт1 последовательно определим из второго и первого уравнений (10). Оставшиеся коэффициенты А3 (3 < г < т, 1 < ] < г — 1) вычислим последовательным решением систем линейных алгебраических уравнений Вквк = с* (2 < к < т — 1), где вк = (вк+1;1, •••, вк+1,к)Т, ск = (ск1, • • • , скк)Т

В качестве примера приведем коэффициенты 10-стадийной схемы (4) второго порядка точности с согласованными областями устойчивости промежуточных численных схем:

Р10,1

Р10,3 Р10,5 Р10,7 Р10,9 в2,1 вз,2

в4,2 в5,1

в5,3 вб,1

вб,3

вб,5

в7,2

в7,4

в7,б

в8,2

в8,4

в8,б

в9,1

в9,3 =

в9,5 =

в9,7 =

в10,1

в10,3

в10,5

в10,7

^10,9

а3 =

«5 =

а7 =

«9 =

= -1.8196042548247 = +0.62780912355711 = +0.52697647868521 = +0.25850897771127 = +0.10582824603966 = -7.5165266543482 = -0.40442926460761 ■ 10 = -0.86161426365635 ■ 10 = -0.15541344297494 = +0.24288745824190 = -0.37816232408515 = +0.41790691370223 = +0.55277464597430 ■ 10 = +0.41579093026965 ■ 10 = +0.24947672376381

-4 -4

-1

-3

= +0.53002565495431 ■ 10-1 = +0.83091373687116 ■ 10-3 = +0.36693156810609 = +0.11566112969376 = -1.2883182174482 = +0.77379662163441 = +0.30819901982081 = +0.11081342034211 = -1.5783549552468 = +0.75090999599718 = +0.41555458184504 = +0.17963250597238 = +0.50730034923075 ■ 10-1 +2.46572640299832 ■ 10-2 +1.48519331295003 ■ 10-1 +3.51419025544931 ■ 10-1 +6.35203175855605■ 10

-1

Р10,2 :

Р10,4 :

Р10,б :

Р10,8 :

Р10,10

в3,1

в4,1

в4,3

в5,2

в5,4

вб,2

вб,4

в7,1

в7,3 в7,5 в8,1

в8,3

= +0.26171232237173 ■ 10 = +0.70107890176425 = +0.37388421552143 = +0.17246666567217 = +0.50434522649909■10

-2

-1

в8,5 в8,7 в9,2 в9,4 = в9,б = в9,8 = в10,2 в10,4 в10,б в10,8 а2 = а4 = аб = «8 = «10 =

= +0.24697706956444 = -0.17271889464125 = +0.94543917346749 = -0.43222611482215 = +0.61088538639525 = +0.12174369114793 = +0.14473316234684 = -0.66049210371349 = +0.58451948281918 = +0.12449656624973 = -0.97345739728368 = +0.71164946366367 = +0.20973020417453 = +0.51842795293118

-1 -1

10-1 10-4 10

10

-1

-3

-2

= +0.13506048429757 ■ 10 = +0.48747322823252 = +0.19072487421537 = +0.51163624215609 ■ 10-1 = +0.19537733055761 ■ 10-2 = +0.60172655326385 = +0.27750005508315 = +0.10781983872087 -7.51652665434820 +7.71858664562584 ■ 10-2 +2.39876960252498 ■ 10-1 +4.83188677384359 ■ 10-1 = +8.07472383864321 ■ 10-1

В качестве многочленов устойчивости (19) использовались полиномы с максимальным интервалом устойчивости [4], при этом длина интервала устойчивости основной схемы равна 81.112.

3. Контроль точности и устойчивости

С применением (14) видим, что локальная ошибка $п,т т-стадийной схемы (4) имеет вид

7т = ^1/6 - Х3 63^ ^/7 + 0(^4) = (1/6 - 43^3/'7 + 0(^4), (21)

где 43 = ст3 — коэффициент при г3 многочлена устойчивости (20). Тогда согласно [1] для контроля точности вычислений можно использовать оценку ошибки £п,т вида

£п,т =(1/6 - О^2//п + 0(^3).

(22)

т

Величину ега , т можно оценить многими способами. Разлагая кг (1 < г < т) в ряды Тейлора, нетрудно видеть, что

кг — к,- = (аг — «3 )к2/ / + О(к3), г =

где а5 (2 < 5 < т) определены формулой (5). Отсюда следует, что £га,т можно вычислить по приближенной формуле

£п,т ~ [(1/6 — ст3)/(«г — «3)](кг — к^-), 1 < г,^ < т, г = ^ (23)

Для жестких задач характерным является быстрое изменение решения на небольших промежутках. Для того чтобы в этом случае избежать неудовлетворительной точности, естественно в (23) использовать стадии, вычисленные в крайних точках отрезка [¿га,^га+1]. Стадия к1 вычисляется в точке ¿п, и поэтому в (23) положим ] = 1 (а1 = 0), г > ]. Значение а2 = Сц выбирается из условия совместности уравнений (18), и поэтому может принимать, вообще говоря, произвольное значение. Для методов с согласованными областями устойчивости имеет место неравенство а3 < а4 < • • • < 1. Отсюда следует, что ни одна стадия в точке ¿п+1 не вычисляется. С другой стороны, применение в (23) стадий с достаточно большими номерами будет приводить к существенному увеличению вычислительных затрат в случае повторных вычислений решения вследствие невыполнения требуемой точности расчетов. С целью повышения эффективности алгоритма интегрирования поступим следующим образом.

В качестве предварительной оценки точности вычислений будем применять величину £гт, определенную по формуле

/ _ 1/6 — Ст3 /7 / \

£п,т = -(к2 — к1^

а2

Так как к1 зависит от размера шага линейно, то нарушение неравенства

11^га,т 11< £ (24)

приведет всего лишь к одному дополнительному вычислению правой части задачи (3). Здесь £ — требуемая точность расчетов, || ■ || — некоторая норма в . Учитывая, что

к/(уга+1) — к1 = к2 / / + О(к3),

окончательное решение по точности примем на основе неравенства

1к77,т11< £, (25)

где

£",т = (1/6 — 43)/[к/(Уп+1) — к1]^

Выбор шага по точности будем осуществлять следующим образом. Пусть приближенное решение в точке ¿п вычислено с шагом кп. Учитывая, что £гт = О(к2) и £»гт = О (к2), определим число по формуле

^К.тН = е

Если < 1, т. е. нарушено неравенство (24) и требуемая точность не выполняется, то положим Л,га равным Л,га и повторим вычисление к2. В противном случае определим число из соотношения

,т11 =

Если д2 < 1, то нарушено неравенство (25). Тогда Л,га полагается равным д2Л.га и стадии к (2 < г < т) вычисляются заново. При д2 > 1 вычисляется приближенное решение, а новый шаг Л.га+1 задается формулой

Л,га+1 = шт(дь ф)^. (26)

Выбор числа стадий будем осуществлять с использованием неравенства для контроля устойчивости численной схемы [1]. Для построения данного неравенства применим метод (4) для решения (3) при / (у) = Ау + Ь, где А и Ь соответственно матрица и вектор с постоянными коэффициентами размерности N. Нетрудно видеть, что в этом случае выполнены соотношения = &2 = 7 + «2^77, кз = / + ОзЛ2// + а2вз2Лз/2/™, где / = Ауп + Ь, / = А. Выберем д» и д'' (1 < г < 3) из условия выполнения соотношений

з з

£ д»к» = , £ = /(27)

г=1 г=1

Нетрудно видеть, что равенства (27) будут выполнены, если

= -, ь2 = -0^, ь3 1

1 «2вз2 ' 2 «2вз2 ' 3 «2^32 '

6'/ = --, 62' = -, &з' = 0.

«2 «2

Теперь оценку модуля максимального собственного числа Л^ах матрицы Якоби д/(уп)/ду задачи (3) можно вычислить степенным методом по приближенной формуле

йЛГх = 1—1—г шах

|«2^з2| 1<?<М

"га

[а2^з + аз^2 - («2 + «з)^

[к2 - ^

Тогда неравенство для контроля устойчивости т-стадийного метода (4) имеет вид

^ЛГх < 7т, (28)

где 7т — длина интервала устойчивости т-стадийной схемы.

4. Алгоритм интегрирования с переменным числом стадий

Пусть имеется набор т-стадийных методов, длины интервалов устойчивости которых обозначим через 7т (3 < т < М, М > 3). На каждом шаге выбор числа стадий будем осуществлять следующим образом. Пусть приближенное решение в точке ¿п вычислено с шагом Л,га по т-стадийной схеме. Определим Л4+1 по формуле (26).

Если т < М и кП+1АГх > 7т, то т положим равным (т +1), число д3 вычислим из условия устойчивости по формуле

?3кП+1 А^аХ = 7т,

а новый шаг кга+1 перевычислим следующим образом:

кга+1 = ш1п(дь ^2, (29)

Если т > 3 и кга+1АГХ < 7т-1, то т положим равным (т — 1). В остальных случаях число стадий не изменяется. При т = М новый шаг задается формулой (29), при т =3 — соотношением (26). В случае расчетов с максимальным числом стадий применение (29) позволяет удержать шаг в пределах области устойчивости.

5. Результаты расчетов

В состав алгоритма интегрирования переменного шага и с переменным числом стадий включены методы второго порядка с числом стадий т (3 < т < 14). Для определенности выбраны численные схемы с максимальным интервалом устойчивости. Коэффициенты многочленов устойчивости (см. ниже) получены с помощью алгоритма из [1]:

т = 3 73 = —6.260700000

С33 = = 0.6250000000

т= 4 74 = —0.1204670000 ■ 102

С43 = = 0.7808448345 ■ 10-1 С44 = 0.3608453922 ■ 10-2

т= 5 75 = —0.1945690000 ■ 102

С53 = = 0.8460849927 ■ 10-1 С54 = 0.5527124819 ■ 10-2 С55 = 0.1221964350 ■ 10-3

т= 6 76 = —0.2850430000 ■ 102

С63 = = 0.8799401907 ■ 10-1 С64 = 0.6616916777 ■ 10-2 С65 = 0.2217607053 ■ 10-3

С66 = = 0.2731155893 ■ 10-5

т= 7 77 = —0.3919240000 ■ 102

С73 = = 0.8998502098 ■ 10-1 С74 = 0.7287754889 ■ 10-2 С75 = 0.2929815057 ■ 10-3

С76 = = 0.5723750735 ■ 10-5 С77 = 0.4336798850 ■ 10-7

т= 8 78 = —0.5152260000 ■ 102

С83 = = 0.9125773964 ■ 10-1 С84 = 0.7728176610 ■ 10-2 С85 = 0.3436678727 ■ 10-3

С86 = = 0.8297336203 ■ 10-5 С87 = 0.1029826713 ■ 10-6 С88 = 0.5148094796 ■ 10-9

т= 9 79 = —0.6549570000 ■ 102

С93 = = 0.9212164140 ■ 10-1 С94 = 0.8032277127 ■ 10-2 С95 = 0.3804328437 ■ 10-3

С96 = = 0.1037334639 ■ 10-4 С97 = 0.1627525710 ■ 10-6 С98 = 0.1365234306 ■ 10-8

С99 = = 0.4743117465■ 10-11

т= 10 710 = —0.8111200000 ■ 102

С10,3 = 0.9273532641 ■ 10- 1 С10,4 = 0.8250827248 ■ 10-2 С10,5 = 0.4077305837 ■ 10-3

с10,6 = 0.1202172903 ■ 10- 4 С10,7 = 0.2165863427 ■ 10-6 с10,8 = 0.2337894537 ■ 10-8

с10,9 = 0.1388784147 ■ 10- 10 с10,10 = 0.3490928048 ■ 10-13

т= 11 711 = —0.9837160000 ■ 102

С11,3 = 0.9318712290 ■ 10- 1 С11,4 = 0.8413065880 ■ 10-2 С11,5 = 0.4284624834 ■ 10-3

С11,6 = 0.1333201614 ■ 10- 4 С11,7 = 0.2630173525 ■ 10-6 С11,8 = 0.3304691889 ■ 10-8

С11,9 = 0.2562757224 ■ 10- 10 с11,10 = 0.1118194634 ■ 10-12 С11,11 = 0.2099977764 ■ 10-15

т= 12 712 = —0.1172747000 ■ 103

с12,3 = 0.9352947408 ■ 10- 1 С12,4 = 0.8536760476 ■ 10-2 С12,5 = 0.4445343203 ■ 10-3

С12,6 = 0.1438143468 ■ 10- 4 С12,7 = 0.3023697970 ■ 10-6 С12,8 = 0.4204580146 ■ 10-8

с12,9 = 0.3838519723 ■ 10- 10 с12,10 = 0.2212616523 ■ 10-12 с12,11 = 0.7302820006 ■ 10-15

С12,12 = 0.1051890200 ■ 10 -17

m = 13 Yi3 = —0.1378213000 ■ 103

c13,3 = 0.9379514494 ■ 10-i Ci3,4 = = 0.8633199686 ■ 10-2 Ci3,5 = 0.4572230222 ■ 10-3

c13,6 = 0.1523025589 ■ 10-4 Ci3,7 = 0.3355378847 ■ 10-6 c13,8 = 0.5014834871 ■ 10-8

ci3,9 = 0.5112962591 ■ 10-i0 Ci3,i0 = 0.3502954352 ■ 10-i2 с13,11 = 0.1542745108 ■ 10-14

ci3,12 = 0.3946094014 ■10-i7 с13,13 = 0.4455721670 ■ 10-20

m= 14 Yi4 = -0.1600115000 ■ 103

Ci4,3 = 0.9400547623 ■ 10-i c14,4 = 0.8709829298 ■ 10-2 ci4,5 = 0.4674036548 ■ 10-3

Ci4,6 = 0.1592403480 ■ 10-4 Ci4,7 = = 0.3635021510 ■ 10-6 Ci4,8 = 0.5732072002 ■ 10-8

c14,9 = 0.6328016128 ■ 10-i0 c14,10 = 0.4879793010 ■ 10-i2 c14,11 = 0.2575379337 ■ 10-14

Ci4,i2 = 0.8865299187 ■ 10-i7 Ci4,13 = 0.1793358233 ■ 10-i9 Ci4,14 = 0.1617028584 ■ 10-22

Соответствующие области устойчивости "почти" многосвязные. В [1] описан алгоритм расчета коэффициентов многочленов устойчивости, при которых область устойчивости имеет заданные, естественно "разумные" форму и размер. Поэтому при необходимости полученные коэффициенты можно перевычислить.

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

Построенный алгоритм интегрирования применялся для решения уравнения Ван-дер-Поля

у! = У2, У2 = 100(1 - у2)у2 - yi, yi(0) = 2, y2(0) = 0, h0 = 2 ■ 10-2 (0 < t < 1000). (30)

Для (30) характерно многократное чередование переходных участков с участками установления, что позволяет не только оценить эффективность построенного алгоритма интегрирования, но и проверить работу алгоритма выбора числа стадий.

Для решения задачи (30) с точностью 10-2 с использованием построенного алгоритма потребовалось 78 734 вычислений if правой части. Если области устойчивости основной и промежуточных численных схем не согласованы, то решение тем же самым алгоритмом задачи (30) приводит к if = 87 311, причем дополнительные затраты возникают за счет повторных вычислений решения (возвратов) вследствие невыполнения требуемой точности расчетов, а требуемая точность нарушается вследствие неустойчивости промежуточных численных формул.

Численный эксперимент проводился с целью продемонстрировать повышение эффективности за счет согласования областей устойчивости. Дальнейшее более эффективное использование данных методов мы видим в составе более мощного алгоритма переменных порядка, шага и с переменным числом стадий [1]. Однако данный алгоритм эффективнее многих известных явных методов и может быть использован самостоятельно, в частности, при решении задачи (30) методом Мерсона [4] (if = 363 195), методом Фельберга [5] (if = 366 384), с помощью алгоритма интегрирования с переменными порядком и шагом на основе методов Адамса (if = 396 927).

Список литературы

[1] Новиков Е. А. Явные методы для жестких систем. Новосибирск: Наука, 1997.

[2] Golushkü M. I., NoviKüv E. A. Explicit fourth-order methods for stiff system // Russ.

J. Numer. Anal. Math. Modelling. 1999. Vol. 4, No. 1. P. 71-85.

[3] Новиков В. А., Новиков Е. А. О построении явных методов типа Рунге — Кутта с расширенными областями устойчивости. Красноярск: ВЦ СО РАН, 1998.

[4] MersüN R. H. An operational methods for integration processes // Proc. Symp. on Data Proc. Weapons Research Establishment. Salisbury, Australia, 1957.

[5] FEHLBERG E. Klassische Runge - Kutta - Formeln fünfter und siebenter Ordnung mit Schrittweitenkontrolle // Computing. 1969. No. 4. P. 93-106.

Поступила в редакцию 20 февраля 2001 г.