Построение и анализ гибридного РК-метода тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

Введение.

0.1. Стратегия метода.

0.2. Порядок аппроксимации и упрощающие условия

Бутчера.

0.3. Устойчивость методов Рунге-Кутты.

0.4. Обзор методов Рунге-Кутты.

0.4.1. Суперточные РК-методы.

0.4.2. Нильпотентные РК-методы.

0.5. Обзор диссертации по

главам.

Глава 1.

Решение нелинейных РК-систем методом предобуславливания.

1.1. Реализация РК-метода.

1.2. Алгоритм Бутчера для однократно неявного РК-метода.

1.3. Идея предобуславливания.

1.4. Исследование предобуславливателя.

1.5. Численные эксперименты.

Глава 2.

Построение обобщенного преобразованного РК-метода.

2.1. Построение преобразованного РК-метода.

2.2. Обобщение преобразованного РК-метода.

2.3. Алгоритм построения обобщенного преобразованного РК-метода.

2.4. Численные эксперименты.

Глава 3.

Гибридный РК-метод - новая стратегия решения систем ОДУ.

3.1. Описание стратегии.

3.2. Описание драйвера.

3.3. Обоснование цели создания гибридного РК-метода.

3.4. Численные эксперименты.

3.5. Выводы из численных экспериментов.

0.1. Стратегия метода Рунге-Кутты.

Пусть дано обыкновенное дифференциальное уравнение = У (?)), te[0,Tl y(0) = y0,

1) и требуется найти его приближенное решение уп в дискретных точках tn, tn g [0,Т]. В методе Рунге-Кутты ([21],[29],[32]), если решение уп в точке *„ известно, то для нахождения приближенного решения v„, в точке (п„л = tn - т последовательно находятся приближенные решения 77, в промежуточных узлах £ с

Таким образом, если формально проинтегрировать уравнение (1) на интервале [/„ ,£.],/ = \(\)m +1, и квадратурную формулу построить на узлах £ с весами Д., / = 1(1 )т, то решение уп+1 будет описываться системой нелинейных алгебраических уравнений п ' 'h+i y&) = y„+jf(t,y(t))dt,

Ъ =УИ,

3=1

Уп+1 — Vm+\~> с обозначениями ft = tn+A,T, i = 0(1)/я +1, A0=0, <ЛМ, Лт+1 - 1, где т- число разбиений интервала [t„,t„+J. Будем говорить, что задание параметров Ду, Л,,, т определяет т -стадийную РК-схему. Введем матричные и векторные обозначения: fin fin . А.

J2 fill fill • • • filn fi ml fi ml •■• fimr, (fim+1,1 J ■•• > fim+\,m )'

Г] = (ъ,.,?1т)Т , / = (/„.,/„ )7 e=(l,.,l)r.

Тогда система (2) запишется в виде: ri = yne + rBf, Уп+1 = Л + dm+J

На практике интерес представляет векторная форма уравнения (1), в которой y(t) и f(t,y) являются векторами порядка N :

Я0 = (/°(0,.,.у(Л,)(0)г,

Поэтому в (2) под величинами уп, rji, /( также следует понимать векторы:

3) fj=(fj°\.,fr)r

Таким образом, векторы ц = ) и / = /(7) = (//,.JmT) имеют порядок mN. Тогда в (2) к -е компоненты векторов 7,, /, / = l(l)w, одновременно являются ((/'-!)# + £)-ми компонентами векторов т/, /.

Записать систему уравнений (2) непосредственно относительно векторов rj и / можно, воспользовавшись прямым произведением матриц (или произведением Кронекера). Кронекеровым произведением [22] двух квадратных матриц U = (иг]) и V = {vtj) порядка к и I соответственно, обозначаемым (U ®V), называется матрица

F) = unV . порядка kl. В этой записи матрица {U®V) представлена как блочная, в которой блок, расположенный в / -й строке и j -м столбце, есть матрица и/.

Тогда система (3) запишется в виде: ri = e<8>yn+T(B®EN)f(Tj), (4)

5) где En - единичная матрица порядка N.

В дальнейшем будем использовать следующее свойство произведений Кронекера [22]:

А <8> В) ■ (С <8> £>) = АВ ® CD .

0.2. Порядок аппроксимации и упрощающие условия Бутчера.

Система алгебраических уравнений (2) является аппроксимацией задачи (1). Обозначив через y(tn+x) точное решение уравнения (1) в точке tn+l, а через j)n+1 решение системы (2) после замены уп на y(tn), определим величину называемую локальной погрешностью аппроксимации, которая характеризует точность метода решения.

Метод Рунге-Кутты называют имеющим порядок аппроксимации М, если М - наибольшее натуральное число, для которого выполняется

0(гЛ/+1) при г —О для любого достаточное число раз дифференцируемого решения уравнения (1). Величина М определяется методом и не зависит от каких-либо свойств дифференциального уравнения. Порядок конкретного метода (причем вместе с системой уравнений для определения его параметров Ду и А,) можно найти с помощью разложения в ряд Тейлора локальной погрешности аппроксимации. Однако такой способ приводит к слишком громоздким выражениям. Для большинства методов порядок аппроксимации может быть найден с помощью упрощающих условий, введенных Бутчером ([9],[13]).

Говорят, что т -стадийный метод Рунге-Кутты удовлетворяет упрощающему условию

А(к), если 1п+1=0(тм), m 1

В(к), если =-, \<г<к,

1 Г m 1

С(к), если У ДД/ =.Я/ , 1 < г < к, 1 < i < m . j=i r

Условие A(k) означает, что метод Рунге-Кутты имеет порядок аппроксимации не менее к.

Условие В(к) означает, что порядок точности квадратурной формулы, аппроксимирующей интеграл n+1 y(tn+x)=yn + \f(t,y(t))dt, tn равен к.

Условие C(k) означает, что порядок точности каждой из квадратурных формул, аппроксимирующих интегралы

4, y(Z) = yn+jf(t,y(№, f„ больше или равен к (следовательно, величины % аппроксимируют решение в промежуточных узлах с порядком к). Введя обозначения f\ Л, .

W =

1 /L

1 1

Л.

Л. т-1 то-1

- матрица Вандермонда, а ^ - ее к -й столбец, перепишем условия В(к) и С(к) в матричном виде: В(к):

С(к): b^W.T=-eTt \<г<к г

BW.r=-AW.r, 1 <г<к. г

Упрощающие условия позволяют вывести различные методы Рунге-Кутты [16]. Пусть заданы т различных узлов Л,,.,Лт. Тогда матрицы Л и W невырождены, и условие В(ш) образует систему т линейных уравнений для весов /?т+и,.,/?т+иг. Следовательно, для них существует единственное решение. Потребовав, чтобы метод удовлетворял условию С(т), получаем систему линейных уравнений, полностью определяющую матрицу В.

Таким образом, если заданы различные узлы Л,,.,Лт, то упрощающие условия В(т) и С(т) однозначно определяют некоторую схему Рунге-Кутты.

0.3. Устойчивость методов Рунге-Кутты.

В практически важных задачах имеют дело с функциями достаточно узкого класса, среди которых особое место занимают диссипативные функции.

Для определения этого понятия [16] рассмотрим кусочно-непрерывную функцию v(t): [О, Т] -> R такую, что при t е [0, Т] справедливо неравенство г,) - / 01, z2 \zx-z2)< v( t)\zx -z2f \/zx ,z2 еД, где Dt - область определения функции /(/,•): Dt ~> RN, (*,*) - скалярное произведение в RN, ||*|| - евклидова норма.

Пусть у it), y(t) - два решения уравнения (1) с начальными условиями у0

Если v(0 < 0, то есть

-^2)<о, то [16] j/(f2) - y(t2 )|| < \\y(t, )->'(/,)[, о < t, < t2 < Т.

Следовательно, с увеличением t два любых решения «притягиваются» друг к другу. Это свойство называется контрактивностью. Дифференциальные уравнения, обладающие этим свойством, называются диссипативными. Соответственно векторная функция / (/,•) также называется диссипативной. Простейшее из диссипативных уравнений. = 0, t> 0, у(О) = у0, at используют в качестве модели для исследования устойчивости численных методов решения нелинейных систем. Его решение на интервале [tn,tn^] имеет вид y(tn+l) = eTy(tn). Численная модель, с другой стороны, дает yn+l = R(r)yn, где R: С С - полином или рациональная функция с вещественными коэффициентами, причем [21]

R(r) = \ + Tbm^E-rBy'e. Таким образом, возникает проблема аппроксимации экспоненты некоторой функцией R(t) , которая называется функцией устойчивости. Ее выражение через коэффициенты характеристического многочлена От(л), т i=i матрицы Бутчера В имеет вид [21]: т ( 1 J а Л

Vl- Чк 1 т Л к=1

Выражения для коэффициентов q,, j = 0(l)m, характеристического многочлена матрицы В будут даны в разделе 0.4.2.).

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

Метод Рунге-Кутты называется абсолютно устойчивым для некоторого zgC, если для этого z выполняется неравенство

ВД| < 1, из которого следует, что |у„+,|<|.у„|. Множество всех точек геС, для которых справедливо неравенство |/?(z)| < 1, называется областью абсолютной устойчивости метода. Абсолютная устойчивость - это минимальное свойство, которому должна удовлетворять РК-схема, чтобы подавлять вычислительные возмущения.

Если область абсолютной устойчивости содержит левую полуплоскость Re(z) < О, то метод называется А-устойчивым. А-устойчивый метод с условием R(-оо) = 0 называется L-устойчивым. А-устойчивость была впервые введена Далквистом [14], а L-устойчивость - Эйле [36].

0.4. Обзор методов Рунге-Кутты.

Метод, предложенный К. Рунге, усовершенствованный В. Куттой и названный методом Рунге-Кутты (РК-методом), развивается с 1895 года. Сначала это был только явный метод, ориентированный на гладкие скалярные функции и на очень маленький шаг интегрирования.

Во второй половине XX века во многих областях (в химической кинетике, в управлении горением в соплах ракет и самолетов, в каталитических реакторах и др.) возникла потребность в решении систем нелинейных обыкновенных дифференциальных уравнений. Это привело к появлению неявных методов, в которых решение получалось итерационно.

Одношаговый РК-метод подвергся модификации в двух направлениях.

1) Розенброк [26] модифицировал явный РК-метод, введя в него матрицу Якоби. Позднее было осознано, что этот метод эквивалентен одной итерации DIRK-метода (см. O.4.2.), по этой причине он налагает очень большие ограничения на шаг. Этот метод широко используется и сейчас [25], а его эффективность обеспечивается нетривиальным «замораживанием» матрицы Якоби и правой части.

2) С 1963 года благодаря Бутчеру ([9]-[13]), Александеру [2], Бурриджу ([4]-[8]), Эйле [36], Чипману ([30]-[31]) и другим начинает развиваться теория неявного РК-метода, опирающаяся на замену дифференциального уравнения на интервале интегральным и значительно расширяющая область применения методов Рунге-Кутты. Бутчер предложил неявную

РК-схему, основанную на квадратуре Гаусса. Впоследствии на основе квадратур Радо и Лобатто был предложен ряд других РК-схем. Все вместе они представляют суперточный РК-метод, отличительными особенностями которого являются А-устойчивость, а в некоторых схемах даже L-устойчивость, и высокая аппроксимация, поскольку А,,, / = 1(1)/??, являются корнями соответствующих ортогональных многочленов.

Суперточные РК-схемы обладают серьезным недостатком: они требуют решения системы линейных алгебраических уравнений порядка mN. При численном решении это приводит к большому объему арифметических действий и к большим погрешностям (что воспринимается даже как потеря устойчивости).

Поэтому интенсивно развивались и так называемые нильпотентные методы ([20],[21]). К ним относятся три группы методов: явный, диагонально неявный и однократно неявный РК-методы. Все они основаны на квадратуре Ньютона-Котеса. Нильпотентные РК-методы уступают суперточным в порядке аппроксимации и устойчивости, однако они значительно проще в реализации, так как для них существуют алгоритмы, обеспечивающие эффективное решение линейной системы.

Компактное представление метода Рунге-Кутты дает так называемая таблица Бутчера [16]:

К А, - А*, ~

Ля Ал 1 ••• Р mm

Рт±\,1 •■• Рт+\,т

В зависимости от выбора параметров Д., Я и т получаются различные РК-схемы.

Остановимся подробнее на свойствах двух основных классов методов Рунге-Кутты: суперточных (методах повышенного порядка точности) и нильпотентных.

Ае В

0.4.1. Суперточные РК-методы.

К суперточным РК-методам относятся методы Гаусса, Радо и Лобатто, основанные на квадратурных формулах высокого порядка. Узлами Д , i = 1(1 )т, суперточных РК-методов являются корни характеристического многочлена матрицы А ([19]-[21]):

V' М т-со.

ZH)J о м

J J м-j т —со- j

Л1 m-j

Здесь и далее используется обозначение

U. С/=

7! йО-О!

Его представление в форме Родриго имеет вид: m! J

•М-т

-{X

М-т+со

А-1Г-),

М! где М - порядок аппроксимации РК-метода, а со - кратность корня Л = 0 многочлена пт(Л). Чтобы матрица Вандермонда была невырожденной, необходимо наложить ограничение со- 0 или со = 1.

Рассмотрим основные свойства каждого из суперточных РК-методов. 1) Порядок аппроксимации т -стадийного метода Гаусса М = 2т [9]. Все его узлы лежат внутри интервала [o,l], следовательно со = 0. В этом случае характеристический многочлен матрицы Л, nmG (Л) = \лт (Л - \)т 1,

2 т)\с1Лт1 J представляет собой многочлен Лежандра, определенный на отрезке [од] ([1],[15]). Метод Гаусса удовлетворяет условию С(т). Функция устойчивости метода Гаусса имеет вид: т к Z

R(z)

-к f2 f0k\K т I к-0

ЫУ к\

2т - к тп ■

Метод Гаусса является А-устойчивым, но не является L-устойчивым, так как |Л(-оо)| = 1.

2) Порядок аппроксимации т -стадийных методов Радо 1 и Радо 2 М - 2т -1. Эти методы характеризуются тем, что первый или последний узел совпадает с соответствующим концом интервала интегрирования .

В методе Радо 2 Л1 = 0, следовательно со = 1, а Я, являются корнями многочлена хнт (Я) = т[ [г (Я -1)- ],

2т -1)! dX 1 J который может быть выражен через характеристический многочлен РК-метода Гаусса [1]: xfift^J'W + xJ'W

В методе Радо 1 Хт = 1, со = 0, а Я,, / = l(l)m -1, являются корнями многочлена

А) = {Х "1Г1= 0 " Я) + (1" А) • Узлы методов Радо 1 и Радо 2 связаны соотношениями

V =1- i = \(\)m.

Метод Радо 1 удовлетворяет условию С(т), а метод Радо 2 условию C(m-l) ([31],[36]). Функция устойчивости методов Радо имеет вид: г* (2т-\-к^ Z

Ti к\

R(r) т у к\ т и хорошо аппроксимирует экспоненту при отрицательных г, поскольку д(-оо) = о. Методы Радо L-устойчивы.

3) Методы Лобатто имеют порядок аппроксимации М = 2т-2. Первый и последний узлы совпадают с концами интервала интегрирования: Я, =0,

Ят = 1, то есть со = 1, а остальные узлы Я,., / = 2(l)m -1, являются корнями многочлена я! dm~2 лт (Л) =

2т-2)!с/Ат-2 который выражается через многочлен РК-метода Гаусса [1]:

Узлы методов Лобатто симметрично расположены относительно центра интервала [од]. Для метода Лобатто 1 [36] выполняется условие С(ш), а для метода Лобатто 2 [31] - условие C(m-l). Функция устойчивости метода Лобатто 1 имеет вид

1 Тк Г2/Я-2-АЛ

R(T) =

V J y(zй к к\

2т - 2 - к

V т~х J ' а метода Лобатто 2

R(T) =

2т-2-к^ т I к(2т-2-к^ ' т -- 2 0 ^ откуда видно, что метод Лобатто 1 А-устойчивый, а метод Лобатто 2 -L-устойчивый.

Функция устойчивости каждого из суперточных РК-методов представляет собой аппроксимацию Паде экспоненты [23], которая является наилучшей аппроксимацией экспоненты не только в вещественной, но и в мнимой области.

Достоинством суперточных РК-методов является высокий порядок аппроксимации и А-устойчивость. Однако их матрицы коэффициентов В содержат очень мало нулевых элементов, а это означает, что все (или почти все) уравнения в нелинейной системе (4) РК-метода связаны. Следовательно, при реализации этих методов надо решать систему из mN нелинейных уравнений, из-за чего эти методы значительно проигрывают в сравнении с другими.

0.4.2. Нильпотентные РК-методы.

Большая группа методов Рунге-Кутты основана на свойстве нильпотентности матриц ([19]-[21]). Будем рассматривать случай, когда матрица (В - //£) - нильпотентная с индексом нильпотентности т, то есть jj- т-кратное собственное число матрицы В. Матрица В с таким свойством называется полной, нильпотентный РК-метод с полной матрицей В так же называется полным.

Теория нильпотентных РК-методов связана с многочленами JIareppa At (/О [1], определяемыми формулой j=° Jи удовлетворяющими рекуррентным соотношениям

Lx= 0, L0 = 1, к +1 )4+1 (Л) = (Л - 2к -1 )Lk (Л) - кЬкЛ (Л). Порядок аппроксимации полного нильпотентного РК-метода М - т, если Ljr>(~)^ 0, и М = т +1, если выполнено еще и дополнительное м условие Lm+lm(—) = 0. В нильпотентном РК-методе порядок аппроксимации м не может быть больше т +1 [12]. Здесь т 1

I (1) (Л) = У (-1)"' J р (./-/)!

- первая производная многочлена JTareppa порядка т .

Функция устойчивости нильпотентного РК-метода [6] имеет вид m i

1-/«■)"

А-устойчивость зависит от выбора так как |Л(г)| не всегда меньше единицы.

Если в нильиотентном РК-методе матрица В - нижняя треугольная, то метод называется диагонально неявным (DIRK-метод) ([2],[24]). Его недостатком является то, что в промежуточных узлах порядок аппроксимации не может быть больше 1, то есть для него выполнено упрощающее условие Бутчера С(1). Если в DIRK-методе положить // = О, то в этом случае Ду = 0 при j > i для всех / и вектор / может быть вычислен явным образом по значениям Поэтому такие формулы называются явными. Для явного РК-метода дробно-рациональная функция устойчивости вырождается в многочлен: т -j у=1 fоткуда следует плохая устойчивость явных РК-схем. Более того, многочлен может аппроксимировать экспоненту лишь в малой окрестности, поэтому шаг т не может быть большим.

Необходимо отметить, что и дробно-рациональная аппроксимация нильпотентных РК-схем, в отличие от суперточных, не является наилучшей, поскольку коэффициенты характеристического многочлена подчинены другому условию, а именно, условию нильпотентности.

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

В общем случае матрица В с единственным собственным значением /л не является нижней треугольной. Такие методы называются однократно неявными (SIRK-методы) ([4]-[8]). В отличие от DIRK-методов, они удовлетворяют упрощающему условию Бутчера С(т-1), что обеспечивает в промежуточных узлах порядок аппроксимации т -1.

Условие С(ш-1), 1

BW.k = —AW,k, к к = \{\)т-\, определяет только т-1 первых столбцов матрицы BW. Неизвестный нам вектор BW.m, в предположении |Л| ф 0, может быть представлен [21] в виде mBW.m=DAW.m, где D- некоторая диагональная матрица. С помощью этого представления условие С(т) выражается в матричном виде т-\ к=1 т-1 1111 +~AW-»> +~DAw.m = k=\ к m m m

6) AWh + A(D-E)W.meJh, где h = diag(\,\/2,.,l/m), em = (0,.Д1)Г .

Вместо того, чтобы определять элементы самой матрицы В, часто удобнее рассматривать преобразованную матрицу, такую, например, как матрица В, полученная в результате преобразования подобия [28]

В = W^BW = W~xAWh + W~lMP - E)W.meJh . (7)

Это преобразование возможно в силу невырожденности матрицы W. Заметим, что

W~lAW = Fr, где о 1

F =

О 1 Pi Рх

- матрица Фробениуса [22] с характеристическим многочленом т

Pf (X) = X"' -^Р/""' ■ i= 1

Введя обозначение вектора z - W~l A(D - E)W,m, перепишем выражение (7) в виде

В = FT h + zejh = (FT + zej )h = FzTh. (8)

Через Fz обозначим матрицу, полученную из F заменой pi на р, + zmM :

О 1

F. О 1

Pn,+Zi ••• Рг+г«,-1 Л+2, Характеристический многочлен матрицы B = F2Th, а следовательно и В вычисляется рекуррентно и равен т - /)!

Р, + •

1 >»'

С другой стороны, характеристический многочлен нильпотентной матрицы В индекса нильпотентности т может быть представлен в виде: т т 1 м

Приравнивая в этих двух представлениях одного и того же многочлена коэффициенты при одинаковых степенях Л, найдем m-i)\ т\ откуда 14/1 Г! г I

Р, + Zm-M = ("I) 7--Т7 т И , т-г)\ i = 1(1 )т

9)

10)

Из (10) видно, что фиксируя тот или иной вектор 2, мы определяем коэффициенты характеристического многочлена матрицы В, а следовательно, задаем узлы Л{, i = 1(1 )m. В наиболее простом случае z = 0 выражения (9) и (10) принимают вид: ^ZlAp 5 7 = 1(1)W; (11) я!

А-НУ'-^^У, i=KO». 02) да -/)!

Соотношения (6) при этом сводятся к виду BW = AWh, т.е. к условию С(т), которое означает, что величины ц. аппроксимируют решение дифференциального уравнения (1) в промежуточных узлах , /' = \(\)т, с порядком т.

РК-метод, удовлетворяющий условию С(т), называется коллокационным [5]. Для него определяющими являются соотношения (11), связывающие характеристические многочлены матриц В и А. Иными словами, в коллокационном РК-методе матрица В полностью определяется упрощающим условием С(ш) по заданному набору узлов, Я являющихся корнями характеристического многочлена кт{Х) = т\jumLm(—). М

Характеристический многочлен матрицы Л однократно неявного РК-метода в самом общем случае на основании (9) записывается в виде

2 т У- ,=1

По своим свойствам однократно неявные методы уступают суперточным, так как они условно устойчивы и имеют порядок аппроксимации т или т +1. Однако они обладают важным преимуществом, поскольку для них существует эффективный и экономичный алгоритм решения линейной системы, впервые предложенный Бутчером [11], в котором число арифметических операций пропорционально 0{Nl), в отличие от 0(ni3Nl) в суперточных методах Рунге-Кутты.

Представляет интерес создание нового РК-метода, который объединил бы в себе достоинства как суперточных, так и однократно неявных методов Рунге-Кутты: А-устойчивость, высокий порядок аппроксимации и экономичное решение нелинейной системы.

0.5. Обзор диссертации по главам.

В первой главе рассматривается решение системы нелинейных уравнений в методе Рунге-Кутты, которая методом Ньютона-Рафсона сводится к линейной порядка mN. Основная идея первой главы заключается в том, что линеаризованная система трактуется как схема решения нелинейной системы методом предобусдавливания. В этом случае невязка строится по матрице В суперточного РК-метода, что гарантирует высокий порядок аппроксимации и хорошую устойчивость. В состав предобуславливателя входит матрица В однократно неявного РК-метода, для которого существует алгоритм Бутчера, снижающий число арифметических операций при решении линеаризованной системы с 0(тьЫъ) до 0(jV3). Матрица однократно неявного РК-метода, в свою очередь, строится на узлах суперточного. Таким образом, введением идеи предобуславливания удается расширить возможности применения однократно неявных и суперточных РК-схем, объединив их преимущества в одном методе, названном методом Ньютона-Рафсона с предобуславливанием. Доказывается возможность построения матрицы SIRK-метода на узлах суперточных РК-схем Гаусса и Радо 1, а также исследуется близость предобусловленной матрицы к исходной.

Во второй главе предлагается расширение класса однократно неявных РК-методов. В преобразованном РК-методе Бурриджа при специальном выборе константы с получаются схемы «с забеганием вперед», в которых последний промежуточный узел выходит за пределы интервала интегрирования. Решение в этом узле вычисляется дважды, как в последнем промежуточном узле на предыдущем шаге интегрирования и как в первом промежуточном узле - на текущем шаге интегрирования. Сравнение полученных результатов дает возможность проанализировать поведение решения и на основе этого анализа гибко скорректировать стратегию вычислений. Порядок аппроксимации преобразованного РК-метода Бурриджа равен т . Повысить порядок аппроксимации до М - т +1 позволяет введение дополнительного параметра v, в результате чего получается обобщенный преобразованный РК-метод. Таким образом, при специальном выборе параметров /j, v и с получается новый РК-метод порядка т +1 со свойством «забегания вперед». В этой главе исследуется выбор параметра v, предлагается алгоритм построения обобщенного преобразованного РК-метода и численно изучается выбор «забегающего» узла.

Третья глава посвящена построению гибридного РК-метода, который является комбинацией метода Ньютона-Рафсона с предобуславливанием и схемы «с забеганием вперед», описанных в предыдущих главах. Гибридный РК-метод представляет собой новую стратегию решения систем ОДУ в методе Рунге-Кутты. Новизна заключается в следующем: предлагается основные вычисления вести методом Ньютона-Рафсона с предобуславливанием, а на контрольных шагах применять схему «с забеганием вперед». Метод Ньютона-Рафсона с предобуславливанием обеспечит высокий порядок аппроксимации и А-уСтойчивость, а схема «с забеганием вперед» даст возможность проанализировать поведение решения и при необходимости изменить шаг и скорректировать параметры счета.

3.5. Выводы из численных экспериментов.

Анализ результатов экспериментов позволяет сделать следующие выводы:

- если начальный шаг интегрирования выбран удачно, то гибридный РК-метод проигрывает в экономичности по сравнению с методом Ньютона-Рафсона с предобуславливанием на 5-10%;

- если начальный шаг меньше оптимального для данной задачи, то гибридный РК-метод оказывается значительно экономичнее;

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

Таким образом, построенный гибридный РК-метод может быть рекомендован как эффективный метод решения систем ОДУ, обеспечивающий хорошую устойчивость, высокую точность и экономичность вычислений.

Заключение

В диссертации предлагается новый метод решения системы нелинейных уравнений в методе Рунге-Кутты. В основу метода легла идея трактовки линеаризованной по методу Ньютона-Рафсона системы как схемы решения нелинейной системы методом предобуславливания. В этом случае невязка строится по матрице В суперточного РК-метода, что гарантирует высокий порядок аппроксимации и А- устойчивость. В состав предобуславливателя входит матрица В однократно неявного РК-метода, для которого существует экономичный алгоритм, предложенный Бутчером, уменьшающий число арифметических операций при решении линеаризованной системы с 0{тъЫ3) до 0(ЫЪ). Матрица SIRK-метода, в свою очередь, строится на узлах суперточного РК-метода. Таким образом, введением идеи предобуславливания удается расширить возможности применения однократно неявных и суперточных РК-схем, объединив их преимущества в одном методе.

Предлагается обобщение преобразованного РК-метода Бурриджа, когда введением дополнительного параметра v удается повысить порядок аппроксимации с m до и + 1 и получить схемы «с забеганием вперед», что дает возможность осуществлять контроль решения на шаге интегрирования. Таким образом, путем построения обобщенного РК-метода порядка т +1 со свойством «забегания вперед» удается расширить класс однократно неявных РК-методов.

На основе описанных методов построен гибридный РК-метод, позволяющий гибко менять стратегию вычислений. Основные вычисления ведутся методом Ньютона-Рафсона с предобуславливанием, а на контрольных шагах применяется схема «с забеганием вперед», которая позволяет проанализировать поведение решения и при необходимости изменить шаг и скорректировать параметры счета.

1. Абрамович М., Стиган И. Справочник по специальным функциям. -Москва: Наука, 1979.

2. Александер (Alexander R.) Diagonally implicit Runge-Kutta methods for stiff ODEs. SIAM J. Numer. Anal. 14, 1977, p.1006-1021.

3. Артемьев С.С., Шкурко И.О. Простой, быстрый, надежный алгоритм переменного порядка и шага, основанный на методах типа Розенброка. -Новосибирск: Изд. ВЦ СО АН СССР, 1984.

4. Бурридж (Burrage К.) Stability and efficiency properties of implicit Runge-Kutta methods. Thesis, Dept. of Mathematics, University of Auckland, 1978.

5. Бурридж (Burrage K.) A special family of Runge-Kutta methods for solving stiff differential equations. BIT, 18,1978, p.22-41.

6. Бурридж (Burrage K.) Hight order algebraically stable Runge-Kutta methods. BIT, 18, 1978, p.373-383.

7. Бурридж (Burrage K.) Iteration schemes for singly-implicit Runge-Kutta methods. Computer Science Report № 24, University of Auckland, 1981.

8. Бурридж (Burrage K.) Efficiently iinplementable algebraically stable Runge-Kutta methods. SIAM J. Numer. Anal. 19, 1982, p.245-258.

9. Бутчер (Butcher J.C.) Implicit Runge-Kutta processes. Math. Сотр. 18, 1964, p.50-64.

10. Бутчер (Butcher J.C.) An algebraic theory of integration methods. Math. Сотр. 26, 1972, p.79-106.

11. Бутчер (Butcher J.C.) On the implemantation of implicit Runge-Kutta methods. BIT 16, 1977, p.237-240.

12. Бутчер (Butcher J.C.) A transformed implicit Runge-Kutta method. -Report Series № 111, Dept. of Mathematics, University of Auckland, 1977.

13. Бутчер (Butcher J.C.) The Numerical Analysis of Ordinary Differential Equations. Runge-Kutta and General Linear Methods, Wiley, Chichester-New York-Brisbane-Toronto-Singapore, 1987.

14. Далквист (Dahlquist G.) A special stability problem for linear multistep methods. BIT 3, 1963, p.27-43.

15. Далквист, Бьерк (Dahlkquist G., Bjorck A.) Numerical methods. Prentice Hall, Inc., Englewood Cliffs, New Jersey, 1974.

16. Деккер К., Вервер Я. Устойчивость метода Рунге-Кутты для жестких нелинейных дифференциальных уравнений. Москва.: Мир. 1988.

17. Ильин В.П., Кузнецов Ю.И. Алгебраические основы численного анализа. НСО Наука, 1986.

18. Кузнецов Ю.И. Алгоритмы тройственной структуры. Инф.бюл. «Алгоритмы и программы», 50870001327, ВНИТЦ, 1988, № 9, стр.4.

19. Кузнецов Ю.И. Построение и анализ РК-схем. Новосибирск: Изд. ВЦ СО РАН, 1990.

20. Кузнецов (Kuznetsov Yu.I.) Runge-Kutta method of advanced accuracy: the new point of view// Bull. Nov. Сотр. Centr, Num. Anal. 1994. - № 5 -p.35-53.

21. Кузнецов Ю.И. Алгебраические основы РК метода численного решения ОДУ. - Новосибирск: Изд. ВЦ СО РАН, 1995.

22. Ланкастер П. Теория матриц. 2-е изд. Москва: Наука, 1982.

23. Нерсетт (Norsett S. P.) Multiple Pade approximations to the exponential function. Report Mathematics and Computation, № 4, 1974, Dept. of Mathematics, University of Trondherm.

24. Нерсетт (Norsett S. P.) Semi-explicit Runge-Kutta methods. Report Mathematics and Computation, № 6, 1974, Dept. of Mathematics, University of Trondheim.

25. Новиков E. А. Явные методы для жестких систем. Новосибирск: Наука, Сиб. предприятие РАН, 1997.

26. Розенброк (Rosenbrock Н. Н.) Some general implicit processes for the numerical solution of differential equations. Сотр. J. 5, 1963, p.329-330.

27. Фадеев С.И. Численное исследование модели теплового взрыва. -Новосибирск, ИМ СО АН СССР, 1990.

28. Хайрер, Ваннер (Hairer Е., Wanner G.) Algebraically stable and implementable Runge-Kutta methods of high order. SIAM J. Nuiner. Anal. 18, 1981, p. 1098-1108.

29. Холл Дж., Уайт Дж. Современные численные методы решения обыкновенных дифференциальных уравнений. Москва: Изд. Мир, 1979.

30. Чипман (Chipman F. Н.) Numerical solution of initial value problems using A-stable Runge-Kutta processes, Research Report CSRR 2042, Dept. of AACS, University of Waterloo, 1971.

31. Чипман (Chipman F. H.) A-stable Runge-Kutta processes. BIT 11, 1971, p.384-388.

32. Штеттер X. Анализ методов дискретизации для обыкновенных дифференциальных уравнений. Москва: Мир, 1978.

33. Щур О.Ф. Решение нелинейных РК-систем методом предобусдавливания. Новосибирск: Труды конференции молодых ученых, 1998, с.246-254.

34. Щур (Shchur O.F.) The generalization of the transformed RK-method. Novosibirsk: Bulletin of the Novosibirsk Computing Center, 1999, p.133-139.

35. Щур О.Ф. Новая стратегия решения систем ОДУ в методе Рунге-Кутты. Новосибирск: Труды конференции молодых ученых, 2001, с.289-297.

36. Эйле (Ehle В. L.) On Pade approximations to the exponential function and A-stable methods for the numerical solutions of initial value problems. Thesis, Univ. Waterloo, 1969.