CC BY
30
Библиографический список
1. Бурлуцкая М. Ш., Корнев В. В., Хромов А. П. Система Дирака с недифференцируемым потенциалом и периодическими краевыми условиями // Журн. вычисл. мат. и мат. физ. 2012. Т. 52, № 9. С. 1621-1632.
2. Марченко В. А. Операторы Штурма-Лиувилля и их приложения. Киев : Наук. думка, 1977. 340 с.
3. Djakov P., Mityagin B. Bari-Markus property for Riesz projections of ID periodic Dirac operators // Math. Nachr. 2010. Vol. 283, № 3. P. 443-462.
4. Джаков П. Б., Митягин Б. С. Зоны неустойчивости одномерных периодических операторов Шрёдингера и Дирака // УМН. 2006. Т. 61, № 4(370). С. 77-182. 001: 10.4213/гш2121.
5. Бурлуцкая М. Ш., Хромов А. П. Метод Фурье в смешанной задаче для уравнения с частными производными первого порядка с инволюцией // Журн. вычисл. мат. и мат. физ. 2011. Т. 51, № 12. С. 2233-2246.
Dirac System with Undifferentiable Potential and Antiperiodic Boundary Conditions
V. V. Kornev, A. P. Khromov
Saratov State University, Russia, 410012, Saratov, Astrahanskaya st., 83, KornevVV@info.sgu.ru, KhromovAP@info.sgu.ru
The object of the paper is Dirac system with antiperiodic boundary conditions and complex-valued conditions potential. A new method is suggested for investigating spectral properties of this boundary problem. The method is based on the formulas of the transform operators type. It is rather elementary and simple. Using this method asymptotic behaviour of eigenvalues is specificated and it is proved that eigen and associated functions form Riesz basis with brackets in the space of quadratic summerable two-dimensional vector-functions since eigenvalues may be multiple. The structure of Riesz projection operators is also studied. The results of the paper can be used in spectral problems for equations with partial derivatives of the 1 -st order containing involution.
Key words: Dirac system, spectrum, asymptotics, Riesz basis. References
1. Burlutskaya M. Sh., Kornev V. V., Khromov A. P. Dirac system with non-differentiable potential and periodic boundary conditions. Zh. Vychisl. Mat. Mat. Fix., 2012, vol. 52, no. 9, pp. 1621-1632 (in Russian).
2. Marchenko V. A. Operatory Shturma-Liuvillia i ikh priloxheniia [Sturm-Liouville operators and their applications]. Kiev, Naukova Dumka, 1977, 340 p. (in Russian).
3. Djakov P., Mityagin B. Bari-Markus property for YAK 519.622
Riesz projections of ID periodic Dirac operators. Math. Nachr., 2010, vol. 283, no. 3, pp. 443-462.
4. Djakov P., Mityagin B. S. Instability zones of periodic 1-dimensional Schrodinger and Dirac operators. Russian Math. Surveys, 2006, vol. 61, no. 4, pp. 663-766. DOI: 10.4213/rm2121.
5. Burlutskaya M. Sh., Khromov A. P. Fourier method in an initial-boundary value problem for a first-order partial differential equation with involution. Comput. Math. Math. Phys, 2011, vol. 51, no. 12, pp. 2233-2246.
АЛГОРИТМ ПЕРЕМЕННОГО ПОРЯДКА, ШАГА И ПЕРЕМЕННОЙ КОНФИГУРАЦИИ ДЛЯ РЕШЕНИЯ ЖЕСТКИХ ЗАДАЧ
Е. А. Новиков
Доктор физико-математических наук, главный научный сотрудник отдела вычислительной математики, Институт вычислительного моделирования СО РАН, novikov@icm.krasn.ru
Построено неравенство для контроля устойчивости схемы Ческино второго порядка точности. На основе стадий этого метода построена численная формула первого порядка с расширенным до 32 интервалом устойчивости. На основе ¿-устойчивой (2,1)-схемы и численной формулы Ческино разработан алгоритм переменной структуры, в котором эффективная численная формула выбирается на каждом шаге по критерию устойчивости. Алгоритм предназначен для решения как жестких, так и не жестких задач. Приведены результаты расчетов, подтверждающие эффективность построенного алгоритма.
Ключевые слова: жесткая задача, схема Ческино, (2,1)-метод, контроль точности и устойчивости.
ВВЕДЕНИЕ
Во многих важных приложениях возникает проблема численного решения задачи Коши для системы обыкновенных дифференциальных уравнений. Часто уравнения в частных производных приводятся к обыкновенным дифференциальным уравнениям после дискретизации по пространственным переменным. Полученная таким образом задача, как правило, жесткая и большой размерности. Для решения жестких задач в основном применяются неявные методы, в которых основные затраты приходятся на декомпозицию матрицы Якоби. Для повышения эффективности расчетов в ряде алгоритмов используется замораживание матрицы Якоби, т.е. применение одной матрицы на нескольких шагах интегрирования [1]. Наиболее успешно этот подход применяется в алгоритмах на основе многошаговых численных формул [2]. Эта проблема решается достаточно просто для методов, в которых стадии вычисляются с участием матрицы Якоби в некотором итерационном процессе. В алгоритмах интегрирования на основе известных безытерационных методов, к которым относятся методы типа Розенброка [3] и их различные модификации [1], вопрос о применении одной матрицы на нескольких шагах интегрирования более сложный. В таких алгоритмах матрица Якоби включена в численную схему и ее аппроксимация может приводить к понижению порядка точности. В [4] эта проблема изучается применительно к методам типа Розенброка. Доказано, что максимальный порядок точности данных методов равен двум, если в алгоритме интегрирования одна матрица Якоби применяется на нескольких шагах интегрирования. Это означает, что применение методов типа Розенброка будет эффективным при решении задач небольшой размерности или при небольшой точности расчетов.
Некоторым аналогом замораживания матрицы Якоби является применение в расчетах алгоритмов интегрирования на основе явных и ¿-устойчивых методов с автоматическим выбором численной схемы [5]. В этом случае эффективность алгоритма может быть повышена за счет расчета переходного участка, соответствующего максимальному собственному числу матрицы Якоби, по явному методу. В качестве критерия выбора эффективной численной формулы естественно применять неравенство для контроля устойчивости [6-8]. Отметим, что применение таких комбинированных алгоритмов полностью не снимает проблемы замораживания матрицы Якоби [9].
Здесь построено неравенство для контроля устойчивости метода Ческино второго порядка точности. На основе стадий данной численной формулы построена схема первого порядка с расширенной до 32 единиц по вещественной оси областью устойчивости. На основе ¿-устойчивой (2,1)-схемы и рассмотренных явных численных формул разработан алгоритм переменной структуры, в котором эффективная численная формула выбирается на каждом шаге по критерию устойчивости. Алгоритм предназначен для решения жестких и не жестких задач с точностью расчетов порядка 1%. Приведены результаты вычислений, подтверждающие работоспособность и эффективность построенного алгоритма.
Для численного решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений
где у и / — вещественные Ж-мерные вектор-функции, £ — независимая переменная, в [10] предложен класс численных (т,к)-схем. С точки зрения реализации эти методы столь же просты, как и схемы типа Розенброка. Однако для них значительно проще решаются проблемы численного вычисления и замораживания матрицы Якоби. Кроме того, (т, к)-методы обладают более хорошими свойствами точности и устойчивости при незначительном увеличении вычислительных затрат. В традиционных методах число стадий т полностью описывает численную формулу. В (т, к)-методах для описания численных схем требуется две постоянные т — число стадий и к — количество вычислений правой части системы (1) на шаге интегрирования.
Для решения задачи (1) рассмотрим (2,1)-схему вида
1. (2,1)-МЕТОД РЕШЕНИЯ ЖЕСТКИХ ЗАДАЧ
у' = /Ы, У(£о) = уо, £0 < t < tk,
(1)
Уп+1 = Уп + Р1 кх + Р2к2, кх = / (уп), к2 = кх,
(2)
где к1 и к2 — стадии метода; Бп = Е — аЬАп; Ап — некоторая матрица, представимая в виде Ап = /П + ЬБп + О(Ь2); Е — единичная матрица; Ь — шаг интегрирования; Щ = д/(уп)/ду — матрица Якоби задачи (1); Бп — независящая от шага произвольная матрица; а, р1 и р2 — числовые коэффициенты. Такая формулировка метода позволяет применять (2) с замораживанием как аналитической, так и численной матрицы Якоби [11]. В случае использования в расчетах матрицы Якоби /'п-и, вычисленной к шагов назад, имеем Бп = —кЩ/п, Щ/п = д2/(уп)/ду2. Если матрица Якоби вычисляется численно с шагом г^ = с^Ь, то элементы матрицы Бп имеют вид Ьп^ = 0.5^д2/г(уп)/ду2. Получим коэффициенты ¿-устойчивой численной схемы (2) второго порядка и неравенство для контроля точности вычислений. Разложение точного решения у(Ьп+1) задачи (1) в ряд Тейлора имеет вид
у(1п+1) = у{Ьп) + ЬI + 2 Ь2 /' / + 1ь3 [I'21 + I'' 12 ] + О(Ь4),
где элементарные дифференциалы //'2/ и /"/2 вычислены на точном решении у(Ьп). Для нахождения коэффициентов а, р1 и р2 схемы (2) запишем разложения стадий к1 и к2 в ряды Тейлора и подставим в (2). Получим
уп+1 = уп + (Р1 + Р2)Ь/п + а(р1 + 2р2 )Ь2/'п /п + а2(р1 + 3р2 )Ь3/'2 /п + а(рц + 2р2)Ь3 Бп /п + О(Ь4),
где элементарные дифференциалы /п, /п, /'2/п, Ц и Бп/п вычислены на приближенном решении уп. Полагая уп = у(Ьп) и сравнивая полученные ряды до членов с Ь2, получим условия второго порядка точности схемы (2), т.е.
Р1 + Р2 = 1, а(р1 +2р2 ) = \, •
Исследуем устойчивость схемы (2). Применяя ее к задаче у' = Ху, Ие (Л) < 0, получим уп+1 = Я(х)уп, где х = ЬЛ, а функция устойчивости Q(x) имеет вид
П( ) = 1 + (р1 + р2 — 2а)х + а(а — р1)х2
Q(X) = (1—Щ2 .
Тогда схема (2) будет ¿-устойчивой, если р1 = а. Подставляя это соотношение в условия порядка, имеем:
р1 = а, р2 = 1 — а, (3)
где а определяется из условия ¿-устойчивости а2 — 2а + 0.5 = 0. Сравнивая ряды для точного и приближенного решений до членов с Ь3, получим, что локальная ошибка 5п численной схемы (2), (3) имеет вид
5п = Ь3 [(а — 3)/'2/ + 1 /''/2 — 2Бп/] + О(Ь4).
г2
а = а1, так как в этом случае меньше коэффициент в главном члене Ь3/'2/ локальной ошибки. Рассмотрим одновременно схему типа Розенброка с двумя вычислениями правой части задачи (1) на каждом шаге, т. е.
уп+1 = уп + р1к1 + р2к2, Б пк1 = Ь/(уп), Б пк2 = Ь/(уп + вк\). (4)
Согласно [4] при в = а коэффициенты (3) обеспечивают второй порядок точности схемы (4), а условие а2 — 2а + 0.5 = 0 — ее ¿-устойчивость. Из [1, 4] следует, что численная формула (4) с коэффициентами (3) является одной из наиболее удачных среди двухстадийных методов Розенброка. Локальная ошибка (4) имеет вид
Уравнение устойчивости а2 — 2а + 0.5 = 0 имеет два корня а1 = 1 — 0.5^2 и а2 = 1 + 0.5л/2. Выберем
¿Г = Ь3 [(а — 1)/'2/ + (6 + а)/''/2 — аБп/] + О(Ь4).
Построенная здесь схема (2), (3), так же как и (4) с коэффициентами (3), обладает вторым порядком точности и ¿-устойчивостью, а их локальные ошибки различаются незначительно. В то же время (2)
требует на каждом шаге на одно вычисление функции / меньше (4) при прочих равных затратах, что делает ее предпочтительнее.
Контроль точности вычислений численной схемы (2) построим по аналогии [4]. Для этого введем обозначение
где и вычисляются по формулам (2). Тогда согласно [4] для контроля точности на каждом шаге нужно проверять неравенство
НЭ)|| < е, 1 < Эп < 2,
где е — требуемая точность расчетов, || ■ || — некоторая норма в Ян, а переменная Эп выбирается наименьшей, при которой выполняется данное неравенство.
Отметим одну важную особенность построенной оценки ошибки г>(э'п). Схема (2) является ¿-устойчивой, т. е. для ее функции устойчивости ф(х) имеет место соотношение ф(х) ^ 0 при х ^ —го. Так как для точного решения у(Ьп+1 ) = ехр(х)у(Ьп) задачи у' = Ау, Ие (Л) < 0, выполняется аналогичное свойство, то естественным будет требование стремления к нулю оценки ошибки при х ^ —го. Однако для г>(1) это свойство не выполняется — данная оценка ведет себя А-устойчивым образом. С целью исправления асимптотического поведения ошибки вместо г>(1) введена оценка г>(э'п), 1 < Эп < 2. В этом случае поведение оценки ошибки при Эп = 2 будет согласовано с поведением точного решения тестовой задачи при х ^ —го. Подчеркнем, что в смысле главного члена оценки г>(1) и г>(2) совпадают. Использование г>(э'п) фактически не приводит к увеличению вычислительных затрат. Это связано с тем, что г>(э'п) при Эп = 2 проверяется только в том случае, если оно нарушено при Эп = 1. Такая ситуация встречается достаточно редко, в основном при быстром росте величины шага интегрирования. Однако это позволяет уменьшить количество неоправданных повторных вычислений решения (возвратов).
Оценку максимального собственного числа и>п,о = ЬАп,тах матрицы Якоби системы (1), необходимую для перехода на явную формулу, оценим по формуле
Шп,0 = Ь
д/
дУ
= Ь ■ /ду
Ниже данная оценка будет применяться для автоматического выбора численной схемы — явная или ¿-устойчивая численная формула.
2. МЕТОД ЧЕСКИНО
Для решения (1) рассмотрим явную формулу типа Рунге-Кутта вида
Уп+1 = Уп + Рт1 &1 + Рт2 ^2 + Рш3&3 + Рш4&4,
= ь/(Ьп ,Уп), ^2 = Ь/(Ьп + 0.25Ь,уп + 0.25^1), = Ь/(Ьп + 0.5Ь,уп + 0.5^), (5)
^4 = Ь/(Ьп + Ь, Уп + — 2^2 + 2^3),
где Ь — шаг интегрирования, 1 < г < 4, — стадии метода, 1 < г < 4, — числовые коэффициенты, т — порядок точности метода. При коэффициентах
Р21 = 1, Р22 = —2, Р23 = 2, Р24 = 0, (6)
схема (5), (6) имеет второй порядок точности [12]. Схема (5) с коэффициентами р41 = р44 = 1/6, р42 = 0 и р43 = 2/3 имеет четвертый порядок. Тогда для контроля точности схемы второго порядка можно использовать оценку ошибки <5п,2 вида
¿п,2 => (,Р4г — Р2г )&г . ^—'г=1
В результате для контроля точности вычислений применяется неравенство ||^п,2|| < е, где || ■ || — некоторая норма в Ян, е — требуемая точность расчетов. Учитывая, что имеет место соотношение
$п,2 = 0(Л3), шаг Лас по точности выбирается по формуле Лас = дЛ, где q находится из уравнения д3 ||5п,2|| = е. Если д < 1, то происходит повторное вычисление решения (возврат) с шагом Л, равным дЛ. В противном случае вычисляется приближенное решение, а прогнозируемый шаг Лп+1 вычисляется по формуле Лп+1 = дЛ. Неравенство ||5п,2|| < е хорошо зарекомендовало себя при решении многих практических задач, ниже оно будет использоваться здесь.
3. КОНТРОЛЬ УСТОЙЧИВОСТИ МЕТОДА ЧЕСКИНО
Построим неравенство для контроля устойчивости схемы (2). Для этого применим (2) для решения линейной задачи у' = Ау с постоянной матрицей А. Первые три стадии к1, к2 и к3 применительно к данной задаче имеют вид
к1 = Хуп, к2 = [х + 0.25Х2] у,п, кз = [х + 0.5Х2 + 0.125Х3] где X = ЛА. Нетрудно видеть, что имеют место соотношения
к1 - 2к2 + к3 = 0.125Х3уп, 0.5(^2 - ) = 0.125Х2уп.
Теперь оценку максимального собственного числа wn,i матрицы Якоби системы (1) можно вычислить степенным методом [7]. Введем обозначение
9 Г |(ki - 2k2 + кз)i| \ (7)
wn,i = 2 ■ ma^ -—----} . (7)
, i<i<N { |(k2 - ki)i| J
Тогда для контроля устойчивости метода Ческино можно применять неравенство wn,i < D, где число D ограничивает интервал устойчивости.
Устойчивость методов типа Рунге-Кутта обычно исследуется на скалярном тестовом уравнении y' = Ay, где А есть произвольное комплексное число, Re (А) < 0. Смысл А — некоторое собственное число матрицы Якоби задачи (1). Применяя (5), (6) для решения y' = Ay, получим, что функция устойчивости Q2(x) метода второго порядка точности имеет вид Q2(x) = 1 + x + 0.5x2 + 0.125x3, x = hA. Интервал устойчивости метода второго порядка равен двум, а метода четвертого порядка приблизительно равен 2.8. Поэтому в неравенстве wn,i < D положим D = 2. Учитывая, что wn,i = O(h), шаг hst по устойчивости можно выбирать по формуле hst = rh, где r вычисляется из равенства rwn,i = 2.
Оценка (7) является грубой, потому что максимальное собственное число не обязательно сильно отделено от остальных, в степенном методе применяется мало итераций и дополнительные искажения вносит нелинейность задачи (1). Поэтому контроль устойчивости используется как ограничитель на размер шага интегрирования. В результате прогнозируемый шаг вычисляется по формуле
hn+i = max |hn, min(hac,hst)J, (8)
где hn есть последний успешный шаг интегрирования.
Формула (8) применяется для прогноза шага интегрирования hn+i после успешного вычисления решения с предыдущим шагом hn, и поэтому фактически не приводит к увеличению вычислительных затрат. Если шаг по устойчивости меньше последнего успешного, то он уменьшен не будет, потому что причиной этого может быть грубость оценки максимального собственного числа. Однако шаг не будет и увеличен, потому что не исключена возможность неустойчивости численной схемы. Формула (8) позволяет стабилизировать поведение шага на участке установления решения, где определяющую роль играет устойчивость.
Из результатов расчетов алгоритмом интегрирования с контролем точности и дополнительным контролем устойчивости следует, что фактическая точность вычисления решения на интервале установления значительно выше задаваемой. Это естественно, потому что старые ошибки подавляются за счет контроля устойчивости, а новые ошибки невелики за счет малости производных решения. В такой ситуации эффективнее проводить вычисления методом более низкого порядка точности с широкой областью устойчивости.
4. МЕТОД ПЕРВОГО ПОРЯДКА ТОЧНОСТИ
Для численного решения задачи (1) рассмотрим схему вида
Уп+1 = Уп + Гхкх + Г2 + Гз кз + Г4&4, (9)
где стадии 1 < г < 4, вычислены по формулам (5), а коэффициенты гг, 1 < г < 4, подлежат определению. Заметим, что при г1 = 1, г2 = -2, г3 = 2 и г4 = 0 численная формула (9) имеет второй порядок точности и совпадает с (5) с коэффициентами (6). Построим менее точную схему с максимальным интервалом устойчивости. Для этого применим (9) для решения скалярного тестового уравнения у' = Ау. Получим уп+1 = ^(х)уп, где функция устойчивости ф(х) имеет вид
ф(х) = 1 + (г1 + г 2 + г 3 + г 4 )х + (0.25г2 + 0.5г3 + г4)х2 + (0.125г3 + 0.5г4)х3 + 0.25г4 х4.
Требование первого порядка точности приводит к соотношению г1 + г2 + г3 + г4 = 1, которое ниже будем считать выполненным. Оставшиеся коэффициенты гг выберем таким образом, чтобы метод (9) имел максимальный интервал устойчивости. Для этого рассмотрим многочлен Чебышева Т4(г) = 8г4 — 8х2 + 1 на промежутке [-1,1]. Проведем замену переменных, полагая г = 1 — 2х/7. Получим:
Т4 (х) = 1 — 32х/7 + 160х2 /72 — 256х3/73 + 128х4 /74,
при этом отрезок [7,0] отображается на [—1,1]. Нетрудно показать, что среди всех многочленов Р4(х) вида Р4(х) = 1 + х + с2х2 + с3х3 + с4х4 для Т4(х) неравенство |Т4(х)| < 1 выполняется на максимальном интервале [7,0], 7 = —32. Потребуем совпадения коэффициентов ^(х) и Т4(х) при 7 = —32. Это приводит к соотношениям
г 1 + г 2 + г 3 + г 4 = 1, 0.25г2 + 0.5г3 + г4 = 5/32, 0.125г3 + 0.5г4 = 1/128, 0.25г4 = 18192.
В результате имеем коэффициенты
г 1 = 895/2048, г2 = 257/512, г3 = 31/512, г4 = 1/2048
метода первого порядка точности с максимальным интервалом устойчивости, локальная ошибка 5п которого имеет вид 5п = 9Л2/'//32+0(Л,3). Для контроля точности численной формулы первого порядка будем использовать оценку локальной ошибки. Учитывая, что имеет место к2—к1 = 0.25Л2//п+0(^3) и вид локальной ошибки, неравенство для контроля точности записывается в виде ||к2 — к^Ц < е, где || ■ У — некоторая норма в Ян, е — требуемая точность расчетов.
Интервал устойчивости численной схемы (9) равен 32. Поэтому для ее контроля устойчивости можно применять неравенство и>пд < 32, где и>пд задается (7).
5. АЛГОРИТМ ИНТЕГРИРОВАНИЯ ПЕРЕМЕННОЙ СТРУКТУРЫ
На основе построенных явных методов первого и второго порядков точности легко сформулировать алгоритм переменного порядка и шага. Расчеты всегда начинаются методом второго порядка как более точным. Переход на схему первого порядка осуществляется при нарушении неравенства ш^д < 2. Обратный переход на метод второго порядка происходит в случае выполнения неравенства и>пд < 2. При расчетах по методу первого порядка, наряду с точностью, контролируется устойчивость, а выбор прогнозируемого шага производится по аналогии с методом второго порядка точности по формуле типа (8).
В случае использования схемы (2) формулировка алгоритма интегрирования также не вызывает трудностей. Нарушение неравенства и>пд < 32 вызывает переход на схему (2). Передача управления явным методам происходит в случае выполнения неравенства и>п,о < 32, где оценка и>п,о вычисляется через норму матрицы Якоби.
Численную формулу (2) без потери порядка точности можно применять с замораживанием матрицы Дп. Отметим, что при замораживании матрицы Якоби величина шага интегрирования остается постоянной с целью сохранения свойства ¿-устойчивости метода (2). Попытка замораживания матрицы Бп осуществляется после каждого успешного шага. Размораживание матрицы происходит в следующих случаях: 1) нарушение точности расчетов, 2) если число шагов с замороженной матрицей достигло заданного максимального числа г^, 3) если прогнозируемый шаг больше последнего
успешного в qh раз. Числами ¿ь и qh можно влиять на перераспределение вычислительных затрат. При ¿ь = 0 и qh = 0 замораживания не происходит, при увеличении ¿ь и qh число вычислений правой части возрастает, а количество обращений матрицы Якоби убывает. В расчетах использовались значения ¿ь = 10 и qh = 2. Норма ||^|| в неравенствах для контроля точности вычисляются по формуле
где i — номер компоненты, r — положительный параметр. Если по i-й компоненте решения выполняется неравенство |уП | < r, то контролируется абсолютная ошибка re, в противном случае — относительная ошибка е.
Ниже построенный алгоритм переменного порядка и шага, а также с автоматическим выбором явной или L-устойчивой численной схемы будем называть MK2CES24.
6. РЕЗУЛЬТАТЫ РАСЧЕТОВ
Расчеты проводились с задаваемой точностью е = 10-2 на PC Intel(R) Core i7-3770S CPU@3.10GHz с двойной точностью. Невысокая точность расчетов связана с тем, что в построенном алгоритме применяются схемы низкого порядка точности, и поэтому данным методом осуществлять расчеты с более высокой точностью нецелесообразно. Сравнение эффективности проводилось с известным методом Гира в реализации А. Хиндмарша DLSODE из коллекции ODEPACK [13] и алгоритмом RKMK2 [11] на основе L-устойчивого и явных двухстадийных методов. В качестве тестовой задачи выбрана модель реакции Белоусова-Жаботинского, для которой на промежутке интегрирования характерно наличие трех переходных процессов. Тестовый пример имеет вид
Расчеты проводились с численной матрицей Якоби. Сравнение эффективности проводилось по числу вычислений правой части и матрицы Якоби уа задачи (10) на интервале интегрирования. Решение (10) алгоритмом МК2СЕ824 вычислено с затратами ¡!и = 1 029 и уа = 49. При вычислениях по ¿-устойчивой схеме (2) затраты ¡!и = 926 и уа = 88. Фактическая точность расчетов в конце интервала интегрирования не хуже задаваемой. Решение (10) удалось вычислить явными методами переменного порядка и шага с затратами ¡!и = 978 524. Данная задача слишком жесткая для явных методов. Однако результаты расчетов приведены здесь с целью демонстрации принципиальной возможности применения явных методов с контролем устойчивости и переменным порядком для решения достаточно жестких примеров. Явные методы на некоторых жестких задачах большой размерности могут быть эффективнее ¿-устойчивых методов. Решение данной задачи алгоритмом РКМК2 [11] вычислено с затратами ¡!и = 1 214 и уа = 65. Решение (10) явными методами переменного порядка и шага из алгоритма РКМК2 вычислено с затратами ¡!и = 2 112 678.
При расчетах программой ЭЬЗОЭЕ требуемая точность е = 10-2 достигается при задаваемой точности е = 10-4 с затратами ¡!и = 1 129 и уа = 107. При более высокой точности расчетов ЭЬЗОЭЕ эффективнее построенного алгоритма. Это является следствием низкого порядка точности построенных численных формул.
При задаваемой точности е = 10-2 алгоритм МК2СЕБ24 более чем в 2 раз эффективнее алгоритма ЭЬЗОЭЕ по числу декомпозиций матрицы Якоби, в то время как количество вычислений правой части (10) для МК2СЕБ24 и ЭЬЗОЭЕ различается незначительно. В случае большой размерности задачи (1) построенный алгоритм интегрирования по времени счета может быть эффективнее ЭЬЗОЭЕ.
Построенный алгоритм МК2СЕБ42 предназначен для расчетов с небольшой точностью — порядка 1% и ниже. В этом случае достигается его максимальная эффективность. В МК2СЕБ42 с помощью признака можно задавать различные режимы: 1) явными методами первого или второго порядков
IMI = max { ^ |/( | + r)},
y1 = 77.27(y2 - У1У2 + У1 - 8.375 ■ 10-6y2), у2 = (-У2 - У1У2 + У3)/77.27, уЗ = °.161(У1 - Уз), t е [0,300], У1 (0) = 4, У2(°) = 1.1, Уз(°) = 4, ho = 2 ■ 10-3.
(10)
ЗАКЛЮЧЕНИЕ
точности с контролем или без контроля устойчивости; 2) явными методами с переменным порядком и шагом; 3) ¿-устойчивым методом с замораживанием или без замораживания аналитической или численной матрицы Якоби; 4) с автоматическим выбором численной схемы. Все это позволяет применять данный алгоритм для решения как жестких, так и нежестких задач. При расчетах с автоматическим выбором численной схемы вопрос о том, является ли задача жесткой или нет, перекладывается на алгоритм интегрирования.
Использование неравенства для контроля устойчивости фактически не приводит к увеличению вычислительных затрат, потому что оценка максимального собственного числа матрицы Якоби системы (1) осуществляется через ранее вычисленные стадии и не приводит к росту числа вычислений функции /. Такая оценка получается грубой. Однако применение контроля устойчивости в качестве ограничителя на рост шага позволяет избежать негативных последствий грубости оценки. Более того, в некоторых случаях это приводит к нестандартно высокому повышению эффективности алгоритма. На участке установления за счет контроля устойчивости старые ошибки стремятся к нулю, а новые невелики за счет малости производных решения. В некоторых случаях вместо оценки максимального собственного числа оценивается следующее по порядку. Шаг интегрирования становится больше максимально допустимого, и с таким шагом осуществляется интегрирование до тех пор, пока не нарушается неравенство для контроля точности. Как правило, число таких шагов невелико. Однако величина шага может на порядок превышать максимальный шаг по устойчивости. После нарушения неравенства для контроля точности шаг уменьшается до максимально возможного. Такой эффект может повторятся многократно в зависимости от длины участка установления. В результате средний шаг интегрирования может превышать максимально допустимый.
Применение на участке установления явного метода первого порядка точности с расширенной областью устойчивости позволяет в 16 раза увеличить размер шага интегрирования по сравнению с явным методом второго порядка без увеличения вычислительных затрат. На переходных участках, где определяющую роль играет точность вычислений, более эффективным является метод второго порядка точности, хотя и с небольшой областью устойчивости. Комбинирование методов низкого и высокого порядков с помощью неравенства для контроля устойчивости позволяет повысить эффективность расчетов.
Работа выполнена при финансовой поддержке РФФИ (проекты 11-01-00106 и 11-01-00224).
Библиографический список
1. Хайрер Э, Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М. : Мир, 1999. 685 с.
2. Byrne G. D., Hindmarsh A. C. ODE solvers : a review of current and coming attractions // J. of Comput. Physics. 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. Новиков Е. А. Алгоритм интегрирования жестких задач с помощью явных и неявных методов // Изв. Сарат. ун-та. Нов. сер. Сер. Математика. Механика. Информатика. 2012. Т. 12, вып. 4. С. 19-27.
7. Новиков В. А., Новиков Е. А. Повышение эффектив-
ности алгоритмов интегрирования обыкновенных дифференциальных уравнений за счет контроля устойчивости // ЖВМ и МФ. 1985. Т. 25, № 7. С. 1023-1030.
8. Новиков Е. А. Явные методы для жестких систем. Новосибирск : Наука, 1997. 197 с.
9. Новиков Е. А., Шорников Ю. В. Компьютерное моделирование жестких гибридных систем. Новосибирск : Изд-во НГТУ, 2012. 450 с.
10. Новиков Е. А., Шитов Ю. А., Шокин Ю. И. Од-ношаговые безытерационные методы решения жестких систем // Докл. АН СССР. 1988. Т. 301, № 6. С. 13101314.
11. Новиков A. E., Новиков E. A. Численное решение жестких задач с небольшой точностью // Математическое моделирование. 2010. Т. 22, № 1. С. 46-56.
12. Ceschino F, Kuntzman J. Numerical solution of initial value problems. New Jersey : Prentice-Hall, Englewood Clis, 1966. 287 p.
13. Hindmarsh A.C. ODEPACK, a systematized collection of ODE solvers // Lawrence Livermore National Laboratory, 1982. Preprint UCRL-88007.
Algorithm Variable Order, Step and the Configuration Variables for Solving Stiff Problems
E. A. Novikov
Institute of Computational Modelling, Siberian Branch of the Russian Academy of Sciences, Russia, 660036, Krasnoyarsk, Akademgorodok, novikov@icm.krasn.ru
An inequality for stability control of a Ceschino's scheme of second order of accuracy is constructed. A numerical formula of order one is developed that is based on the stages of the this method and its stability interval is extended to 32. On a base of L-stable (2,1)-scheme and a numerical Ceschino's formula, an algorithm of alternating structure, in which an efficient numerical formula is chosen on an every step by a stability criterion, is constructed. The algorithm is intended for solving stiff and non-stiff problems. There are shown results of calculations, confirming efficiency of this algorithm.
Key words: stiff problem, Ceschino's scheme, (2,1)-method, accuracy and stability control.
References
1. Hairer E., Wanner G. Solving ordinary differential equations II. Stiff and differential-Algebraic problems. New York, Springer-Verlag, 1996, 601 p.
2. Byrne G. D., Hindmarsh A. C. ODE solvers : a review of current and coming attractions. J. of Comput. Physics, 1987, no. 70, pp. 1-62.
3. Rosenbrock H. H. Some general implicit processes for the numerical solution of differential equations. Computer, 1963, no. 5, pp. 329-330.
4. Novikov V. A., Novikov E. A., Yumatova L. A. Freezing of a matrix of Jacobi in the Rosenbrock method of the second order of accuracy. Zhurnal Vychislitel'noi Matematiki i Matematicheskoi Fiziki, 1987, vol. 27, no. 3, pp. 385-390 (in Russian).
5. Novikov E. A. Construction of algorithm for the integrating stiff differential equations on nonuniform schemes. Soviet Math. Dokl., 1984, vol. 30, no. 2, pp. 358361.
6. Novikov E. A. Algorithm of Integrating Stiff Problems Using the Explicit and Implicit Methods. Izv. Sarat. Univ. N.S. Ser. Math. Mech. Inform., 2012, vol. 12, no. 4, pp.19-27 (in Russian).
7. Novikov V. A., Novikov E. A. Increase of efficiency YAK 514.133
of algorithms of integration of the ordinary differential equations at the expense of stability control. Zhurnal Vychislitel'noi Matematiki i Matematicheskoi Fiziki, 1985, vol. 25, no. 7, pp. 1023-1030 (in Russian).
8. Novikov E. A. Explicit methods for stiff systems. Novosibirsk, Nauka, 1997, 197 p. (in Russian).
9. Novikov E. A., Shornikov Yu. V. Computer modeling of stiff hybrid systems. Novosibirsk, publishing house NGTU, 2012, 450 p. (in Russian).
10. Novikov E. A., Shitov Yu. A., Shokin Yu. I. One-step iteration-free methods of solving stiff systems. Soviet Math. Dokl., 1989, vol. 38, no. 1, pp. 212-216.
11. Novikov A. E., Novikov E. A. Numerical integration of stiff systems with low accuracy. Mathematical Models and Computer Simulations, 2010, vol. 2, no. 4, pp. 443452. DOI: 10.1134/S2070048210040046.
12. Ceschino F., Kuntzman J. Numerical solution of initial value problems. New Jersey: Prentice-Hall, Englewood Clis, 1966, 287 p.
13. Hindmarsh A. C. ODEPACK, a systematized collection of ODE solvers. Lawrence Livermore National Laboratory, 1982, preprint UCRL-88007.
ГИПЕРБОЛИЧЕСКИЕ ПАРАЛЛЕЛОГРАММЫ ПЛОСКОСТИ Я
Л. Н. Ромакина
Кандидат физико-математических наук, доцент кафедры геометрии, Саратовский государственный университет им. Н. Г. Чернышевского, romakinaln@mail.ru
На гиперболической плоскости Н положительной кривизны в модели Кэли-Клейна исследованы гиперболические параллелограммы. Проведена их классификация, получены метрические соотношения между величинами углов и выражения длин ребер через меры углов при вершинах.
Ключевые слова: гиперболическая плоскость Н положительной кривизны; параллелограмм; гиперболический параллелограмм.
ВВЕДЕНИЕ
Гиперболическую плоскость Н положительной кривизны [1-5] рассматриваем в проективной интерпретации Кэли-Клейна как внешнюю относительно овальной линии 7, называемой абсолютом плоскости Н, область проективной плоскости Р2.
