Параметризация задачи Коши для нелинейных дифференциальных уравнений с контрастными структурами Текст научной статьи по специальности «Математика»

ВЕСТНИК МОРДОВСКОГО УНИВЕРСИТЕТА

Том 28, № 4. 2018

http://vestnik.mrsu.ru

ISSN Print 0236-2910 ISSN Online 2313-0636

УДК 517.91

DOI: 10.15507/0236-2910.028.201804.486-510

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

Е. Б. Кузнецов*, С. С. Леонов, Е. Д. Цапко

ФГБОУ ВО «Московский авиационный институт (национальный исследовательский университет)» (г. Москва, Россия) *[email protected]

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

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

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

Для цитирования: Кузнецов Е. Б., Леонов С. С., Цапко Е. Д. Параметризация задачи Коши для нелинейных дифференциальных уравнений с контрастными струк-

© Кузнецов Е. Б., Леонов С. С., Цапко Е. Д., 2018

@ I Контент доступен по лицензии Creative Commons Attribution 4.0 License. ES^H This work is licensed under a Creative Commons Attribution 4.0 License.

турами // Вестник Мордовского университета. 2018. Т. 28, № 4. С. 486-510. DOI: https://doi.org/10.15507/0236-2910.028.201804.486-510

Благодарности: Исследование проведено при финансовой поддержке Российского научного фонда (код проекта 18-19-00474).

The Parametrization of the Cauchy Problem for Nonlinear Differential Equations with Contrast Structures

E. B. Kuznetsov*, S. S. Leonov, E. D. Tsapko

Moscow Aviation Institute (National Research University) (Moscow, Russia)

*[email protected]

Introduction. The paper provides an analysis of numerical methods for solving the Cauchy problem for nonlinear ordinary differential equations with contrast structures (interior layers). Similar equations simulate various applied problems of hydro- and aeromechanics, chemical kinetics, the theory of catalytic reactions, etc. An analytical solution to these problems is rarely obtained, and numerical procedure is related with significant difficulties associated with ill-conditionality in the neighborhoods of the boundary and interior layers. The aim of the paper is the scope analysis of traditional numerical methods for solving this class problems and approbation of alternative solution methods.

Materials and Methods. The traditional explicit Euler and fourth-order Runge-Kutta methods, as well as the implicit Euler method with constant and variable step sizes are used for the numerical solution of the Cauchy problem. The method of solution continuation with respect to the best argument is suggested as an alternative to use. The solution continuation method consists in replacing the original argument of the problem with a new one, measured along the integral curve of the problem. The transformation to the best argument allows obtaining the best conditioned Cauchy problem. Results. The computational difficulties arising when solving the equations with contrast structures by traditional explicit and implicit methods are shown on the example of the test problem solution. These difficulties are expressed in a significant decrease of the step size in the neighborhood of the boundary and interior layers. It leads to the increase of the computational time, as well as to the complication of the solving process for super stiff problems. The authenticity of the obtained results is confirmed by the comparison with the analytical solution and the works of other authors.

Conclusions. The results of the computational experiment demonstrate the applicability of the traditional methods for solving the Cauchy problem for equations with contrast structures only at low stiffness. In other cases these methods are ineffective. It is shown that the method of solution continuation with respect to the best argument allows eliminating most of the disadvantages inherent to the original problem. It is reflected in decreasing the computational time and in increasing the solution accuracy.

Keywords: contrast structures, method of solution continuation, the best argument, ill-conditionality, the Cauchy problem, ordinary differential equation

For citation: Kuznetsov E. B., Leonov S. S., Tsapko E. D. The Parametrization of the Cauchy Problem for Nonlinear Differential Equations with Contrast Structures. Vestnik Mordovskogo universiteta = Mordovia University Bulletin. 2018; 28(4):486-510. DOI: https://doi.org/10.15507/0236-2910.028.201804.486-510

Acknowledgements: This work was supported by the Russian Science Foundation, project no. 18-19-00474.

Введение

В статье рассмотрена задача Коши для обыкновенного дифференциального уравнения п-го порядка, содержащего малый параметр £ при старшей производной [1]:

( _7 _7П-1 \

а у ^ ау а у Л £ —- = F t, у, — ,...,-f , е> 0 (1)

dtn ^ аГ аг1)

с начальными условиями

У{0) = = ,к =п -(2)

В работах А. Н. Тихонова уравнения вида (1) получили название уравнений с малым параметром при старшей производной. Впоследствии они стали известны как сингулярно возмущенные уравнения (в практических задачах их называют уравнениями с пограничным слоем).

Характерной особенностью решений задачи (1)-(2) является наличие областей резкого изменения решения в окрестности начальной (или, в случае краевой задачи, граничной) точки. Одним из первых эффект пограничного слоя отметил в своих работах Л. Прандтль, рассматривавший движение вязкой жидкости с малым трением, описываемое уравнениями Навье-Стокса1. Впоследствии сингулярно возмущенные уравнения нашли множество применений в гидроаэромеханике, химической кинетике, теории каталитических реакций, дифференциальной геометрии, проектировании атомных реакторов и т. д. С середины 70-х гг. ХХ в. помимо задач с пограничными слоями активно изучаются

задачи с внутренними слоями2, получившими в работах А. Б. Васильевой, В. Ф. Бутузова и Н. Н. Нефедова название контрастных структур [2].

Несмотря на широкую область применения, уравнения с пограничным слоем и контрастными структурами остаются для исследователей трудоемкими в использовании. Если функция правой части уравнения (1) является нелинейной, то найти точное аналитическое решение задачи (1)-(2) удается лишь в исключительных случаях. Поэтому приходится использовать численные и приближенно-аналитические (асимптотические) методы решения. Стоит отметить, что явные методы численного решения задачи Коши3 при £ ^ 0 требуют чрезмерного уменьшения шага интегрирования для получения решения высокой точности, что делает их малоэффективными для данного класса задач.

Кроме явных методов, существует большой класс специализированных методов решения жестких задач, к которым относится и задача (1)-(2). Им посвящена монография Э. Хай-рера и Г. Ваннера4. Однако и специализированные методы чаще всего не позволяют устранить все недостатки, присущие явным методам. Более того, поскольку многие методы решения жестких задач основаны на применении неявных схем, возникает проблема решения нелинейных алгебраических или трансцендентных уравнений и их систем, что сопряжено с построением итерационной последовательности и доказательством ее сходимости из выбранного начального приближения. Все это усложняет вычислительный

1 Прандтль Л. Теория несущего крыла. Ч. 1. Движение жидкости с очень малым трением. М. ; Л. : ГНТИ, 1931. С. 5-11.

2 Чанг К., Хауэс Ф. Нелинейные сингулярно возмущенные краевые задачи : теория и приложения. М. : Мир, 1988. 247 с.

3 Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. М. : Мир, 1990. 512 с.

4 Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М. : Мир, 1999. 685 с.

процесс. Также неявные методы могут вызывать затруднения при переходе через внутренние слои.

Другой класс методов заключается в представлении решения задачи (1)-(2) рядом по степеням малого параметра с. К этим методам относятся алгоритм А. Б. Васильевой5-6 и метод регуляризации, предложенный С. А. Ломовым7 [9]. Указанные методы широко известны и востребованы. Однако стоит обратить внимание, что построение асимптотических разложений требует сложных аналитических выкладок. При малых с, т. е. при с ^ 0, для достижения необходимой точности достаточно учитывать лишь несколько слагаемых, но, если это условие не выполняется, может понадобиться учет десятков слагаемых. Таким образом, если для сверхжестких задач асимптотические методы дают хороший результат, то для некоторых классов плохо обусловленных задач их применение может вызывать затруднения.

Все это говорит о потребности в новых эффективных методах и подходах, и такие методы разрабатываются. Начиная со второй половины 90-х гг. ХХ в. в работах Э. И. Григолюка, В. И. Шала-шилина, Е. Б. Кузнецова и их учеников развивается метод продолжения решения по наилучшему аргументу8, известный также как метод длины дуги. Он состоит в замене исходного аргумента на новый, отсчитываемый по касательной вдоль интегральной кривой рассматриваемой задачи. Доказано, что указанный аргумент доставляет исходной задаче наилучшую обусловленность, поэтому он получил название наилучшего. Применительно к рассматриваемому классу задач новый аргумент позволяет обходить трудности, связанные

с плохой обусловленностью исходной задачи, в том числе уменьшить показатель жесткости. Это дает возможность получить вычислительные преимущества при численном решении. В статье А. А. Белова и Н. Н. Калиткина [3] рассматривается ряд начальных задач с контрастными структурами. Они решаются при помощи преобразования к наилучшему аргументу специальными методами с переменным шагом интегрирования, автоматически изменяемым по кривизне интегральной кривой [4]. В другой их работе полученные результаты дополняются сопоставлением с правилом Рунге - Ромберга - Ричардсона [5]. Однако ни в одной из указанных статей нет сравнения с решением исходных непреобразованных задач, хотя и дается ряд рекомендаций относительно методов их решения.

Целью данной статьи является сравнительный анализ решения тестовой начальной задачи с контрастными структурами традиционными явными и неявными методами разного порядка точности, а также методом продолжения решения по наилучшему аргументу. Важным является аргументированное определение границ применимости традиционных методов к решению задач с контрастными структурами. Полученные в данной работе результаты согласуются с выводами, сделанными в статьях А. А. Белова и Н. Н. Калиткина [3; 5], дополняя их.

Обзор литературы

Теория сингулярно возмущенных уравнений берет начало в работах А. Н. Тихонова [1; 6-7]. В них впервые был рассмотрен общий вид нелинейных уравнений и систем с малым

5 Васильева А. Б., Бутузов В. Ф. Асимптотические разложения решений сингулярно возмущенных уравнений. М. : Наука, 1973. 272 с.

6 Васильева А. Б., Бутузов В. Ф. Асимптотические методы в теории сингулярных возмущений. М. : Высшая школа, 1990. 208 с.

7 Ломов С. А. Введение в общую теорию сингулярных возмущений. М. : Наука, 1981. 400 с.

8 Шалашилин В. И., Кузнецов Е. Б. Метод продолжения решения по параметру и наилучшая параметризация в прикладной математике и механике. М. : Эдиториал УРСС, 1999. 224 с.

параметром при старшей производной, дано определение области влияния решения вырожденного уравнения и вырожденной системы, а также доказаны первые общие утверждения о близости решения вырожденного уравнения (системы уравнений) к решению исходной задачи. Полученные результаты были развиты в работах ученицы А. Н. Тихонова А. Б. Васильевой. Совместно с В. Ф. Бутузовым и другими ее учениками и коллегами были получены фундаментальные результаты по представлению решений сингулярно возмущенных задач асимптотическими рядами специальных видов.

Помимо уже упомянутых моногра-фий9-10, стоит отметить более поздние работы А. Б. Васильевой и А. А. Плот-никова11, а также В. Ф. Бутузова12. В конце 80-х - начале 90-х гг. ХХ в., помимо задач с пограничными слоями, в работах А. Б. Васильевой, В. Ф. Бутузова и Н. Н. Нефедова рассмотрены задачи с контрастными структурами (внутренними слоями). Основополагающей работой в этом направлении стала статья «Контрастные структуры в сингулярно возмущенных задачах» [2]. Обзор результатов, связанных с решением задач с контрастными структурами, дан в статье «Асимптотическая теория контрастных структур (обзор)» [8].

В последние пять лет новые результаты в области асимптотических методов решения сингулярно возмущенных задач получены во многих работах. Наиболее важные результаты получены представителями школы А. Н. Тихонова и А. Б. Васильевой. В работах В. Ф. Бутузова и соавторов с использованием метода дифференциальных неравенств рассматривается построе-

ние асимптотики для сингулярно возмущенных краевых задач (Дирихле и Неймана) для эллиптических уравнений с внутренними и угловыми пограничными слоями [9-10], а также для первой начально-краевой задачи для сингулярно возмущенного параболического уравнения [12]. Кроме того, рассмотрены краевые задачи для обыкновенных дифференциальных уравнений с многозонным внутренним слоем (контрастными структурами типа «всплеск» и «ступенька») и построены асимптотические разложения для них [13]. В указанных работах В. Ф. Бу-тузова рассматриваются случаи кратных корней вырожденного уравнения. В данном случае классический алгоритм А. Б. Васильевой не применим. Предложена модификация метода построения асимптотических разложений в случае двухкратного и трехкратного корней. К этой тематике можно отнести и статью М. В. Козлова и В. Н. Щенникова, посвященную асимптотической устойчивости вырожденного решения сингулярно возмущенных систем обыкновенных дифференциальных уравнений с однородной правой частью [14].

В работах Н. Н. Нефедова и соавт. рассматривается решение начально-краевых задач для уравнения реакции-диффузии-адвекции с контрастными структурами. Доказательство условий существования решения ведется методом дифференциальных неравенств. В статье «Существование и асимптотическая устойчивость периодического решения с внутренним переходным слоем в задаче со слабой линейной адвекцией» рассмотрено уравнение реакции-диффузии-адвекции со слабой адвекцией и реактивным членом

9 Васильева А. Б., Бутузов В. Ф. Асимптотические разложения решений сингулярно возмущенных уравнений. М. : Наука, 1973. 272 с.

10 Васильева А. Б., Бутузов В. Ф. Асимптотические методы в теории сингулярных возмущений. М. : Высшая школа, 1990. 208 с.

11 Васильева А. Б., Плотников А. А. Асимптотическая теория сингулярно возмущенных задач. М. : Физический факультет МГУ, 2008. 398 с.

12 Бутузов В. Ф. Асимптотические методы в сингулярно возмущенных задачах. Ярославль : Изд-во ЯрГУ, 2014. 140 с.

с кубической нелинейностью, даны достаточные условия существования решения [15]. При наличии большого адвективного слагаемого уравнение реакции-диффузии-адвекции рассмотрено в работе «Асимптотическое приближение решения уравнения реакция-диффузия-адвекция с нелинейным адвективным слагаемым» [16]; исследовано решение вида движущегося фронта. Существование и устойчивость решения краевой задачи для многомерного стационарного уравнения реакции-диффузии-адвекции с контрастными структурами рассмотрено в статье «Существование и устойчивость контрастных структур в многомерных задачах реакция-диффузия-адвекция в случае сбалансированной нелинейности» [17]. В работе «Асимптотика движения фронта в задаче реакция-диффузия-адвекция» [18] для начально-краевой задачи для уравнения реакции-диффузии-адвекции с граничными и внутренними слоями построено асимптотическое разложение решения с движущимся фронтом. Общая схема построения асимптотических разложений для задач с контрастными структурами и исследование их устойчивости дано в статье «Общая схема асимптотического исследования устойчивых контрастных структур» [19].

Стоит отметить, что на сегодняшний день аппарат дифференциальных неравенств играет одну из ведущих ролей при исследовании существования решения сингулярно возмущенных уравнений. Впервые дифференциальные неравенства были предложены в 1919-1920-м гг. С. А. Чаплыгиным13 для приближенного интегрирования систем дифференциальных уравнений.

В теории сингулярно возмущенных уравнений дифференциальные неравенства впервые использованы М. На-гумо14. Кроме указанных работ, значительное место в исследовании решений задач с пограничными и внутренними слоями занимает уже упомянутая монография К. Чанга и Ф. Хауэса15, в которой изложена основа аппарата дифференциальных неравенств применительно к доказательству существования решений нелинейных сингулярно возмущенных краевых задач, а также рассмотрено множество практических задач математической физики и вычислительной геометрии.

Отметим также вклад в теорию сингулярных возмущений С. А. Ломова и И. С. Ломова, в монографии которых дана математическая теория пограничного слоя для линейных дифференциальных уравнений в одномерном и многомерном случаях для операторов с различными свойствами16. В уже упомянутой монографии С. А. Ломова помимо линейных дифференциальных уравнений рассмотрены также некоторые классы нелинейных уравнений17.

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

13 Чаплыгин С. А. Новый метод приближенного интегрирования дифференциальных уравнений. М. ; Л. : Гостехиздат, 1950. 103 с.

14 Nagumo M. Uber das Verhalten der Integrals von fur ky" + f (x, y, y', k) = 0 fur k ^ 0 // Proceedings of the Physico-Mathematical Society of Japan. 1939. No. 21. P. 529-534.

15 Чанг К., Хауэс Ф. Нелинейные сингулярно возмущенные краевые задачи : теория и приложения. М. : Мир, 1988. 247 с.

16 Ломов С. А., Ломов И. С. Основы математической теории пограничного слоя. М. : Изд-во Москов. ун-та, 2011. 456 с. URL: http://www.rfbr.ru/rffi/ru/books/o_72182

17 Ломов С. А. Введение в общую теорию сингулярных возмущений. М. : Наука, 1981. 400 с.

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

Большинство работ по численному решению жестких и плохо обусловленных задач не рассматривают отдельно класс задач с контрастными структурами. Одними из немногих работ, посвещенных численному решению начальных задач с контрастными структурами, являются статьи Н. Н. Калиткина и А. А. Белова [3; 5]. В этих работах для решения тестовых задач с контрастными структурами предложено использовать комбинацию метода длины дуги со специальными методами на неравномерных сетках при оценке локальной погрешности по правилу Рунге и кривизне в точке интегральной кривой [4]. Полученные результаты применяются к задачам химической кинетики [25].

Идея метода продолжения решения по параметру используется в математике и механике еще в работах А. Пуанкаре и У. Леверье. По сути, метод замены переменной под знаком интеграла и метод малого параметра являются отражением этой общей идеи. В вычислительных целях метод продолжения решения по параметру был впервые использован в работах бельгийского математика М. Лаэя для построения решений трансцендентных уравнений и систем [26-27]. Позднее в работах советского математика Д. Ф. Дави-денко был предложен метод сведения системы нелинейных уравнений к си-

стеме обыкновенных дифференциальных уравнений [28-29]. Однако при построении замкнутых кривых приходится менять параметр продолжения решения в окрестности предельных особых точек. Для решения этой проблемы было предложено использовать параметры продолжения общего вида.

С целью решения нелинейных краевых задач теории упругости в работе И. И. Воровича и В. Ф. Зипаловой был предложен параметр продолжения, отсчитываемый по касательной к кривой множества решений рассматриваемой задачи [30]. В этой же работе впервые показано, что данный параметр продолжения будет в некотором смысле наилучшим для задач, кривые множества решений которых имеют предельные особые точки. Доказательство этого факта было намечено в статье голландского инженера Э. Рикса [31] применительно к исследованию устойчивости прощелкивающихся и выпучивающихся конструкций.

С конца 70-х гг. XX в. началось систематическое исследование применения наилучшего параметра к задачам математики и механики с предельными особыми точками и точками бифуркации. Стоит отметить монографии Э. И. Григолюка и В. И. Шалашилина18, а также К. Георга и Ю. Л. Аллгауэра19, подводящие итог полученным результатам. В уже упомянутой монографии В. И. Шалашилина и Е. Б. Кузнецова20 было доказано, что аргумент, отсчитываемый по касательной к кривой множества решений системы нелинейных уравнений, доставляет задаче наилучшую обусловленность, т. е. является наилучшим. Более того, этот результат был обобщен на системы обыкновенных дифференциальных, дифферен-

18 Григолюк Э. И., Шалашилин В. И. Метод продолжения решения по параметру в нелинейных задачах механики твердого деформируемого тела. М. : Наука, 1988. 232 с.

19 Allgower E. L., Georg K Introduction to Numerical Continuation Methods. Berlin ; Heidelberg : Springer-Verlag, 1990. 388 p.

20 Шалашилин В. И., Кузнецов Е. Б. Метод продолжения решения по параметру и наилучшая параметризация в прикладной математике и механике. М. : Эдиториал УРСС, 1999. 224 с.

циально-алгебраических и функционально-дифференциальных уравнений. В дальнейших работах Е. Б. Кузнецова было рассмотрено обобщение полученных результатов на многомерный случай21, а также дан алгоритм прохождения точек бифуркации различной коразмерности с использованием редукции Ляпунова-Шмидта22. Один из последних результатов связан с развитием нового подхода - метода продолжения решения по модифицированному наилучшему аргументу, отсчитываемому в направлении, близком к касательному, для систем обыкновенных дифференциальных уравнений [32-33]. Кроме применения продолжения решения по наилучшему параметру к решению жестких задач в статьях Н. Н. Калитки-на и соавторов, данный подход успешно использовался А. А. Семеновым при исследовании устойчивости панелей конических оболочек из ортотропных материалов при деформировании [34], а также иностранными исследователями - при решении гиперболических систем с предельными особыми точками и физически нелинейных задач [35-36], и т. д. Специальные виды параметров продолжения используют в своих работах Е. А. Лопаницын при расчете тонких пологих оболочек с учетом конечных прогибов23 и С. С. Гаврюшин -при расчете напряжений и деформаций сложных стержневых и оболочечных элементов конструкций24.

Материалы и методы

Приведем основные численные методы решения задач с контрастными структурами, используемые в статье. Известно, что любое обыкновен-

ное дифференциальное уравнение п-го порядка можно свести к системе п обыкновенных дифференциальных уравнений первого порядка. Поэтому сформулируем методы решения для системы дифференциальных уравнений вида

^ = f(У»У2,•••,Уп) 1 =1 п (3)

с начальными условиями

Уi(0) = Уi,o, i =1 П (4)

где fI (t, у,, у2,..., уп) - заданные функции.

Явный метод Эйлера

Одним из наиболее простых методов численного решения задачи Коши является явный метод Эйлера, схема которого на ^ом шаге имеет вид25:

У и к+1 = У ¡, к + Ь • Л (к, Уи, У2,к, • • •, Уп,к),

= *к + Ь, 1 = 1, •••, п, (5)

где h - шаг интегрирования по аргументу

Таким образом, решение задачи (3)-(4) на каждом шаге сводится к вычислению рекуррентных соотношений (5).

Явный метод Рунге - Кутты четвертого порядка точности

Несмотря на свою простоту, при необходимости высокой точности решения явный метод Эйлера даже для простых нежестких задач может требовать значительного уменьшения шага интегрирования, что приводит к уве-

21 Кузнецов Е. Б. Некоторые приложения метода продолжения решения по наилучшему параметру. М. : Изд-во МАИ, 2013. 160 с.

22 Кузнецов Е. Б. Параметризация краевых задач и прохождение точек бифуркации. М. : Изд-во МАИ, 2016. 160 с.

23 Григолюк Э. И., Лопаницын Е. А. Конечные прогибы, устойчивость и закритическое поведение тонких пологих оболочек. М. : Изд-во МАМИ, 2004. 162 с.

24 Гаврюшин С. С., Барышникова О. О., Борискин О. Ф. Численный анализ элементов конструкций машин и приборов. М. : Изд-во МГТУ им. Н. Э. Баумана, 2014. 479 с.

25 Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. М. : Мир, 1990. 512 с.

личению вычислительных затрат. Для устранения этого недостатка можно использовать методы более высокого порядка точности. Наиболее часто применяется метод Рунге - Кутты четвертого порядка точности, схема которого записывается в виде26:

yi,к+1 yi,к + г ' 6

( + 2 К 2„к + 2 К 3а + К 4, л),

сленного с шагом к, можно записать

где

4+1 = к+А

в виде

А, =р\ + О(к^),

где выражение рк - главная часть локальной погрешности, вычисляемая по формуле

V211 - Vй

Р _ У к У к 2 2р-1 '

(6)

- ква-

К11Л = ^ ((> У1,к Уп,к ) К 2,к = / к +1 уи + '.., упЛ + 2 К!пЛ

К 3,-,к = / {*к + 2' У1,к + 2 •К 21к'•••' Уп,к + 2-К 2п,к у К 4,к = (к + Ь,Уи + к •К 31,к упЛ + А •К Зп,к).

Переменный шаг интегрирования и метод Рунге-Ромберга-Ричардсона

На практике показано, что постоянный шаг интегрирования целесообразно использовать только для нежестких задач. В случае, когда задача (3)-(4) является плохо обусловленной или жесткой, более эффективным становится применение переменного шага интегрирования. Существует множество методов построения неравномерных сеток, некоторые из которых приведены в работах Н. Н. Калиткина и соавт.27, О. Б. Арушаняна и С. Ф. Залеткина28. Традиционно для смены шага используется правило Рунге, известное также как метод Рунге-Ромберга-Ричардсона.

Согласно правилу Рунге, локальную погрешность численного решения задачи (3)-(4) на ^ом шаге, вычи-

точность Рк

Здесь ук = ('• ••.у„,к)Т, || „_

дратичная (евклидова) норма вектора, р - порядок точности численного метода, верхними индексами к и 2к обозначены используемые шаги интегрирования.

Используя значение (6) как оценку локальной погрешности, можно реализовать процедуру смены шага. Задавая численного решения с1, при следует уменьшить вдвое. Процедура уменьшения шага повторяется до момента, когда рк < е1. Если значение рк <е1/ 2р, то шаг увеличивается вдвое.

Неявный метод Эйлера

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

26 Там же.

27 Вычисления на квазиравномерных сетках : монография / Н. Н. Калиткин [и др.]. М. : Физмат-лит, 2005. 224 с.

28 Арушанян О. Б., Залеткин С. Ф. Численное решение обыкновенных дифференциальных уравнений на Фортране. М. : Изд-во МГУ, 1990. 336 с.

не только является одним из наиболее простых методов данного класса, но и обладает значительной устойчиво-стью30. Его схема имеет вид31:

yi,k+1 = yi,k + ^' ^ (к+1, У1,к+1, У2,к+1,' •',Уп,к+1),

К+\ = К + Ь, ' =1,П (7)

Решение задачи (3)-(4) находится на каждом шаге из системы нелинейных уравнений (7). Система уравнений в работе решается двумя методами: простых итераций и Ньютона32.

Структура системы (7) делает удобным в использовании метод простых итераций, который сводится на каждом шаге к итерационному процессу

4+1 = У1, к + ^ ' (к+1, у1Д+1, • • •, УП',\+1 ) =

У7+1) = Ук)

„( i )

- J -

| , ф |

-- Pi ((

(j) k+1''

y( j ) ••'У n, k+1

Здесь у = 0,..., jk; ]к - количество итераций на к-ом шаге.

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

Записывая систему (7) в виде

где ук' - решение задачи (3)-(4) на к-ом шаге и ]-ой итерации; J - матрица Якоби системы (8).

Замечание 1

Важной задачей при использовании итерационных методов решения нелинейных уравнений и систем является выбор начального приближения. Для задачи (3)-(4) начальное приближение приходится выбирать на каждом шаге численного метода. Если на к-ом шаге получено решение у((к) и оно не обращает уравнения системы в бесконечность, то можно рекомендовать выбор начального приближения для к + 1-го шага в форме

^ = ЯА).

В противном случае начальное приближение можно выбирать в виде

yn, k+1!

= У i, k+1 — У i, k —

>1+1 = ykJk)

где S - вектор, компоненты которого выбираются из условий сходимости метода.

Замечание 2

Для задачи (3)-(4) сходимость итерационного процесса зависит от величины шага интегрирования h. Условия сходимости методов простых итераций и Ньютона накладывают ограничения на выбор h. При этом использование правила Рунге в большинстве случаев позволяет получить значения шага интегрирования, удовлетворяющие условиям сходимости. Лишь в окрестностях предельных особых точек и внутренних слоев необходимо дополнительно уменьшать шаг интегрирования. Наилучшая параметризация Неявные методы позволяют найти решения жестких и плохо обусловлен-

29 Там же.

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

31 Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие

и дифференциально-алгебраические задачи. М. : Мир, 1999. 685 с.

-h ■ fi (

, У1„

, — Уп

0 = 0 (8)

и вводя обозначение

Ф = (, Ф 2,..., Ф й),

запишем итерационный процесс по методу Ньютона в виде

ных задач Коши, но не могут сравниться по быстродействию с явными методами, т. к. на каждом шаге неявного метода приходится решать систему нелинейных уравнений. Этот недостаток удается устранить при помощи метода продолжения решения33, который заключается в замене исходного аргумента задачи Коши на новый. Можно выбирать различные аргументы продолжения решения, но наиболее часто используют наилучший аргумент. Данный аргумент отсчитывается по касательной к интегральной кривой рассматриваемой задачи и обладает рядом исключительных свойств. Для задачи (3)-(4) наилучший аргумент X записывается в скалярной форме:

dk2 = dy1 + --- + dy2n + dt2. (9)

Будем полагать, что все переменные у1,у2,...,уп и аргумент t зависят от X. Дополняя систему (3) соотношением (9) и разрешая полученную систему относительно производных по X, получим преобразованную систему

(У _ fi (t, У»:., УП )

dk Q(t,y^ :.,Уп) dt _ 1 dk Q( y^..^ Уп)

i _ 1,..., п,

(10)

где

_Q О, Уl,^^^, Уп )=_

= >/1 + ((,Ух,-;Уп) + - + /П(,У!,-,Уп).

Начальные условия (4) для системы (10) перепишутся в виде

Задача (10)—(11) обладает рядом важных для численного решения свойств34:

1. Квадратичная норма правой части системы (10) равна единице, т. е. устраняются все вычислительные трудности, связанные с неограниченным возрастанием правых частей системы (3).

2. Обусловленность системы (10) является наилучшей.

3. Показатель жесткости преобразованной системы (10) меньше, чем у исходной.

Все отмеченные свойства дают возможность решать преобразованную задачу (10)—(11) любыми численными методами, в том числе и явными.

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

Результаты исследования

В статье А. А. Белова и Н. Н. Калит-кина [3] для проверки эффективности численных методов при решении задач с контрастными структурами предложен ряд тестовых задач. Проанализируем решение одной из этих задач, называемой степенным тестом.

Постановка задачи

Рассмотрим начальную задачу для уравнения следующего вида [3]:

du dt

4(t) (u2 - a2 )2

(u2 + a2)

с начальным условием

ы(0) = 0.

(12)

(13)

При % (t) = %0 • cost решение задачи (12)—(13) имеет вид:

u (t) = ■

2E(t )a2

1 + 71 + 4a2Z2(t)'

(14)

y(0) = y,0, t(0) = 0, i = 1,...,n. (11) где S(t) = £0smt.

32 Калиткин Н. Н. Численные методы. СПб. : БХВ-Петербург, 2011. 592 с.

33 Шалашилин В. И., Кузнецов Е. Б. Метод продолжения решения по параметру и наилучшая

параметризация в прикладной математике и механике. М. 496

Эдиториал УРСС, 1999. 224 с.

Физико-математические науки

В дальнейших расчетах а = п. Жесткость задачи характеризуется множителем при периодической функции. Условно по значению параметра можно разбить задачи вида (12)—(13) на следующие классы [5]: при малых значениях (порядка нескольких единиц) задача является нежесткой, при £0 > 10 -жесткой и при £0 > 1000 - сверхжесткой. Это определяется видом правой части: при приближении значения косинуса к нулю скорость изменения решения, даже при большом параметре |0, стремится к нулю, а при стремлении значения косинуса к единице значение правой части становится наибольшим. Таким образом, происходит смена плавно меняющихся компонент решения с быстро меняющимися переходными зонами, которые тем короче, чем больше значение |0. Данные особенности видны на рис. 1.

Явные методы решения

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

даны выше, представлены на рис. 1 и в табл. 1.

На рис. 1 изображены кривые, соответствующие решениям задачи (12)—(13) для значений = 1; 10; 100, полученные методом Рунге-Кутты четвертого порядка точности с постоянным шагом h = 0,001. Полученные численные решения визуально совпадают с аналитическими. Для явного метода Эйлера кривые аналогичны изображенным на рис. 1.

В табл. 1 приведены данные о средней погрешности с, вычисленной с использованием аналитического решения (14), и времени счета (с решений задачи (12)—(13), полученных явными методами с постоянным шагом интегрирования. Из нее видно, что метод Рунге-Кутты четвертого порядка выигрывает у явного метода Эйлера в точности, хотя и требует большего времени счета. Стоит также отметить, что при повышении значения параметра ¿0 невозможно получить решение при малых шагах интегрирования. Так, для ¿¿0 = 10 невозможно получить решение при

Р и с. 1. Численное решение задачи (12)-(13) при £0 = 1; 10; 100, метод Рунге-Кутты четвертого порядка с постоянным шагом h = 10-3

F i g. 1. Numeric solution of problem (12)-(13) for £0 = 1; 10; 100, fourth order Runge-Kutta method with constant step size h = 10-3

Т а б л и ц а 1 T a b l e 1

Погрешность решения и время счета для задачи (12)-(13), явные методы с постоянным шагом

Computational error and computational time for problem (12)-(13), explicit methods with constant step size

Параметры / Parameters Явный метод Эйлера / Explicit Euler method Метод Рунге-Кутты четвертого порядка / Fourth order Runge-Kutta method

¿0 h £ tc, c £ tc, c

1 0,1 0,4548 0,009 2,1332-10-4 0,0081

0,01 0,0361 0,0095 2,2776-10-8 0,0656

0,001 0,0035 0,0255 2,146110-12 0,1732

10 0,1 - - - -

0,01 0,0826 0,0109 2,3893-10-5 0,0168

0,001 0,008 0,0301 2,5136-10-9 0,0759

100 0,1 - - - -

0,01 - - - -

0,001 0,0127 0,0381 2,4843-10-6 0,1011

h = 0,1, а при ¿0 = 1 000 нельзя получить решение при h > 10-3.

Здесь и далее расчет проведен на персональном компьютере с процессором DualCore Intel Core i5-660 3477 MHz, оперативной памятью 2,00 ГБ DDR3 SDRAM Kingston «ValueRAM» KVR1333D3N9/2G (2 шт.), видеокартой ATI Radeon HD 5770 (Juniper) 1,00 ГБ и 64-разрядной операционной системой Windows 10. Используемые методы реализованы в вычислительной среде Matlab R2012b.

Замечание 3

В рассматриваемом диапазоне шагов интегрирования не удалось получить решение задачи (12)-(13) для ¿0 > 1 000.

Анализируя полученные результаты, отметим, что для решения задачи с фиксированным ¿0 шаг интегрирования h не может превышать величину 1/^0. Это говорит о малоэффективности применения явных методов для решения сверхжестких задач при больших значениях ¿0.

Неявный метод Эйлера

Явные методы с постоянным шагом интегрирования позволяют получить ре-

шение задачи (12)—(13) только при умеренных значениях £0. Уже при £0 = 103 в рассматриваемом диапазоне шагов интегрирования не удается найти решение. Можно ожидать, что неявные методы за счет расширения области устойчивости позволят получить решение при больших шагах интегрирования.

Используем неявный метод Эйлера. На ^ом шаге решение задачи (12)—(13) неявным методом Эйлера находится как решение нелинейного уравнения вида:

cos(tk+1) (+1 )2 - «2 )

(uk+1 )2+ а1

(15)

Уравнение (15) решается методом простых итераций и методом Ньютона. Итерационный процесс метода простых итераций задается соотношением:

U£1) = u - h ■

^cos^+i) (( )2 - a2 )

( j )

+ a

а для метода Ньютона -

(j+i) _ и(j) - F(tk+i, uk+1)

_ F '(t u (j))

Ги V k+1>" k+1)

(16)

В формуле (16)

F(t, u) = u - uk + h

<%0cost (u2 - a a)

и + a

2L,cost (u2 - a2 )u (u2 + 3a2)

F(t,u) = 1-h• \ 2 -К

(u + a )

где FU(t, u) =

dF (t, u) du

Кривые, полученные при решении задачи (12)—(13) неявным методом Эйлера, имеют вид, аналогичный представленному на рис. 1. В табл. 2 даны значения средней погрешности с и времени счета tc решений задачи (12)—(13), полученные неявным методом Эйлера с постоянным шагом интегрирования

при использовании метода простых итераций и Ньютона.

Расчеты показывают, что по сравнению с явным методом Эйлера неявный метод Эйлера позволяет получить результаты с меньшей погрешностью. Однако время счета в связи с усложнением вычислительного процесса как для метода простых итераций, так и метода Ньютона, может увеличиваться на порядок и более. При малой жесткости метод Ньютона позволяет быстрее получить решение задачи по сравнению с методом простых итераций, но с повышением жесткости уступает ему как в точности, так и в быстродействии. Это можно связать и с особенностями реализации метода Ньютона (в частности, с выбором начального приближения на каждом шаге) и с затруднением при переходе через внутренний слой для неявного метода Эйлера. Таким образом, используемый неявный метод не позволяет значительно улучшить результаты явных методов, затрачивая

Т а б л и ц а 2 T a b l e 2

Погрешность решения и время счета для задачи (12)-(13), неявный метод Эйлера с постоянным шагом

Computational error and computational time for problem (12)-(13), implicit Euler method with constant step size

Параметры / Parameters Неявный метод Эйлера / Implicit Euler method

Метод простой итерации / Fixed-point iteration method Метод Ньютона / The Newton's method

¿0 h £ tc, c £ tc, c

1 0,1 0,3137 0,0662 0,3138 0,0255

0,01 0,034 0,1774 0,0345 0,1457

0,001 8,6629T0-4 1,188 0,0035 0,6134

10 0,1 - - - -

0,01 0,074 0,1864 0,0787 0,0878

0,001 0,0029 1,1653 0,008 0,7302

100 0,1 - - - -

0,01 - - - -

0,001 0,0078 0,9977 0,0126 1,5954

при этом большее время счета. Для повышения точности можно, как и в статье А. А. Белова и Н. Н. Калиткина «Численные методы решения задач Коши с контрастными структурами» [5], рекомендовать использование более сложных методов, например, неявных обратных схем Рунге-Кутты [37] или явно-неявных комплексных схем Розенброка [38].

Наилучшая параметризация

Из приведенных выше расчетных данных видно, что ни явные методы, ни неявный метод Эйлера не позволяют получить решение задачи (12)—(13) без значительных вычислительных затруднений. Применим к этой задаче метод продолжения решения по наилучшему аргументу35, описанный выше.

Наилучший аргумент X для рассматриваемой задачи запишется в скалярном виде:

dX2 = du2 + dt2.

(17)

Преобразуя к аргументу (17) систему (12), получим:

du

-%(?) (и2 - а2)

¿X и2

dt

¿Х

+ а2 )2 ) (и2 - а2 )4 (и2 + а2)

(18)

и 2 - а2)

Начальные условия (13) для системы (18) примут вид:

и(0) = 0, I(0) = 0. (19)

Как отмечалось выше, преобразованная задача (18)—(19) обладает рядом преимуществ при численном решении. В табл. 3 приведены результаты расчета преобразованной задачи (18)—(19) методом Рунге-Кутты четвертого по-

рядка точности с постоянным шагом (качественно для метода Эйлера результаты аналогичны), где I - шаг по аргументу X. Кривые, полученные при численном решении, аналогичны приведенным на рис. 1. Для сравнения в табл. 3 приведены результаты решения исходной задачи (12)—(13) с шагом, подбираемым из условия равенства порядков средней ошибки решения исходной и преобразованной задач.

Видно, что для непреобразованной задачи (12)—(13) для достижения того же порядка ошибки приходится брать шаг до нескольких порядков меньше по сравнению с преобразованной. При этом, несмотря на то, что интервал изменения значений аргумента X больше, чем у аргумента I, для задачи (18)—(19) удается сократить время счета от 0,4 (при = 1) до 10 раз (при 4 = 1 000).

В табл. 4 приведены результаты расчета задач (12)—(13) и (18)—(19) методом Рунге-Кутты четвертого порядка с переменным шагом. Процедура изменения шага описана выше. Для преобразованной задачи (18)—(19) параметр точности для правила Рунге 9 = 10-1 . Для исходной задачи (12)—(13) параметр 0 выбирался из условия равенства порядков средней погрешности для обеих задач.

Результаты показывают, что и при переменном шаге интегрирования использование наилучшего аргумента X позволяет получить вычислительные преимущества: удается сократить время счета минимум в два раза. Для сверхжестких задач уменьшение времени счета может достигать нескольких порядков. Отметим также, что, в отличие от преобразованной задачи, решение исходной задачи зависит от начального шага; в особенности это характерно для больших значений £0.

Однако стоит указать и на ряд недостатков наилучшей параметризации. Преобразованная задача (18)—(19)

34 Там же.

Т а б л и ц а 3 T a b l e 3

Средняя ошибка и время счета задач (12)-(13) и (18)-(19), метод Рунге-Кутты четвертого порядка с постоянным шагом Computational error and computational time for problems (12)-(13) and (18)-(19), fourth order Runge-Kutta method with constant step size

¿0 Параметризованная задача / Parameterized problem Непараметризованная задача / Nonparametrized problem

l £ tc, c h e tc, c

1 0,1 5,4369T0-7 0,0213 0,02 3,6307^10-7 0,0418

0,01 2,9799T0-11 0,183 0,002 3,6427^10-11 0,2405

0,001 5,0522T0-14 0,5505 0,0002 6,4592T0-13 2,0987

10 0,1 2,4647T0-4 0,0236 0,02 3,0375^10-4 0,0422

0,01 3,2723T0-9 0,0969 0,001 2,5136^10-9 0,4414

0,001 3,7533T0-12 0,5886 2T0-4 3,662T0-12 2,2979

100 0,1 0,3525 0,0289 0,004 0,0012 0,1303

0,01 4,834-10-6 0,1178 0,001 2,4843-10-6 0,4484

0,001 5,0063T0-11 0,5289 5T0-5 1,707fr10-11 8,2501

1000 0,1 - - - - -

0,01 0,0047 0,2086 4,7T0-4 0,0012 0,9134

0,001 4,0885T0-8 0,5288 7T0-5 5,9638^10-8 5,9466

Т а б л и ц а 4 T a b l e 4

Средняя ошибка и время счета задач (12)-(13) и (18)-(19), метод Рунге-Кутты четвертого порядка с переменным шагом Computational error and computational time for problems (12)-(13) and (18)-(19), fourth order Runge-Kutta method with variable step size

Параметры / Parameters Параметризованная задача / Parameterized problem Непараметризованная задача / Nonparametrized problem

he, / I0 в £ tc, c в e tc, c

1 0,1 10-12 1,1561T0-10 0,2388 10-10 2,1774T0-10 0,4892

0,01 10-12 1,6079T0-10 0,2592 3-10-10 1,1338T0-10 0,571

0,001 10-12 1,7008T0-10 0,1889 10-10 2,1458T0-12 1,3038

10 0,1 10-12 2,365^10-9 0,2065 10-7 1,516410-9 1,6533

0,01 10-12 2,2218T0-9 0,2812 4T0-8 1,1642 10-9 2,025

0,001 10-12 1,6625 10-9 0,2691 10-7 2,5134^10-9 1,3089

100 0,1 10-12 3,497^10-8 0,2982 10-4 5,8346^10-8 3,2128

0,01 10-12 3,1933 10-8 0,3346 10-5 2,3953^10-8 4,0536

0,001 10-12 3,0792T0-8 0,3015 3T0-5 2,3106^10-8 4,7691

1000 0,1 10-12 4,2953T0-7 0,2729 - - -

0,01 10-12 3,2001T0-7 0,2318 10-3 9,1802-10-8 15,794

0,001 10-12 8,7845T0-7 0,1772 10-2 6,4254^10-7 9,8967

имеет более сложный вид и большую размерность по сравнению с исходной (12)—(13). Для рассматриваемой задачи это несущественно, но при решении задач большой размерности преимущества, даваемые наилучшей параметризацией, могут перекрываться сложностью преобразованной задачи.

Еще одним недостатком является заниженная оценка погрешности по правилу Рунге для преобразованной задачи. Уже при ¿¿0 = 1 оценка погрешности на два порядка ниже, чем реальная средняя погрешность. Для сверхжестких задач разница возрастает до пяти порядков. Диаметрально противоположная ситуация характерна для исходной задачи. Здесь правило Рунге дает завышенную оценку погрешности начиная с ¿0 = 10. Если посмотреть на характер изменения шага интегрирования, изображенный на рис. 2 в полулогарифмических координатах, то видно, что уже при значении ¿0 = 30 значение шага интегрирования не поднимается выше 10—5, снижаясь до 10—12 в окрестности пограничных и внутреннего слоев. При этом значение шага изме-

няется плавно, т. е. после прохождения внутреннего слоя шаг начинает плавно возрастать до максимального значения, а при подходе к внутреннему слою -плавно убывать до минимума (при этом ширина переходного слоя уменьшается при увеличении ¿0 ). Данный механизм изменения шага интегрирования и приводит к уменьшению погрешности решения на горизонтальных участках до значений меньших, чем задаваемая точность 0.

Для преобразованной задачи кривая решения в зависимости от аргумента X изменяется более плавно, как это отмечено в статьях А. А. Белова и Н. Н. Ка-литкина [3; 5]. По этой причине такого же падения погрешности не происходит; наоборот, в узловых точках, соответствующих внутренним слоям, концентрируется максимальная погрешность. В заключение отметим, что, исходя из характера изменения шага интегрирования, можно объяснить и увеличение времени счета для исходной задачи. Помимо усложнения процесса решения, связанного со сменой шага, уже при относительно малом значении ¿0 = 30 шаг

ЦОД)

Р и с. 2. Шаг интегрирования в полулогарифмических координатах для задачи (12)-(13) при 4 = 30, метод Рунге-Кутты четвертого порядка, в = 10-4

F i g. 2. Step size in semilogarithmic coordinates for problem (12)-(13), ¿0 = 30, fourth order Runge-Kutta method, в = 10-4

интегрирования не поднимается выше 10—5, при повышении же значения ¿0 шаг интегрирования будет уменьшаться еще сильнее, а вместе с этим будет возрастать время счета.

Замечание 4

В отличие от непреобразованной задачи (12)—(13), переход к наилучшему аргументу позволяет получить хорошие результаты при параметре жесткости ¿0 большем, чем рассмотренные в табл. 3 и табл. 4. При задаваемой точности 0 = 10—12 и начальном шаге h = 0,001 удается построить решение вплоть до ¿0 = 106.

Замечание 5

Использование неявного метода Эйлера с переменным шагом не позволило получить результаты значительно лучшие, чем приведенные в табл. 2. Получить решение со средней погрешностью £ порядка 10—2 удается, но время счета по сравнению с результатами табл. 2 возрастает до 10 раз при ¿0 < 100. При большей точности (£ < 10—3) или большей жесткости (¿0 > 1000) результаты либо не удавалось получить, либо время счета составляло до 1 000 с.

Обсуждение и заключение

В статье рассмотрено применение традиционных явных и неявных методов решения задачи Коши с контрастными структурами. Получены численные решения задачи (12)—(13) с постоянным и переменным шагом, изменяемым по правилу Рунге. Анализируя результаты, приведенные на рис. 1 и в табл. 1-2, отметим следующее.

1. Все полученные расчетные данные хорошо согласуются как с аналитическим решением, так и с результатами, полученными А. А. Беловым и Н. Н. Калиткиным [3; 5].

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

h не может превосходить значения 1/^0. Уже при ¿0 > 1 000 решение в рассматриваемом диапазоне шагов интегрирования построить не удалось.

3. Используемый в статье неявный метод Эйлера дает возможность уменьшить погрешность полученного решения для малых ¿0, но не позволяет построить решение при ¿0 > 1 000. При этом, в силу особенностей реализации, время счета для неявного метода Эйлера значительно превышает аналогичные значения для явных методов.

4. Метод Ньютона, используемый при решении нелинейных уравнений в неявном методе Эйлера, при малых ¿0 позволяет сократить время счета по сравнению с методом простых итераций, но уступает ему в точности.

5. Неявный метод Эйлера вызывает значительные затруднения при переходе через внутренний слой, что сказывается на быстродействии и точности.

Для устранения отмеченных недостатков исходная задача (12)—(13) преобразуется к наилучшему аргументу X. Преобразованная задача (18)—(19) обладает рядом вычислительных преимуществ (см. табл. 3-4).

1. При решении задачи (18)—(19) с постоянным шагом удается получить большую точность. При этом для исходной задачи, в зависимости от значения ¿0, аналогичные результаты можно получить с меньшим (от нескольких раз до нескольких порядков) шагом. Время счета для преобразованной задачи меньше от 0,4 (при ¿0 = 1) до 10 раз (при ¿0 = 1 000).

2. При использовании переменного шага интегрирования задача (18)—(19) также выигрывает в быстродействии. Переход к аргументу X позволяет при одинаковых значениях средней погрешности сократить время счета от двух раз для малых значений ¿0 до нескольких порядков при ¿0 = 1 000. Более того, в отличие от исходной, преобразованную задачу можно решать вплоть до значений ¿0 = 106.

3. Для преобразованной задачи (18)—(19) процедура смены шага слабо зависит от начального шага /0, в то время как для исходной задачи решение существенно зависит от выбора начального шага

Помимо значительных преимуществ задача (18)—(19) имеет и ряд недостатков.

1. Правило Рунге для преобразованной задачи дает заниженную оценку локальной погрешности. Уже при ¿0 = 1 оценка погрешности на два порядка ниже, чем реальная средняя погрешность. Для сверхжестких задач разница возрастает до пяти порядков. Диаметрально противоположная ситуация отмечается с исходной задачей. Для нее правило Рунге дает завышенную оценку погрешности начиная с £ = 10. Это необходимо учитывать при построении решения.

2. В отличие от исходного уравнения (12), преобразованная система (18) имеет большую размерность и усложненный вид. В ряде случаев это может

перекрывать преимущества наилучшей параметризации.

Несмотря на отмеченные недостатки, наилучшая параметризация является одним из наиболее эффективных методов решения плохо обусловленных задач, к которым можно отнести и задачи с контрастными структурами. Показано, что традиционные явные и неявные методы имеют существенные ограничения в использовании при больших ¿0 и малоэффективны. В то же время переход к наилучшему аргументу позволяет получить решение рассматриваемой задачи даже при больших ¿0 быстро и с приемлемой точностью. В дальнейшем для подтверждения полученных результатов необходимо рассмотреть решение еще ряда тестовых и прикладных задач данного класса. Однако, опираясь на уже полученные результаты решения различных задач36—38 [33], можно предположить, что метод продолжения решения по наилучшему аргументу будет также эффективен и для других задач данного класса.

СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ

1. Тихонов А. Н. О зависимости решений дифференциальных уравнений от малого параметра // Матемематический сборник. 1948. Т. 22 (64), № 2. С. 193-204. URL: http://mi.mathnet.ru/msb6075

2. Васильева А. Б., Бутузов В. Ф., Нефедов Н. Н. Контрастные структуры в сингулярно возмущенных задачах // Фундаментальная и прикладная математика. 1998. Т. 4, № 3. С. 799-851. URL: http://mi.mathnet.ru/fpm344

3. Belov A. A., Kalitkin N. N. Features of calculating contrast structures in the Cauchy problem // Mathematical Models and Computer Simulations. 2017. Vol. 9, Issue 3. P. 281-291. DOI: https://doi. org/10.1134/S2070048217030048

4. Belov A. A., Kalitkin N. N. Curvature-based grid step selection for stiff Cauchy problems // Mathematical Models and Computer Simulations. 2017. Vol. 9, Issue 3. P. 305-317. DOI: https://doi. org/10.1134/S207004821703005X

5. Белов А. А., Калиткин Н. Н. Численные методы решения задач Коши с контрастными структурами // Моделирование и анализ информационных систем. 2016. Т. 23, № 5. С. 529-538. DOI: https://doi.org/10.18255/1818-1015-2016-5-529-538

35 Там же.

36 Там же.

37 Кузнецов Е. Б. Некоторые приложения метода продолжения решения по наилучшему параметру. М. : Изд-во МАИ, 2013. 160 с.

38 Кузнецов Е. Б. Параметризация краевых задач и прохождение точек бифуркации. М. : Изд-во МАИ, 2016. 160 с.

6. Тихонов А. Н. О системах дифференциальных уравнений, содержащих параметры // Матемематический сборник. 1950. Т. 27 (69), № 1. С. 147-156. URL: http://mi.mathnet.ru/msb5907

7. Тихонов А. Н. Системы дифференциальных уравнений, содержащие малые параметры при производных // Матемематический сборник. 1952. Т. 31 (73), № 3. С. 575-586. URL: http:// mi.mathnet.ru/msb5548

8. Butuzov V. F., Vasileva A. B., Nefedov N. N. Asymptotic theory of contrast structures (review) // Automatics and Remote Control. 1997. Vol. 58, Issue 7. P. 1068-1091. Available at: http://mi.mathnet.ru/ at2615

9. Butuzov V. F., Levashova N. T., Melnikova A. A. A steplike contrast structure in a singularly perturbed system of elliptic equations // Computational Mathematics and Mathematical Physics Journal. 2013; 53(9):1239-1259. DOI: https://doi.org/10.1134/S0965542513090054

10. Бутузов В. Ф., Денисов И. В. Угловой пограничный слой в нелинейных эллиптических задачах, содержащих производные первого порядка // Моделирование и анализ информационных систем. 2014. Т. 21, № 1. С. 7-31. URL: http://mi.mathnet.ru/mais356

11. Бутузов В. Ф., Белошапко В. А. Сингулярно возмущенная эллиптическая задача Дирихле с кратным корнем вырожденного уравнения // Моделирование и анализ информационных систем. 2016. Т. 23, № 5. С. 515-528. DOI: https://doi.org/10.18255/1818-1015-2016-5-515-528

12. Butuzov V. F., Bychkov A. I. Asymptotics of the solution of the initial boundary value problem for a singularly perturbed parabolic equation in the case of a triple root of the degenerate equation // Computational Mathematics and Mathematical Physics Journal. 2016. Vol. 56, Issue 4. P. 593-611. DOI: https://doi.org/10.1134/S0965542516040060

13. Бутузов В. Ф. О контрастных структурах с многозонным внутренним слоем // Моделирование и анализ информационных систем. 2017. Т. 24, № 3. С. 288-308. DOI: https://doi.org/10.18255/1818-1015-2017-3-288-308

14. Козлов М. В., Щенников В. Н. Асимптотическая устойчивость однородных сингулярных систем // Вестник Мордовского университета. 2017. Т. 27, № 4. С. 546-554. DOI: https://doi.org/10.15507/0236-2910.027.201704.546-554

15. Нефедов Н. Н., Никулин Е. И. Существование и асимптотическая устойчивость периодического решения с внутренним переходным слоем в задаче со слабой линейной адвекцией // Моделирование и анализ информационных систем. 2018. Т. 25, № 1. С. 125-132. DOI: https:// doi.org/10.18255/1818-1015-2018-1-125-132

16. Антипов Е. А., Левашова Н. Т., Нефедов Н. Н. Асимптотическое приближение решения уравнения реакция-диффузия-адвекция с нелинейным адвективным слагаемым // Моделирование и анализ информационных систем. 2018. Т. 25, № 1. С. 18-32. DOI: https://doi.org/10.18255/1818-1015-2018-1-18-32

17. Давыдова М. А., Нефедов Н. Н. Существование и устойчивость контрастных структур в многомерных задачах реакция-диффузия-адвекция в случае сбалансированной нелинейности // Моделирование и анализ информационных систем. 2017. Т. 24, № 1. С. 31-38. DOI: https:// doi.org/10.18255/1818-1015-2017-1-31-38

18. Antipov E. A., Levashova N. T., Nefedov N. N. Asymptotics of the front motion in the reac-tion-diffusion-advection problem // Computational Mathematics and Mathematical Physics Journal. 2014. Vol. 54, Issue 10. P. 1536-1549. DOI: https://doi.org/10.1134/S0965542514100029

19. Нефедов Н. Н. Общая схема асимптотического исследования устойчивых контрастных структур // Нелинейная динамика. 2010. Т. 6, № 1. С. 181-186. URL: http://mi.mathnet.ru/nd65

20. Efstathiou C., Luhar M. Mean turbulence statistics in boundary layers over high-porosity foams // Journal of Fluid Mechanics. 2018. Vol. 841. P. 351-379. DOI: https://doi.org/10.1017/jfm.2018.57

21. Comparison of turbulent boundary layers over smooth and rough surfaces up to high Reynolds numbers / D. T. Squire [et al.] // Journal of Fluid Mechanics. 2016. Vol. 795. P. 210-240. DOI: https:// doi.org/10.1017/jfm.2016.196

22. Swaters G. E. Internal dissipative boundary layers in the cross-equatorial flow of a grounded deep western boundary current // Geophysical & Astrophysical Fluid Dynamics. 2017. Vol. 111, no. 2. P. 91-114. DOI: https://doi.org/10.1080/03091929.2017.1287909

23. Kumar D. A parameter-uniform method for singularly perturbed turning point problems exhibiting interior or twin boundary layers // International Journal of Computer Mathematics. 2018. P. 1-18. DOI: https://doi.org/10.1080/00207160.2018.1458098

24. Xu H., Jin Y. L. The contrast structures for a class of singularly perturbed systems with het-eroclinic orbits // Discrete Dynamics in Nature and Society. 2016. Article ID 6405853. DOI: http:// dx.doi.org/10.1155/2016/6405853

25. Belov A. A., Kalitkin N. N., Kuzmina L. V. Modeling of chemical kinetics in gases // Mathematical Models and Computer Simulations. 2017. Vol. 9, Issue 1. P. 24-39. DOI: https://doi.org/10.1134/ S2070048217010057

26. Lahaye M. E. Une metode de resolution d'une categorie d'equations transcendentes // Comptes rendus hebdomadaires des seances de L'Academie des sciences. 1934. Vol. 198, no. 21. P. 1840-1842.

27. Lahaye M. E. Solution of system of transcendental equations // Académie royale de Belgique. Bulletin de la Classe des sciences. 1948. Vol. 5. P. 805-822.

28. Давиденко Д. Ф. Об одном новом методе численного решения систем нелинейных уравнений // Доклады Академии наук СССР. 1953. Т. 88, № 4. С. 601-602.

29. Давиденко Д. Ф. О приближенном решении систем нелинейных уравнений // Украинский математический журнал. 1953. Т. 5, № 2. С. 196-206.

30. Ворович И. И., Зипалова В. Ф. К решению нелинейных краевых задач теории упругости методом перехода к задаче Коши // Прикладная математика и механика. 1965. Т. 29, вып. 5. С. 894-901.

31. Riks E. The application of Newton's method to the problem of elastic stability // Journal of Applied Mechanics. 1972; Vol. 39, Issue 4. P. 1060-1065. URL: http://appliedmechanics.asmedigitalcollec-tion.asme.org/article.aspx?articleid=1400638

32. Kuznetsov E. B., Leonov S. S. Parametrization of the Cauchy problem for systems of ordinary differential equations with limiting singular points // Computational Mathematics and Mathematical Physics Journal. 2017. Vol. 57, Issue 6. P. 931-952. DOI: https://doi.org/10.1134/S0965542517060094

33. Кузнецов Е. Б., Леонов С. С. Примеры параметризации задачи Коши для систем обыкновенных дифференциальных уравнений с предельными особыми точками // Журнал вычислительной математики и математической физики. 2018. Т. 58, № 6. С. 914-933. DOI: https:// doi.org/10.7868/S0044466918060066

34. Semenov A. A. Strength and stability of geometrically nonlinear orthotropic shell structures // Thin-Walled Structures. 2016. Vol. 106. P. 428-436. DOI: https://doi.org/10.1016/j.tws.2016.05.018

35. May S., Vignollet J., de Borst R. A new arc-length control method based on the rates of the internal and the dissipated energy // Engineering Computations. 2016. Vol. 33, Issue 1. P. 100-115. DOI: https://doi.org/10.1108/EC-02-2015-0044

36. A local pseudo arc-length method for hyperbolic conservation laws / X.Wang [et al.] // Acta Me-chanica Sinica. 2015. Vol. 30, no. 6. P. 956-965. DOI: https://doi.org/10.1007/s10409-014-0091-0

37. Kalitkin N. N., Poshivaylo 1 P. Computations with inverse Runge-Kutta schemes // Mathematical Models and Computer Simulations. 2014. Vol. 6, Issue 3. P. 272-285. DOI: https://doi.org/10.1134/ S2070048214030077

38. Rosenbrock H. H. Some general implicit processes for the numerical solution of differential equations // The Computer Journal. 1963. Vol. 5, no. 4. P. 329-330. DOI: https://doi.org/10.1093/com-jnl/5.4.329

Поступила 05.07.2018; принята к публикации 03.09.2018; опубликована онлайн 28.12.2018

Об авторах:

Кузнецов Евгений Борисович, профессор кафедры моделирования динамических систем, ФГБОУ ВО «Московский авиационный институт (национальный исследовательский университет)» (125993, Россия, г. Москва, Волоколамское ш., д. 4), доктор физико-математических наук, профессор, ResearcherID: S-9576-2018, ORCID: https://orcid.org/0000-0002-9452-6577, [email protected]

Леонов Сергей Сергеевич, доцент кафедры моделирования динамических систем, ФГБОУ ВО «Московский авиационный институт (национальный исследовательский университет)» (125993, Россия, г. Москва, Волоколамское ш., д. 4), кандидат физико-математических наук, ResearcherID: Е-7641-2017, ORCID: https://orcid.org/0000-0001-6077-0435, [email protected]

Цапко Екатерина Дмитриевна, студент кафедры моделирования динамических систем, ФГБОУ ВО «Московский авиационный институт (национальный исследовательский университет)» (125993, Россия, г. Москва, Волоколамское ш., д. 4), ResearcherID: S-8613-2018, ORCID: https:// orcid.org/0000-0002-4215-3510, [email protected]

Заявленный вклад соавторов:

Е. Б. Кузнецов — проведение теоретических исследований по методу продолжения решения по параметру, формулирование основной концепции исследования, доработка текста; С. С. Леонов — обзор и анализ литературы, разработка вычислительных алгоритмов метода продолжения решения по наилучшему аргументу, формулировка выводов, редактирование текста; Е. Д. Цапко — разработка программного комплекса и проведение численных экспериментов, подготовка первоначального варианта рукописи, верстка текста.

Все авторы прочитали и одобрили окончательный вариант рукописи.

REFERENCES

1. Tikhonov A. N. On the dependence of the solutions of differential equations on a small parameter. Matematicheskiy sbornik = Mathematical Collection (New Series). 1948; 22(2):193-204. Available at: http://mi.mathnet.ru/msb6075 (In Russ.)

2. Vasileva A. B., Butuzov V. F., Nefedov N. N. Contrast structures in singularly perturbed problems. Fundamentalnaya iprikladnaya matematika = Fundamental and Applied Mathematics. 1998; 4(3):799-851. Available at: http://mi.mathnet.ru/fpm344 (In Russ.)

3. Belov A. A., Kalitkin N. N. Features of calculating contrast structures in the Cauchy problem. Mathematical Models and Computer Simulations. 2017; 9(3):281-291. DOI: https://doi.org/10.1134/ S2070048217030048

4. Belov A. A., Kalitkin N. N. Curvature-based grid step selection for stiff Cauchy problems. Mathematical Models and Computer Simulations. 2017; 9(3):305-317. DOI: https://doi.org/10.1134/S207004821703005X

5. Belov A. A., Kalitkin N. N. Numerical methods of solving Cauchy problems with contrast structures. Modelirovaniye i analiz informatsionnykh sistem = Modelling and Analysis of Information System. 2016; 23(5):529-538. DOI: https://doi.org/10.18255/1818-1015-2016-5-529-538 (In Russ.)

6. Tikhonov A. N. On systems of differential equations containing parameters. Matematicheskiy sbornik = Mathematical Collection. 1950; 27(1):147-156. Available at: http://mi.mathnet.ru/msb5907 (In Russ.)

7. Tikhonov A. N. Systems of differential equations containing a small parameter in front of the derivatives. Matematicheskiy sbornik = Mathematical Collection. 1952; 31(3):575-586. Available at: http:// mi.mathnet.ru/msb5548 (In Russ.)

8. Butuzov V. F., Vasileva A. B., Nefedov N. N. Asymptotic theory of contrast structures (review). Automatics and Remote Control. 1997; 58(7):1068-1091. Available at: http://mi.mathnet.ru/at2615

9. Butuzov V. F., Levashova N. T., Melnikova A. A. A steplike contrast structure in a singularly perturbed system of elliptic equations. Computational Mathematics and Mathematical Physics Journal. 2013; 53(9):1239-1259. DOI: https://doi.org/10.1134/S0965542513090054

10. Butuzov V. F., Denisov I. V. Corner boundary layer in nonlinear elliptic problems containing derivatives of first order. Modelirovaniye i analiz informatsionnykh sistem = Modeling and Analysis of Information Systems. 2014; 21(1):7-31. Available at: http://mi.mathnet.ru/mais356 (In Russ.)

11. Butuzov V. F., Beloshapko V. A. Singularly perturbed elliptic Dirichlet problem with a multiple root of the degenerate equation. Modelirovaniye i analiz informatsionnykh sistem = Modeling and Analysis of Information Systems. 2016; 23(5):515-528. DOI: https://doi.org/10.18255/1818-1015-2016-5-515-528 (In Russ.)

12. Butuzov V. F., Bychkov A. I. Asymptotics of the solution of the initial boundary value problem for a singularly perturbed parabolic equation in the case of a triple root of the degenerate equation. Computational Mathematics and Mathematical Physics Journal. 2016; 56(4):593-611. DOI: https://doi. org/10.1134/S0965542516040060

13. Butuzov V. F. On contrast structures with a multizonal interior layer. Modelirovaniye i analiz informatsionnykh sistem = Modeling and Analysis of Information Systems. 2017; 24(3):288-308. DOI: https://doi.org/10.18255/1818-1015-2017-3-288-308 (In Russ.)

14. Kozlov M. V., Shchennikov V. N. Asymptotic stability of homogeneous singular systems. Vestnik Mordovskogo universiteta = Mordovia University Bulletin. 2017; 27(4):546-554. DOI: https://doi.org/10.15507/0236-2910.027.201704.546-554 (In Russ.)

15. Nefedov N. N., Nikulin E. I. Existence and stability of the periodic solution with an interior transitional layer in the problem with a weak linear advection. Modelirovaniye i analiz informatsionnykh sistem = Modeling and Analysis of Information Systems. 2018; 25(1):125-132. DOI: https://doi.org/10.18255/1818-1015-2018-1-125-132 (In Russ.)

16. Antipov E. A., Levashova N. T., Nefedov N. N. Asymptotic approximation of the solution of the reaction-diffusion-advection equation with a nonlinear advective term. Modelirovaniye i analiz informatsionnykh sistem = Modeling and Analysis of Information Systems. 2018; 25(1):18-32. DOI: https:// doi.org/10.18255/1818-1015-2018-1-18-32 (In Russ.)

17. Davydova M. A., Nefedov N. N. Existence and stability of the solutions with internal layers in multidimensional problems of the reaction-diffusion-advection type with balanced nonlinearity. Modelirovaniye i analiz informatsionnykh sistem = Modeling and Analysis of Information Systems. 2017; 24(1):31-38. DOI: https://doi.org/10.18255/1818-1015-2017-1-31-38 (In Russ.)

18. Antipov E. A., Levashova N. T., Nefedov N. N. Asymptotics of the front motion in the reaction-diffusion-advection problem. Computational Mathematics and Mathematical Physics Journal. 2014; 54(10):1536-1549. DOI: https://doi.org/10.1134/S0965542514100029

19. Nefedov N. N. General scheme of asymptotic investigation of stable contrast structures. Nelineynaya dinamika = Nonlinear Dynamics. 2010; 6(1):181-186. Available at: http://mi.mathnet.ru/ nd65 (In Russ.)

20. Efstathiou C., Luhar M. Mean turbulence statistics in boundary layers over high-porosity foams. Journal of Fluid Mechanics. 2018; 841:351-379. DOI: https://doi.org/10.1017/jfm.2018.57

21. Squire D. T., Morrill-Winter C., Hutchins N., Schultz M. P., Klewicki J. C., Marusic I. Comparison of turbulent boundary layers over smooth and rough surfaces up to high Reynolds numbers. Journal of Fluid Mechanics. 2016; 795:210-240. DOI: https://doi.org/10.1017/jfm.2016.196

22. Swaters G. E. Internal dissipative boundary layers in the cross-equatorial flow of a grounded deep western boundary current. Geophysical & Astrophysical Fluid Dynamics. 2017; 111(2):91-114. DOI: https://doi.org/10.1080/03091929.2017.1287909

23. Kumar D. A parameter-uniform method for singularly perturbed turning point problems exhibiting interior or twin boundary layers. International Journal of Computer Mathematics. 2018, p. 1-18. DOI: https://doi.org/10.1080/00207160.2018.1458098

24. Xu H., Jin Y. The contrast structures for a class of singularly perturbed systems with het-eroclinic orbits. Discrete Dynamics in Nature and Society. 2016. Article ID 6405853. DOI: http:// dx.doi.org/10.1155/2016/6405853

25. Belov A. A., Kalitkin N. N., Kuzmina L. V. Modeling of chemical kinetics in gases. Mathematical Models and Computer Simulations. 2017; 9(1):24-39. DOI: https://doi.org/10.1134/S2070048217010057

26. Lahaye M. E. Une metode de resolution d'une categorie d'equations transcendentes. Compter Rendus Hebdomataires des Seances de L'Academie des Sciences. 1934; 198(21):1840-1842.

27. Lahaye M. E. Solution of system of transcendental equations. Académie Royale de Belgique. Bulletin de la Classe des Sciences. 1948; 5:805-822.

28. Davidenko D. F. On a new numerical method for solving systems of nonlinear equations. Doklady Akademii nauk SSSR = Reports of the Academy of Science of the USSR. 1953; 88(4):601-602. (In Russ.)

29. Davidenko D. F. Approximate solution of nonlinear equation systems. Ukrainskiy matematicheskiy zhurnal = Ukrainian Mathematical Journal. 1953; 5(2):196-206. (In Russ.)

30. Vorovich I. I., Zipalova V. F. On solution of nonlinear boundary value elasticity problems by method of transformation to Cauchy problem. Prikladnaya matematika i mekhanika = Applied Mathematics and Mechanics. 1965; 29(5):894-901. (In Russ.)

31. Riks E. The application of Newton's method to the problem of elastic stability. Journal of Applied Mechanics. 1972; 39(4):1060-1065. Available at: http://appliedmechanics.asmedigitalcollection.asme.org/ article.aspx?articleid=1400638

32. Kuznetsov E. B., Leonov S. S. Parametrization of the Cauchy problem for systems of ordinary differential equations with limiting singular points. Computational Mathematics and Mathematical Physics Journal. 2017; 57(6):931-952. DOI: https://doi.org/10.1134/S0965542517060094

33. Kuznetsov E. B., Leonov S. S. Examples of parametrization of the Cauchy problem for systems of ordinary differential equations with limiting singular points. Zhurnal vychislitelnoy matematiki i matematiches-koy fiziki = Computational Mathematics and Mathematical Physics Journal. 2018; 58(6):914-933. DOI: https://doi.org/10.7868/S0044466918060066 (In Russ.)

34. Semenov A. A. Strength and stability of geometrically nonlinear orthotropic shell structures. Thin-Walled Structures. 2016; 106:428-436. DOI: https://doi.org/10.1016Zj.tws.2016.05.018

35. May S., Vignollet J., de Borst R. A new arc-length control method based on the rates of the internal and the dissipated energy. Engineering Computations. 2016; 33(1):100-115. DOI: https://doi.org/10.1108/ EC-02-2015-0044

36. Wang X., Ma T.-B., Ren H.-L., Ning J.-G. A local pseudo arc-length method for hyperbolic conservation laws. Acta Mechanica Sinica. 2015; 30(6):956-965. DOI: https://doi.org/10.1007/s10409-014-0091-0

37. Kalitkin N. N., Poshivaylo I. P. Computations with inverse Runge-Kutta schemes. Mathematical Models and Computer Simulations. 2014; 6(3):272-285. DOI: https://doi.org/10.1134/S2070048214030077

38. Rosenbrock H. H. Some general implicit processes for the numerical solution of differential equations. The Computer Journal. 1963; 5(4):329-330. DOI: https://doi.org/10.1093/comjnl/5A329

Received 05.07.2018; revised 03.09.2018; published online 28.12.2018

About authors:

Evgenii B. Kuznetsov, Professor, Department of Modelling of Dynamic Systems, Moscow Aviation Institute (National Research University) (4 Volokolamskoe Shosse, Moscow 125993, Russia), D.Sc. (Physics and Mathematics), ResearcherlD: S-9576-2018, ORCID: https://orcid.org/0000-0002-9452-6577, [email protected]

Sergey S. Leonov, Associate Professor, Department of Modelling of Dynamic Systems, Moscow Aviation Institute (National Research University) (4 Volokolamskoe Shosse, Moscow 125993, Russia), Ph.D. (Physics and Mathematics), ResearcherlD: E-7641-2017, ORCID: https://orcid.org/0000-0001-6077-0435, [email protected]

Ekaterina D. Tsapko, Student, Department of Modelling of Dynamic Systems, Moscow Aviation Institute (National Research University) (4 Volokolamskoe Shosse, Moscow 125993, Russia), ResearcherID: S-8613-2018, ORCID: https://orcid.org/0000-0002-4215-3510, [email protected]

Contribution of the authors:

E. B. Kuznetsov - conducting theoretical research on the method of solution continuation with respect to a parameter, formulation of the basic research concepts, final editing the text; S. S. Leonov - reviewing and analysing the literature, development of computational algorithms of the solution continuation method with respect to the best argument, formulation of conclusions, editing the text; E. D. Tsapko - development of software and numerical experiments, writing the draft, word processing.

All authors have read and approved the final version of the paper.