CC BY
55
d
Тогда дифференциальный оператор L = — — Q(t) : W),(
Н) с Ь2 ^ Ь2 обратим и норма обратного оператора имеет соответствующую (каждому случаю теоремы) оценку:
х{Я)
1) IIL II <
2) ||L
-i i
<
xo(Q)x(Q) — sup II ddQ (t)|
te r dt
x(Q)2
Xo(Q)x(Q)2 — sup IIQ(t)dQ(t) — dQQ(t)I
dt
dt
Работа выполнена при финансовой поддержке РФФИ (проект 10-01-00276). Библиографический список
1. Баскаков А. Г. Дихотомия спектра несамосопряженных операторов // Сиб. мат. журн. 1991. Т. 32, № 3. С. 24-30.
2. Баскаков А. Г., Юргелас В. В. Круговая дихотомия
спектра одного класса несамосопряженных операторов // Изв. вузов. 1994. № 10. С. 12-18. 3. Годунов С. К. Современные аспекты линейной алгебры. М. : Научная книга, 1997. 407 с.
УДК 519.622
АЛГОРИТМ ИНТЕГРИРОВАНИЯ ЖЕСТКИХ ЗАДАЧ С ПОМОЩЬЮ ЯВНЫХ И НЕЯВНЫХ МЕТОДОВ
Е. А. Новиков
Институт вычислительного моделирования СО РАН, Красноярск
E-mail: [email protected]
Построены L-устойчивый (3,2)-метод третьего порядка и явная трехстадийная схема типа Рунге-Кутты первого порядка точности. Создан алгоритм интегрирования переменного порядка и шага, в котором выбор эффективной численной схемы осуществляется на каждом шаге с применением неравенства для контроля устойчивости. Приведены результаты расчетов, подтверждающие эффективность построенного алгоритма.
Ключевые слова: жесткие задачи, явный и неявный методы, контроль точности и устойчивости, переменный порядок.
Algorithm of Integrating Stiff Problems Using the Explicit and Implicit Methods
E. A. Novikov
An L-stable (3,2)-method order 3 and an explicit three-stage Runge-Kutta scheme order 1 are constructed. An integration algorithm of variable order and step is constructed that is based on of the two schemes The most effective numerical scheme is chosen for each step by means of stability control. The results are given that confirm the effectiveness of the algorithm.
Key words: stiff problems, explicit and implicit methods, stability and accuracy control, variable order.
ВВЕДЕНИЕ
Во многих важных приложениях возникает необходимость решения жестких задач. Основные тенденции при построении численных методов связаны с расширением их возможностей для решения систем все более высокой размерности. В современных методах решения жестких задач при вычислении стадий применяется Ьи-разложение некоторой матрицы, размерность которой совпадает с размерностью исходной системы дифференциальных уравнений. Разложение осуществляется с выбором главного элемента по строке или столбцу, а иногда и по всей матрице. В случае достаточно большой размерности быстродействие алгоритма интегрирования, фактически, полностью определяется временем декомпозиции этой матрицы. Для повышения эффективности расчетов в ряде алгоритмов используется замораживание матрицы Якоби, т. е. применение одной матрицы на нескольких шагах интегрирования [1]. Наиболее успешно этот подход применяется в многошаговых методах [2]. Не вызывает эта проблема особых трудностей и при построении алгоритмов интегрирования на основе других численных схем, если в них стадии вычисляются с участием матрицы Якоби в некотором итерационном процессе. В этом случае она не влияет на порядок точности численной схемы, а определяет только сходимость итераций. Поэтому необходимость в ее пересчете возникает при существенном замедлении скорости сходимости итерационного процесса.
В алгоритмах интегрирования на основе известных безытерационных методов, к которым относятся методы типа Розенброка [3] и их различные модификации [1], проблема замораживания является более трудной. Следует отметить, что с точки зрения реализации безытерационные методы существенно проще алгоритмов на основе численных формул, в которых стадии вычисляются с применением итераций. Однако в методах вида [3] матрица Якоби влияет на порядок точности численной схемы, и поэтому возникают трудности с ее замораживанием. В [4] доказано, что максимальный порядок точности методов типа Розенброка равен двум, если в алгоритме интегрирования одна матрица Яко-би применяется на нескольких шагах интегрирования. Там же построен алгоритм с замораживанием матрицы Якоби на основе L-устойчивой численной формулы второго порядка, и приведены результаты расчетов, подтверждающие его высокую эффективность при невысокой точности расчетов. Если для таких методов вопрос об использовании одной матрицы на нескольких шагах интегрирования оставить не решенным, то нужно ограничиться либо задачами небольшой размерности, либо невысокой точностью расчетов.
Некоторым аналогом замораживания матрицы Якоби является применение в расчетах алгоритмов интегрирования на основе явных и L-устойчивых методов с автоматическим выбором численной схемы. В этом случае эффективность алгоритма может быть повышена за счет расчета переходного участка, соответствующего максимальному собственному числу матрицы Якоби, при помощи явного метода [5]. В качестве критерия выбора эффективной численной формулы естественно применять неравенство для контроля устойчивости [6-7]. Следует отметить, что применение таких комбинированных алгоритмов полностью не снимает проблемы замораживания матрицы Якоби, потому что явным методом можно просчитать, вообще говоря, только переходный режим, соответствующий максимальному собственному числу матрицы Якоби.
Здесь на основе явного трехстадийного метода типа Рунге-Кутты первого порядка и L-устойчивого (3,2)-метода третьего порядка точности построен алгоритм переменной структуры, в котором допускается замораживание как численной, так и аналитической матрицы Якоби. Приведены результаты расчетов, подтверждающие эффективность алгоритма интегрирования.
1. КЛАСС (М,К)-МЕТОДОВ
Рассмотрим задачу Коши для жесткой системы дифференциальных уравнений:
у' = f (У), V(to) = Уо, to < t < tk, (1)
где y и f — вещественные N-мерные вектор-функции, t — независимая переменная, которая меняется на заданном интервале [t0,tk]. Зададимся целыми числами m и k, k < m, и введем в рассмотрение множества
Mm = {1, 2,..., m}, Mk = {mi e Mm|1 = mi < m2 < ■ ■ ■ < mk < m}, Mm-k = Mm \ Mk, Ji = {mj-1 e Mm|j > 1, mj e Mk, mj < i}, 2 < i < m.
Для решения задачи (1) будем применять (m, к)-схемы вида [8]
Vn+1 = Vn + Pi ki + P2k2 + . . . + Pmkm, Dn = E - ükAn, i-1
Dnki = kf (yn + ^^ eij kj
aij kj, i e Mk, (2)
j=i jeJi
Dnki = ki-1 + ^ ] aij kj , i e Mm-k,
jeJi
где матрица An представима следующим образом:
An = fn + kBn + O(k2), (3)
E — единичная матрица, fn = df (yn)/dy — матрица Якоби системы (1), k - шаг интегрирования, Bn — некоторая матрица, не зависящая от размера шага интегрирования, a, pi, в-ij ,aij — вещественные константы, определяющие свойства точности и устойчивости (2). В отличие от других известных
одношаговых численных схем, в (m, ^-методах правая часть задачи (1) вычисляется не для всех стадий. Для описания вычислительных затрат на шаг интегрирования в традиционных одношаговых методах достаточно одной константы m — числа стадий. В данных схемах затраты определяются двумя постоянными — m и k. В (m, к)-методах на каждом шаге один раз вычисляется матрица Якоби и осуществляется декомпозиция матрицы Dn, k раз вычисляется функция f, m раз осуществляется обратный ход в методе Гаусса. Условие (3) позволяет применять (2) с использованием одной матрицы Якоби на нескольких шагах, которая может быть вычислена как аналитически, так и численно [9]. В отличие от схем типа Розенброка, в рамках (m, к)-методов просто решаются проблемы замораживания матрицы Якоби и ее численной аппроксимации.
2. L-УСТОЙЧИВЫЙ (3,2)-МЕТОД ТРЕТЬЕГО ПОРЯДКА
Для решения задачи (1) рассмотрим численную схему:
Уп+i = Уп + Piki + Р2k2 + Рзкз, Dn = E - ahAn, Dn ki = hf (yn), Dn k2 = ki, Dn кз = hf (Уп+i) + «32 k2, (4)
где матрица An удовлетворяет условию (3). Численную формулу
Уп+i = Уп + 03iki + вз2 k2 (5)
называют внутренней или промежуточной схемой метода (4). Для построения метода третьего порядка разложим стадии k¿, 1 < i < 3, в ряды Тейлора в окрестности точки yn до членов с h3 включительно и подставим в первую формулу (4). Получим
Уп+i = Уп + [pi + Р2 + (1 + «32 )рз ]hfn + [api + 2ap2 + (a + 03i + вз2 + 3a«32)p3]h2 fn fn+ + [a2pi + 3a2p2 + (a2 + 2a^3i + 3a^32 + 6а2аз2 )p3]h3f2fn+
+1(e3i + ^32)2P3h3 fn'fn + [api + 2ap2 + (a + 3ааз2 )рз ]h3 Bn fn + O(h4),
где элементарные дифференциалы вычислены на приближенном решении уп . Разложение точного решения y(tn+i) в ряд Тейлора в окрестности точки tn до членов с h3 включительно имеет вид
y(tn+i) = y(tn) + hf + 2 h2 f 'f + 1 h3 (f '2f + f'' f2) + O(h4),
где элементарные дифференциалы вычислены на точном решении y(tn). Сравнивая эти ряды при yn = y(tn), получим условия третьего порядка схемы (4), т.е.
pi + p2 + (1 + «32)p3 = 1, api + 2ap2 + (a + 03i + вз2 + 3a«32)p3 = 2,
a2pi + 3a2p2 + (a2 + 2a^3i + 3a^32 + 6a2«32 )p3 = 1, (6)
b
(03i + вз2)2 p3 = 3, api + 2ap2 + (a + 3a«32 )p3 = 0.
Последнее соотношение (6) возникает за счет применения «испорченной» матрицы Якоби. Исследуя совместность системы (6), получим коэффициенты
5 3 1 3 3 д 2 д
pi = 4 + 4"32, p2 = -1 - ^"32, p3 = 4, в32 = 3 - P3i,
2 3 3 2 3 1 . .
-a + -a + -a аз2 -7a^3i = 77, (7)
2 4 4 b
где e3i и a32 — свободные параметры.
3. ИССЛЕДОВАНИЕ УСТОЙЧИВОСТИ
Исследуем устойчивость схемы (4) на линейной тестовой задаче
y' = Ay, y(0) = У0, t > 0, (8)
где Л — произвольное комплексное число, Ив (Л) < 0. Смысл Л — некоторое собственное число матрицы Якоби задачи (1). Применяя (4) для решения (8), получим уп+1_ = Яз(х)уп, где х = ЛН, а функция устойчивости Яз (х) имеет вид
Яз(х) = {(-а3 + а2р1 + а2рз - авз1 Рз)х3 + [За2 - 2ар1 - ар2+ +(вз1 + вз2 - 2а)рз]х2 + [-За + р1 + р2 + (1 + «32)рз]х + 1}/(1 - аx)3.
Из вида Яз(х) следует, что для ¿-устойчивости основной схемы (4) необходимо выполнение следующего равенства:
а2 - а(р1 + рз) + АлРз =0. (9)
Запишем условие ¿-устойчивости промежуточной схемы (5). Применяя ее для решения (8), получим уп+1 = Я2(х)уп, где функция устойчивости Я2(х) имеет вид
1 + (вз1 + вз2 - 2а)х + а(а - вз1)х2
Я2(х) =-(г-ах^-•
Из вида Я2(х) следует, что схема (5) будет ¿-устойчивой, если вз1 = а. Учитывая соотношение (9), из (7) имеем, что коэффициенты метода (4), обладающего ¿-устойчивостью основной и промежуточной численных схем, записываются следующим образом: р1 = в31 = а, р2 = -2а + 3/2, рз = 3/4, в32 = 2/3 - а, а32 = 4а/3 - 5/3, где а определяется из условия ¿-устойчивости:
6аз - 18а2 + 9а - 1 = 0. (10)
Уравнение (10) имеет три корня: а1 = 2,40514957850286, а2 = 0,158983899988677 и аз = 0,435866521508459. Согласно [10], схема (4) будет А-устойчива, если 1/3 < а < 1,0685790. Поэтому выбираем корень а = аз = 0,435866521508459.
4. КОНТРОЛЬ ТОЧНОСТИ ВЫЧИСЛЕНИЙ
Для контроля точности вычислений используем идею вложенных методов. Для этого введем дополнительную стадию Бпк4 = кз. Теперь на основе вычисленных ранее стадий к1 , к2 , кз и к4 построим численную схему второго порядка
Уп+1,2 = Уп + Ь1 к1 + Ь2к2 + Ьз кз + 64 к4 (11)
и будем контролировать точность расчетов с применением приближений, вычисленных по формулам (4) и (11), т.е.
||Уп+1 - Уп+1,2 II < е, (12)
где е — требуемая точность расчетов, || ■ || — некоторая норма в Ян. Для нахождения коэффициентов Ь1, Ь2, Ьз и Ь4 запишем разложение стадий к1, к2, кз и к4 в ряды Тейлора и подставим в (11). Сравнивая ряды для точного и приближенного решений, запишем условия второго порядка точности схемы (11):
Ь1 + Ь2 + (1 + аз2)(Ьз + Ь4 ) = 1,
аЬ1 + 2аЬ2 + (а + вз1 + вз2 + 3ааз2 )Ьз + (2а + вз1 + вз2 + 4ааз2 )Ь4 = 2. (13)
Напомним, что схема (4) построена с учетом возможности замораживания матрицы Якоби и ее численной аппроксимации. Ошибки, которые вносятся использованием «испорченной» матрицы Якоби, устраняются за счет последнего соотношения (6). Аналогичное требование к формуле (11) приводит к дополнительному уравнению на коэффициенты Ь^ т.е.
аЬ1 + 2аЬ2 + (1 + 3аз2 )аЬз + (2 + 4аз2 )аЬ4 = 0. (14)
Решая совместно (13) и (14), имеем
1 /4 2\ 4 2 3
Ь1 = 2а - ^ -( 3а - 3) Ьз, Ь2 = 2 - 3а + (3а - 3^3, Ь4 = 4 - Ьз,
где Ьз — свободный параметр. Положив Ьз = 0, получим коэффициенты Ь1 = 2а - 1/2, Ь2 = 2 - 3а, Ь4 = 4/3.
Применение неравенства (12) не приводит к дополнительным вычислениям правой части или обращениям матрицы Якоби. Введение стадии к4, используемой только для контроля точности, увеличивает вычислительные затраты на один обратный ход метода Гаусса на шаг интегрирования. Отметим одну важную особенность построенного неравенства (12). Схема (4) ¿-устойчивая. Следовательно, д3(х) ^ 0 при х ^ —го. Так как для точного решения у(£п+1) = ехр(х)у(£п) задачи (8) выполняется аналогичное свойство, то естественным будет требование стремления к нулю оценки ошибки. Однако для разности [уп+1 — Уп+1,2] это свойство не выполняется, потому что схема (11) ¿-устойчивой не является. С целью исправления асимптотического поведения оценки ошибки вместо (12) будем контролировать неравенство
Р1-" (Уп+1 — Уп+1,2) || < е, 1 < ¿п < 2. (15)
В этом случае поведение оценки ошибки при ¿п = 2 будет согласовано с поведением точного решения задачи (8) при х ^ —го. Подчеркнем, что в смысле главного члена оценки ошибок, применяемые в неравенствах (12) и (15), совпадают при любом ¿п. Использование неравенства (15) вместо (12), фактически, не приводит к увеличению вычислительных затрат. Это связано с тем, что (15) при ¿п = 2 проверяется только в том случае, если оно нарушено при ¿п = 1. Такая ситуация встречается достаточно редко, в основном при быстром росте величины шага. Однако это позволяет выбирать шаг более правильно и тем самым уменьшить количество неоправданных повторных вычислений решения, что повышает надежность расчетов и быстродействие алгоритма интегрирования.
Оценку максимального собственного числа г>п,0 = НЛп,тах матрицы Якоби системы (1), необходимую для перехода на явную схему, оценим по формуле г>п,0 = Н||д//ду||.
5. ЯВНЫЙ МЕТОД ПЕРВОГО ПОРЯДКА
Для численного решения задачи (1) рассмотрим схему вида
Уп+1 = Уп + Г1&1 + Г2 ^2 + Гзкз, = Н/ (Уп), ^2 = Н/ (уп + 021^),
кз = Н/(Уп + вз1 + вз2^2), (16)
где Н — шаг интегрирования, , к2 и кз — стадии метода, г1, г2, гз, 021, вз1, и 0з2 — числовые коэффициенты, определяющие свойства точности и устойчивости (16). Получим соотношения на коэффициенты метода (16) первого порядка точности. Для этого разложим стадии к1, к2 и кз в ряды Тейлора до членов с Н2 и подставим в первую формулу (16). В результате получим
Уп+1 = Уп + (Г1 + Г 2 + Гз)Н/п + [в21Г2 + (вз1 + $л)Гз ]Н2 /П /п + 0(Нз).
Здесь элементарные дифференциалы вычислены на приближенном решении уп. Точное решение у(^п+1) в окрестности точки £п имеет вид
у(*п+1) = у(*п) + Л/ + 1Н2/' / + 0(Нз),
где элементарные дифференциалы вычислены на точном решении у(£п). Построим схему первого порядка с максимальным интервалом устойчивости. Для этого применим (16) для решения тестового уравнения (8). Получим уп+1 = ^(х)уп, где функция устойчивости ф(х) имеет вид
д(х) = 1 + (Г1 + Г 2 + Гз )х + [в21Г2 + (вз1 + вз2)Гз ]х2 + 0210з2^з хз, х = НЛ.
Требование первого порядка точности приводит к соотношению г1 + г2 + гз = 1, которое далее будем считать выполненным. Выберем г2 и гз так, чтобы метод (16) имел максимальный интервал устойчивости. Для этого рассмотрим многочлен Чебышева Тз(г) = 4гз — 3г на промежутке [—1,1]. Проведем замену переменных, полагая г = 1 — 2х/7. Получим
() — 18 48 2 — 32 з
Т з (х) — 1 х \ х ^ х , ^ ^2 ^з
при этом отрезок [7, 0] отображается на отрезок [-1,1]. Среди всех многочленов вида
Р3 (х) = 1 + х + с2х2 + сзх3
для Тз(х) неравенство |Т3(х)| < 1 выполняется на максимальном интервале [7,0], 7 = -18 [7]. Потребуем совпадения коэффициентов Я(х) и Тз(х). Это приводит к соотношениям
44 Г1 + Г2 + Гз = 1, в21 Г2 + (вз1 + вз2 )Гз = ^, 021 ^32^3 = 729 ^
В результате имеем коэффициенты
Г3 = 749(в21в32)-1, Г2 = [27 - (в31 + в32)г3]в2-1' Г1 = 1 - Г2 - Г3
метода первого порядка с максимальным интервалом устойчивости. Положим в21 = 0.5, в31 = -1 и в32 = 2. Тогда на каждом шаге приращения к1, к2 и кз вычисляются в точках Ьп, £п + 0.5Н и £п + Н соответственно. Так как стадии к1 и кз вычисляются в крайних точках интервала интегрирования, то повышается надежность расчетов. Более того, при данных ¡3^ легко построить метод (16) третьего порядка точности [11]. В результате имеем коэффициенты метода (16), т.е.
1 517 208 4
в21 = 2, вз1 = -1> вз2 = 2, Г1 = 729' Г2 = 729' Г3 = 729' (17)
а его локальная ошибка <5пдимеет вид 5п,1 = 19Н2/'//54 + 0(Н3). Метод первого порядка предполагается использовать на участке установления, где ошибки невелики. Поэтому для контроля точности численной формулы (16), (17) можно использовать оценку локальной ошибки. Учитывая, что имеет место
к2 - к1 = 1Н2/п/п + 0(Н3)
и вид локальной ошибки, неравенство для контроля точности записывается в виде
19
^ ||к2 - М < е,
где || ■ || — некоторая норма в Ян, е — точность расчетов.
6. КОНТРОЛЬ УСТОЙЧИВОСТИ
Неравенство для контроля устойчивости численной формулы (16) построим предложенным в [7] способом. Запишем стадии к1, к2 и кз применительно к задаче у' = Ау, где А есть матрица с постоянными коэффициентами. В результате получим
к1 = Хуп, к2 = (X + 0.5Х2 )уп, кз = (X + X2 + X3 )у,п,
где X = НА. Легко видеть, что имеют место соотношения
кз - 2к2 + к1 = X3уп, 2(к2 - к1)= X2уп.
Тогда согласно [7] оценку максимального собственного числа г>пд = НЛптах матрицы Якоби системы (1) можно вычислить по формуле
^ 1 = 0.5 шах |к3 -.2к2 + к11. (18)
Щ1 1<г<^ |к2 - к11
Интервал устойчивости схемы (16), (17) равен 18 [7]. Поэтому для ее контроля устойчивости используется неравенство г>пд < 18.
В случае применения полинома Чебышева в качестве многочлена устойчивости область устойчивости «почти» многосвязная [11]. Появление мнимых частей собственных чисел матрицы Якоби за счет ошибок округлений может приводить к ее сокращению. Алгоритмом из [7] построим многочлен устойчивости
Я(х) = 1 + х + с2 х2 + сз х3,
для которого в экстремальных точках имеет место |ф(х)| = 0.9. Получим с2 = 0.15625736489384 и сз = 0.0061526400319, а коэффициенты г метода первого порядка определяются формулами г3 = сз, г2 = 2(с2 - сз), г1 = 1 - 2с2 + сз .В итоге имеем
г1 = 0.69363791024424, Г2 = 0.30020944972383, Гз = 0.0061526400319238.
(19)
Область устойчивости метода (16), (19) приведена на рисунке, где изображены линии уровня |ф(х)| = 5 при 5 равном 0.5, 0.7 и 1. Длина интервала устойчивости (16), (19) равна примерно 17 (см. рисунок, 7 = -16.93). В расчетах применялись коэффициенты (19), а неравенство для контроля устойчивости имеет вид г>пд < 17, где г>пд определяется по формуле (18).
Оценка (18) является грубой, потому что вовсе не обязательно максимальное собственное число сильно отделено от остальных, в степенном методе применяется мало итераций и дополнительные искажения вносит нелинейность задачи (1). Поэтому контроль устойчивости используется как ограничитель на размер шага интегрирования. Прогнозируемый шаг Нп+1 задается следующим образом. Новый шаг Нас по точности определяется по формуле Нас = q1Нn, где Нп есть последний успешный шаг интегрирования, а д1, учитывая соотношение к2 - к1 = О(Нп), задается уравнением
—2
19 27
?2 ||к2 - к11 = е.
Область устойчивости метода (16), (19)
Шаг Наг по устойчивости задается формулой Наг = q2Нn, где д2, учитывая равенство г>пд = 0(Нп), определяется из уравнения q2г>пд = 17. В результате Нп+1 вычисляется по формуле
Нп+1 = шах[Нп, шт(Нас , Н^)].
(20)
Отметим, что формула (20) применяется для прогноза величины шага интегрирования Нп+1 после успешного вычисления решения с предыдущем шагом Нп, и поэтому, фактически, не приводит к увеличению вычислительных затрат. Если шаг по устойчивости меньше последнего успешного, то он уменьшен не будет, потому что причиной этого может быть грубость оценки максимального собственного числа. Однако шаг не будет и увеличен, потому что не исключена возможность неустойчивости численной схемы. Если шаг по устойчивости должен быть уменьшен, то в качестве следующего шага будет применяться последний успешный шаг Нп. В результате для выбора шага и предлагается формула (20). Данная формула позволяет стабилизировать поведение шага на участке установления решения, где определяющую роль играет устойчивость. Собственно говоря, именно наличие данного участка существенно ограничивает возможности применения явных методов для решения жестких задач.
7. АЛГОРИТМ ПЕРЕМЕННОЙ СТРУКТУРЫ
Первый шаг интегрирования всегда выполняется по явной схеме, потому что она не использует матрицы Якоби. Нарушение неравенства г>пд < 17 вызывает переход на схему (4). Передача управления явному методу происходит в случае выполнения неравенства г>п,0 < 17, где при расчетах ¿-устойчивым (3,2)-методом оценка максимального собственного числа г>п,0 вычисляется через норму матрицы Якоби. Численную формулу (4) без потери порядка точности можно применять с замораживанием матрицы Бп. Отметим, что при замораживании матрицы Якоби величина шага интегрирования остается постоянной. Попытка замораживания матрицы Бп осуществляется после каждого успешного шага. Размораживание матрицы происходит в следующих случаях: 1) нарушение точности расчетов, 2) если число шагов с замороженной матрицей достигло заданного максимального числа ^Н, 3) если прогнозируемый шаг больше последнего успешного в qН раз. Числами ^Н и qН
можно влиять на перераспределение вычислительных затрат. При ¿дН = 0 и дН = 0 замораживания не происходит, при увеличении ¿дН и дН число вычислений функции / возрастает, а число декомпозиций матрицы убывает.
Норма ||£|| в неравенствах для контроля точности вычисляется по формуле
|уп| + г'
где г — номер компоненты, г — положительный параметр. Если по г-й компоненте решения выполняется неравенство |уп| < г, то контролируется абсолютная ошибка ге, в противном случае — относительная ошибка е. Ниже построенный алгоритм переменной структуры будем называть МК3КК1. Алгоритм переменного шага на основе только (3,2)-метода назовем МК3. Отметим, что МК3 является частью МК3КК1 и подключается по признаку. В алгоритмах МК3КК1 и МК3 матрица Якоби перевычислялась в том случае, если прогнозируемый шаг в полтора раза превышает успешный (дН = 1.5). Максимальное число шагов с замороженной матрицей ограничивалось числом десять (гдН = 10).
8. РЕЗУЛЬТАТЫ РАСЧЕТОВ
В качестве тестовой задачи рассмотрено уравнение Ван-дер-Поля для имитации осциллирующих физических процессов [1]. В ней переходные режимы чередуются с участками установления, что позволяет исследовать возможности построенного алгоритма. Уравнение Ван-дер-Поля изучается в виде системы двух уравнений у1 = у2, у2 = [(1 — у2)у2 — у1]/д с начальными условиями у1(0) = 2 и у2(0) =0, Ь е [0,11]. Изменением параметра д варьируется жесткость модели. Сравнение эффективности построенного алгоритма МК3КК1 проводилось с МК3 и методом Гира в реализации А. Хиндмарша ЭЬЗОЭЕ из коллекции ОЭЕРАСК [2]. В таблице приведены результаты расчетов при различных значениях д. Вычислительные затраты приведены в форме г/(¿¿'), где через г/ и обозначены соответственно число вычислений правой части и декомпозиций матрицы Якоби на интервале интегрирования. Расчеты проводились таким образом, чтобы в точном и приближенном решениях совпадали три значащие цифры, где под точным понимается решение при точности вычислений е = 10-9. В алгоритмах МК3КК1, МК3 и ЭЬЗОЭЕ матрица Якоби вычислялась численно.
Результаты расчетов
10-1 10-2 10-3 10-4 10-5 10-6
МК3ЯК1 1297(0) 2964(0) 3243(338) 4362(430) 5047(532) 5809(631)
МК3 1056(84) 1462(241) 3148(373) 4343(487) 5037(536) 5844(685)
ОЬБООЕ 916(92) 1974(234) 2996(368) 4236(535) 4974(562) 5823(723)
При небольшой жесткости задачи в построенном алгоритме МК3КК1 работает только явный метод. Начиная с д = 10^ в зависимости от устойчивости в процессе вычислений комбинируются ¿-устойчивая и явная схемы. Изучение пошаговых вычислений показывает, что на переходных участках расчеты осуществляются по явной формуле, а на участках установления — по ¿-устойчивой схеме. Из сравнения результатов расчетов алгоритмами МК3КК1 и МК3 следует при примерно одинаковом числе вычислений правой части сокращение числа декомпозиций за счет расчетов на части интервала по явному методу. По числу декомпозиций матрицы Якоби построенный алгоритм МК3КК1 при всех значениях д эффективнее алгоритма ЭЬЗОЭЕ.
ЗАКЛЮЧЕНИЕ
В МК3КК1 с помощью признака можно задавать различные режимы расчета: 1) явным методом первого порядка точности с контролем или без контроля устойчивости; 2) ¿-устойчивым методом с замораживанием или без замораживания как аналитической, так и численной матрицы Якоби; 3) с автоматическим выбором численной схемы. Все это позволяет применять данный алгоритм для решения как жестких, так и нежестких задач. При расчетах с автоматическим выбором численной схемы вопрос о том, является ли задача жесткой или нет, перекладывается на алгоритм интегрирования. Наиболее эффективное применение данного алгоритма предполагается при точности вычислений 10-4 и грубее.
Использование неравенства для контроля устойчивости, фактически, не приводит к увеличению вычислительных затрат, потому что оценка максимального собственного числа матрицы Якоби системы (1) осуществляется через ранее вычисленные стадии и не приводит к росту числа вычислений функции /. Такая оценка получается грубой. Однако применение контроля устойчивости в качестве ограничителя на рост шага позволяет избежать негативных последствий грубости оценки. Более того, в некоторых случаях это приводит к нестандартно высокому повышению эффективности алгоритма [11].
Работа выполнена при финансовой поддержке РФФИ (проекты 11-01-00106 и 11-01-00224).
Библиографический список
1. Hairer E., Wanner G. Solving ordinary differential equations II. Stiff and differential-algebraic problems. Berlin : Springer-Verlag, 2004. 614 p.
2. Byrne G. D., Hindmarsh A. C. ODE solvers: a review of current and coming attractions // J. of Comp. Phys. 1987. № 70. P. 1-62.
3. Rosenbrock H. H. Some general implicit processes for the numerical solution of differential equations // Computer. 1963. № 5. P. 329-330.
4. Новиков В. А., Новиков Е. А., Юматова Л. А. Замораживание матрицы Якоби в методе типа Розенброка второго порядка точности // Журн. вычисл. мат. и мат. физ. 1987. Т. 27, № 3. С. 385-390.
5. Новиков Е. А. Построение алгоритма интегрирования жестких систем дифференциальных уравнений на неоднородных схемах // Докл. АН СССР. 1984. Т. 278, № 2. С. 272-275.
6. Новиков В. А., Новиков Е. А. Повышение эффективности алгоритмов интегрирования обыкновенных дифференциальных уравнений за счет контроля устойчиво-
УДК 519.71
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ДИНАМИЧЕСКОГО ХАОСА
В. А. Подчукаев
Саратовская государственная юридическая академия E-mail: [email protected]
Поставлена и решена задача аналитического конструирования по заданной математической модели динамической системы в пространстве состояний сопровождающей её математической модели в фазовом пространстве. Показано, что изображающая точка всякого решения динамической системы общего вида в пространстве состояний принадлежит гиперсфере со смещённым центром в фазовом пространстве (или эквивалентной ей центральной гиперсфере переменного радиуса). Сконструировано аналитическое представление центра смещения, объясняющее происхождение динамического хаоса бесконечными разрывами второго рода в координатах центра смещения. Показано, что эти разрывы порождены переходом через ноль соответствующих компонент вектора состояний.
Ключевые слова: динамическая система, обыкновенные однородные дифференциальные уравнения, гиперсфера, смещённый центр, динамический хаос.
сти // Журн. вычисл. мат. и мат. физ. 1985. Т. 25, № 7. С. 1023-1030.
7. Новиков Е. А. Явные методы для жестких систем. Новосибирск : Наука, 1997. 197 с.
8. Новиков Е. А., Шитов Ю. А., Шокин Ю. И. Од-ношаговые безытерационные методы решения жестких систем // Докл. АН СССР. 1988. Т. 301, № 6. С. 13101314.
9. Новиков A. E., Новиков E. A. Численное решение жестких задач с небольшой точностью // Математическое моделирование. 2010. Т. 22, № 1. С. 46-56.
10. Демидов Г. В., Юматова Л. А. Исследование точности неявных одношаговых методов. Препринт № 11. ВЦ СО АН СССР. Новосибирск, 1976. 22 с.
11. Новикова Е. А. Алгоритм переменного порядка и шага на основе явного трехстадийного метода типа Рунге-Кутта // Изв. Сарат. ун-та. Нов. сер. 2011. Т. 11. Сер. Математика. Механика. Информатика, вып. 3, ч. 1. С. 46-53.
Mathematical Model of Dynamic Chaos V. A. Podchukaev
The problem of analytical designing on the set mathematical model of dynamic system in space of states of mathematical model accompanying it in phase space is put and solved. It is shown, that the representing point of any decision of dynamic system of a general view in space of states conditions belongs to hypersphere with the displaced centre in phase space (or to central hypersphere of variable radius equivalent to it). Analytical representation of the centre of the displacement, an explaining origin of dynamic chaos by infinite ruptures of the second sort in co-ordinates of the centre of displacement is designed. It is shown, that these ruptures are generated by transition through a zero corresponding a component of a vector of states.
Key words: dynamic system, ordinary homogeneous differential equations, hypersphere, displaced centre, dynamic chaos.
