CC BY
41
УДК 519.622
© Е.А. Новиков, М.В. Рыбков
ЧИСЛЕННЫЙ АЛГОРИТМ ПОСТРОЕНИЯ МНОГОЧЛЕНОВ УСТОЙЧИВОСТИ МЕТОДОВ ПЕРВОГО ПОРЯДКА
Работа выполнена при финансовой поддержке РФФИ (проект 14-01-00047)
Построен алгоритм получения коэффициентов многочленов устойчивости до степени m=27, соответствующих явным методам типа Рунге-Кутты первого порядка точности. Показано, что выбором значений многочлена в экстремальных точках можно повлиять на размер и форму области устойчивости. Приведены результаты расчетов.
Ключевые слова: жесткие задачи, явные методы, многочлены устойчивости.
E.A. Novikov, M.V. Rybkov
A NUMERICAL ALGORITHM FOR CONSTRUCTING POLYNOMIALS STABILITY FOR METHODS OF THE FIRST ORDER
An algorithm of stable polynomial coefficients obtaining up to degree m=27 for Runge-Kutta explicit methods of the first order of accuracy is constructed. The choice ofpolynomial's values at the points of extremum can influence on the size and shape of stability domain. The numerical results are submitted. Keywords: stiff problems, explicit methods, stability polynomials.
Введение
В ряде случаев для численного решения жестких задач применяются алгоритмы интегрирования на основе явных методов с расширенными областями устойчивости [1-3]. В [4] описан способ построения многочленов, обеспечивающих максимальную длину интервала устойчивости. С помощью алгоритма [5] получены коэффициенты многочленов устойчивости до степени m=13. Здесь разработан алгоритм получения коэффициентов многочленов устойчивости до степени m=27 для явных методов типа Рунге-Кутты первого порядка. Показано, что форма, размер и структура области устойчивости зависят от расположения корней многочлена устойчивости в комплексной плоскости.
1. Явные методы типа Рунге-Кутты
Для численного решения жестких задач вида
y' = f (t, у), y(0 = y>, f0 < t < tk,
где y и f - гладкие вещественные N-мерные вектор-функции, t - независимая переменная, в [2] предлагается применять явные методы вида
y n+1 = Уп + V ,11 pmiki , К = hf (tn +a,h Уп + v ) ,
где ki, 1<i<m, - стадии метода; h - шаг интегрирования; pmi, a. и fa - коэффициенты, определяющие свойства точности и устойчивости численной схемы. Для упрощения выкладок далее рассмотрим задачу Коши для автономной системы обыкновенных дифференциальных уравнений
У = f(у), у(0 = У0, f0 <t<tk, (1)
для решения которой применим методы вида
УЩ = Уп + V ^/=Д + 1,А , 1 < i < m - ^ У n+1 = Уп + V m=1 pmiki , (2)
где ki=hf(yni-i), 1<i<m, yn0=yn. Все полученные ниже результаты можно использовать для неавтономных задач, если положить
a1 = 0, a = V' 1 В , 2 < i < m .
1 ' i /—tj=V У '
Устойчивость одношаговых методов обычно исследуется на линейном скалярном уравнении y'=Xy, y(0)=y0, t>0 с комплексным числом X, Re(X)<0. Применяя вторую формулу (2) для решения y'=Xy, получим
У +1 = Q (z) У , Q (z) = 1+ V m c z , c . =V m b p ., 1 < i < m ,
s n+1 ■Z-'m^ ' у n ' * / j i=1 mi ' mi / j /=i ij r m j ' '
где г^Х. Отсюда следует, что функция устойчивости явного га-стадийного метода типа Рунге-Кутты представляет собой полином Qm(z) степени т. В [2] приведены условия порядка для методов вида (2). Метод будет иметь первый порядок точности, если рт1+...+ртт=ст1=1. Далее рассмотрим задачу получения таких коэффициентов многочлена устойчивости, чтобы область устойчивости схемы имела заданную, естественно «разумную», форму и размер.
2. Многочлены устойчивости на интервале [ут, 0]
Пусть заданы два числа k и т, к<т. Рассмотрим многочлены вида
О (х) = 1+ Ек е.х' + Е т с.х' , (3)
где коэффициенты с., 1<1<к заданы, а с., к+1<Кт, — свободные. Обычно с., 1<1<к определяются из требований аппроксимации. Поэтому для определенности ниже будем предполагать, что срМ!, 1<1 <к.
Обозначим экстремальные точки (3) через хь..., х„_1, причем х!>х2>...> хт_\. Неизвестные коэффициенты с., к+1«т определим из условия, чтобы многочлен (3) в экстремальных точках х., к<1<т—1 принимал заданные значения, то есть Qm,к(xI)=FI, к<.<т_1, где F(x) _ некоторая заданная функция FI=F(xI). Для этого на х., к<.<т—1 и су, к+1<у<т рассмотрим алгебраическую систему уравнений
отк (х)=F, о:., (х)=о, к <. < т —1, о:.,=1 Г=1 с*-1. (4)
Перепишем (4) в виде, удобном для расчетов на ЭВМ. Для этого обозначим через у, z, g и г векторы с координатами
у = хк+.-1, 2 . = ск+., = ^+.-1 — 1 — Еу=1 суУ. ,
г =—Еу=1 усМ—11 <.< т—к,
через Е\,..., Е5 — диагональные матрицы с элементами на диагонали вида
< = к +., ¿2 = 1/у., < =Ек=1 Усу—1 (к + У)2уУк+У—1,
<=Е к=2 у ( У—1Ну/—2+ Е ^ (к+У)(к+У—1) 2]ук;у—2 ,
= (—1)к+", 1 <|■< т — к,
а через А _ матрицу с элементами ау=у ¡к+], К у'<га_к. С использованием введенных обозначений задачу (4) можно записать в виде
Аг — § = 0, Е2АЕ1 г — г = 0 . (5)
Система (5) плохо обусловлена, что приводит к определенным трудностям при использовании для ее решения метода простой итерации. Для сходимости метода Ньютона требуются хорошие начальные условия, что в данном случае является трудновыполнимой задачей. Если в (4) положить ^=(_1У, к^Кга-1, то есть поставить задачу нахождения полинома с максимальным размером интервала устойчивости, то вопрос о вычислении начального условия у0 решается с использованием значений экстремальных точек многочлена Чебышева на отрезке [_2га2,0], где т _ степень полинома (3). Их можно вычислить по формуле
у{ = m2[cos(iж / т) — 1], 1 <. < т — 1. (6)
Подставляя (6) в первую систему (5), получим коэффициенты полинома Чебышева, для которого юи1(х)|<1 при хе[_2га2,0]. При любом к в качестве начальных условий можно взять (6) и, как показывают расчеты, имеется хорошая сходимость. Если же FIФ(_1)1, к^Кт— 1, то выбор начальных условий является, вообще говоря, искусством.
Опишем способ решения (5), который не нуждается в хороших начальных условиях. Для численного решения (5) используем метод установления, который заключается в том, что для стационарной задачи строится нестационарный процесс, решение которого устанавливается к решению исходной задачи. Итак, рассмотрим задачу Коши:
у' = Е5 (Е2 АЕ А'1 § — г), у(0) = у0. (7)
Ясно, что после определения стационарной точки (7) коэффициенты полинома устойчивости можно вычислить из первой системы (5). Заметим, что при использовании матрицы Е5 все собственные числа матрицы Якоби задачи (7) имеют отрицательные вещественные части, то есть задача (7) устойчивая. Из результатов расчетов следует, что (7) является жесткой задачей. Методы решения таких за-
дач предполагают вычисление матрицы Якоби, что в случае (7) связано с трудностями. Поэтому для ее решения используем метод второго порядка точности с численным вычислением и замораживанием матрицы Якоби [6].
3. Многочлены устойчивости на интервале [-1, 1]
Нетрудно показать, что с ростом т коэффициенты многочлена устойчивости стремятся к нулю. С использованием алгоритма [2] коэффициенты с, к+1<'<т получены до степени т=13. Рассмотрим алгоритм построения многочленов с заданными свойствами на промежутке [-1, 1]. В этом случае коэффициенты с, растут не так быстро, и можно построить многочлены при значениях т>13. Обозначим через |/т| длину интервала устойчивости т-стадийной явной формулы типа Рунге-Кутты, то есть на интервале [ут, 0] выполняется условие |От,к(х)<1. Тогда заменой переменных х=1-22/ут переведем интервал [/т, 0] в [-1, 1]. Получим многочлен
вт (*) = 1 т=0 . (8)
Тогда коэффициенты d', 0<'<т многочлена (8) связаны с коэффициентами с,, 0<'<т многочлена (3) связаны соотношением
с = и¥ё, (9)
где d=(d0,..., dm) , с=(с0,..., ст) , и - диагональная матрица с элементами на главной диагонали вида и''=(_2/ут)'_1, 1<'<т+1, а элементы V1 матрицы V задаются формулами
V1 = 1, 1 < 1 < т + 1; у'1 = V1,1 -1 + у1'"1,1'"1, 2<' < 1 < т +1; у'1 = 0, ' > 1.
Нетрудно видеть, что V представляет собой «треугольник Паскаля», элементы которого вычисляются по рекуррентной формуле. После построения многочлена (8) на интервале [-1,1], с помощью (9) легко вычислить коэффициенты (3).
Теперь перейдем к построению многочлена (8). Обозначим экстремальные точки (8) через Zl,..., zm_l, причем zl>z2>...>zm_l. Коэффициенты 0<'<т определим из условия, чтобы многочлен (8) в экстремальных точках 1<'<т-1 принимал заданные значения, то есть
О (г.) = Е, 1 <' < т -1,
где Е(2) _ некоторая заданная функция Е=Е(2'). Для этого на 1<'<т-1 и dj, 0<1<т рассмотрим алгебраическую систему уравнений
От (2,) = Е, О'т(2г) = 0, 1 <' < т -1; от (2) = 1 ,"=1 (10)
причем выполнены условия нормировки От(-1)=( -1)т и От(1)=1.
Перепишем (10) в виде, удобном для расчетов на ЭВМ. Для этого обозначим через у, w, g и г векторы с координатами
у1 = , г1 = 0, 1 < 1 < т -1; м>1 = di-1, 1 <' < т +1;
g¡ = Е, 1 <' < т -1; g¡ = 1, ' = т ; g1 = (-1)т , ' = т +1; через Е1 и Е2 _ соответственно матрицы размерности (т+1)х(т+1) и (т-1)х(т+1) с элементами на главной диагонали вида
г3! = 1 -1, 1 < 1 < т +1; в^ = 1/ уг, 1 <' < т -1, а через А _ матрицу размерности (т+1)х(т+1) с элементами
а1 = у/-1, 1 <' < т -1, 1 < 1 < т +1; ат,1 = 1, ат+1,1 = (-1) 1+1, 1 < 1 < т + 1.
Теперь задача (10) имеет вид
Aw - g = 0 , Е2AE1w - г = 0. (11)
Для численного решения (11) используется метод установления [2]. Определив коэффициенты многочлена (8), с помощью соотношения (9) вычисляем коэффициенты многочлена (3). Значение ут находим из условия, что искомый многочлен соответствует методу интегрирования первого порядка точности, то есть имеет место с1=1. Выписав вторую строку соотношения (11) и сделав необходимые преобразования, получим
Хт фЕт:^}/^ ^^ = -2т2.
4. Результаты расчетов
Из расчетов следует, что коэффициент ст многочлена (3) убывает с ростом т, и, в частности, при т=13 и к=1 величина ст порядка 10-26. Из-за ошибок округления продолжить расчет задачи (7) при т>13 проблематично. При расчетах задачи (11) установлено, что коэффициенты di, 0</<т многочлена (8) по абсолютной величине растут с ростом т и, в частности, при т=13 величина тах0<г-<т|^| порядка 105, а при т=25 - порядка 109, то есть коэффициенты di растут медленнее. Переход от коэффициентов многочлена (8) к коэффициентам (3) по формуле (9) происходит после решения задачи (11), что позволяет вычислить коэффициенты многочленов устойчивости до степени т=27.
Опишем влияние выбора значений функции F на размер и форму области устойчивости. Если положить ^=(-1)г, кЯ<т-1, то длина интервала устойчивости известна и определяется как |^т|=2т2. В этом случае для заданного т получаем максимальную длину области устойчивости метода вдоль вещественной оси. Область устойчивости таких методов является почти многосвязной, и в результате ошибок округлений, которые приводят к появлению небольших мнимых частей собственных чисел матрицы Якоби, область устойчивости сокращается.
Чтобы ошибки округления не приводили к сужению области устойчивости, ее нужно «растянуть» вдоль мнимой оси в точках соприкосновения эллипсов. Для этого можно положить Fi=(-1)1 /, 1</<т-1, 0<и<1. Расчеты показывают, что если выбрать /=0,9, то длина интервала устойчивости метода сокращается на 5-8% по сравнению с максимально возможной длиной, равной 2т2. При этом область устойчивости растягивается вдоль мнимой оси на указанных участках. Это обеспечивает лучшие свойства устойчивости метода к ошибкам округления при несущественном сокращении длины интервала устойчивости. Если же положить /=0,95, то длина интервала устойчивости сократится всего на 3-4%. Область устойчивости пятистадийного метода при ¡и=0,9 приведена на рис. 1. Длина интервала устойчивости такого метода составляет |^т|=46,79. Для лучшей визуализации корней многочлена (3) на всех рисунках в комплексной плоскости {Ы} приведены линии уровней |Qm,k(x)|=1, |Qm,k(x)|=0,8, |Qm,k(x)|=0,6, |Qm,k(x)|=0,4 и |Qm,k(x)|=0,2.
С уменьшением / от единицы до нуля корни многочлена (3) располагаются все ближе друг другу на вещественной оси, поэтому длина интервала устойчивости соответствующего метода закономерно уменьшается. Эллипсы, которые хорошо выражены при / =1, сближаются, не обеспечивая при этом достаточно большого растяжения области устойчивости вдоль мнимой оси. Поэтому в зависимости от решаемой задачи целесообразно использовать значение / в интервале от 0,8 до 0,95.
При решении задач, собственные числа матрицы Якоби которых имеют большие мнимые части и решения которых носят осциллирующий характер, часто не требуется значительное расширение интервала устойчивости. В этом случае шаг из условия точности выбирается достаточно малым, и поэтому расширение области устойчивости требуется в основном по мнимой оси. В случае наличия чисто мнимых собственных чисел нужно, чтобы на некотором участке мнимой оси выполнялось условие |Qm,k(x)|=1. При повышении порядка точности, то есть с ростом к, это условие выполняется са-
т=5, к=1, F={-0,9, 0,9, -0,9, 0,9}, Ы=46,79
В случае методов первого порядка точности, то есть при к=1, выполнения данного условия можно добиться выбором значений функции F. Например, при т=4, к=1, F={0,85, 0,95, 0,85} получим многочлен, удовлетворяющий этому условию. Поскольку т четное, и все значения функции F положительные, то график многочлена не пересекает ось абсцисс, а многочлен имеет две пары комплексно сопряженных корней. Поэтому область устойчивости растягивается вдоль мнимой оси, захватывая ее участок. При этом длина области устойчивости метода вдоль вещественной оси невелика и составляет |хт|=2,18. При уменьшении значений функции F длина области устойчивости вдоль действитель-
ной оси возрастает и, в частности, при Е={0,55, 0,65, 0,55} становится практически прямоугольной с |^„|=5,30. С дальнейшим уменьшением значений Fдлина интервала устойчивости |/т| также растет, но область устойчивости захватывает все меньший участок мнимой оси. Таким образом, при построении методов первого порядка для интегрирования задач, решения которых носят осциллирующий характер, разумно выбирать многочлены устойчивости, которые вблизи начала координат в комплексной плоскости {Ы} имеют пару комплексно сопряженных корней. При этом значения функции F, которым соответствуют эти корни, необходимо выбирать близкими к единице, чтобы область захватывала участок мнимой оси.
В случае нечетного т хотя бы один корень многочлена (3) будет действительный. Использовать это свойство можно различными способами. Например, если выбором значений функции F расположить действительный корень многочлена пятой степени посередине интервала, а остальные корни оставить комплексными, то можно построить область устойчивости в форме, близкой к прямоугольнику, вытянутому вдоль мнимой оси. Область устойчивости пятистадийного метода, многочлен устойчивости которого в экстремальных точках принимает значения F={0,2, 0,5, -0,5, -0,2}, приведена на рис. 2. Отметим, что чем ближе значения F2 и F3 к единице, тем больше получается длина интервала устойчивости, но область превращается в почти многосвязную.
Для значений функции F, влияющих на появление комплексно сопряженных корней, ситуация обратная - чем меньше Fl и тем длиннее интервал устойчивости метода. Однако комплексно сопряженные корни располагаются все ближе друг к другу и к действительной оси, что приводит к изменению формы области устойчивости - прямоугольник превращается в эллипс.
В общем случае, в зависимости от требуемых свойств многочлена, варьированием параметров т и k выбором значений функции F можно конструировать области устойчивости различной формы и размера.
Рис. 2. Область устойчивости при значениях параметров
т=5, £=1, F={0,2, 0,5, -0,5, -0,2}, Ы=17,21
Заключение
С помощью алгоритма построения многочленов с заданными свойствами на промежутке [-1, 1] получены коэффициенты многочленов устойчивости методов первого порядка до степени т=27. Показано, что выбором значений функции F можно влиять на размер и форму области устойчивости. Предложенный алгоритм позволяет повышать эффективность явных методов и создавать алгоритмы переменного порядка и шага для решения умеренно жестких задач.
Литература
1. Новиков Е.А., Шорников Ю.В. Компьютерное моделирование жестких гибридных систем. Новосибирск: Изд-во НГТУ, 2012. 451 с.
2. Новиков Е.А. Явные методы для жестких систем. Новосибирск: Наука, 1997. 195 с.
3. Новиков А.Е., Новиков Е.А. Алгоритм переменного порядка и шага на основе стадий метода Дорманда-Принса восьмого порядка точности // Вычислительные методы и программирование. 2007. Т.8, №2. С. 317-325.
4. Скворцов Л.М. Простой способ построения многочленов устойчивости для явных стабилизированных методов Рунге-Кутты // Матем. моделирование. 2011. Т.23, №1. С. 81-86.
5. Новиков Е.А., Рыбков М.В. Численный алгоритм конструирования областей устойчивости явных методов // Системы управления и информационные технологии. 2014. №1.1(55). С. 173-177
6. Новиков А.Е., Новиков Е.А. L-устойчивый (2,1)-метод решения жестких неавтономных задач // Вычислительные технологии. 2008. №13. С. 477-482.
Новиков Евгений Александрович, доктор физико-математических наук, профессор, главный научный сотрудник ИВМ СО РАН, заведующий кафедрой СФУ, тел. (391)2494724, е-mail novikov@icm.krasn.ru.
Рыбков Михаил Викторович, ассистент кафедры СФУ, тел. (905)0860115, E-mail: mixailrybkov@yandex.ru
Novikov Evgeny Alexandrovich, doctor of physical and mathematical sciences, professor, senior researcher, Institute of Computational Modeling SB RAS, head of the department, Siberian Federal University. Rybkov Mikhail Viktorovich, teaching assistant, Siberian Federal University.
УДК 517.948
© Г.А. Шишкин
КРАЕВЫЕ ЗАДАЧИ ИНТЕГРОДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ВОЛЬТЕРРА
ЗАПАЗДЫВАЮЩЕГО ТИПА
Используя новую модификацию функции гибкой структуры, в статье исследуется возможность решения краевых задач для линейных уравнений запаздывающего типа.
Ключевые слова: краевая задача, интегродифференциальные уравнения Вольтерра, разрешающее уравнение, функция гибкой структуры, запаздывающий тип уравнений.
G.A. Shishkin
THE BOUNDARY VALUE PROBLEMS OF VOLTERRA INTEGER-DIFFERENTIAL EQUATIONS OF RETARDING TYPE
In the article a possibility of solution the boundary value problems for linear equations of retarding type is studied using new updating of flexible structure function.
Keywords: boundary value problem, the Volterra integer-differential equations, resolving equation, function of flexible structure, retarding type of equations.
Введение
В работе [2] доказано, что задача Коши для всех интегродифференциальных уравнений Вольтерра с отклоняющимся аргументом запаздывающего типа с помощью функции гибкой структуры [1] преобразуется к разрешающему интегральному уравнению типа Вольтерра с обыкновенным аргументом, решение которого существует и притом единственное при выполнении условий ограниченности функций, входящих в уравнение. В данной работе исследуем вопрос о возможности аналогичных преобразований линейной краевой задачи для интегродифференциальных уравнений Вольтерра с отклоняющимся аргументом запаздывающего типа.
Постановка задачи и ее решение
Выпишем общий вид уравнений запаздывающего типа
' n4 f(x)y(')(u](x)) + a}k7(x,r)yil)(u](r))dr + h\Kn0(x,r)y(n\r)dr = f(x), (1)
У0«+ЕЕ
У=0 i=0
где и0 (х) = X, и; (х) < X У; = 1,1 и иу (х) х, функции (х), и■ (х) и / ( X ) - непрерывны, ядра К у ( х, ц ) - регулярны в квадрате а < х, Ц < Ь , с начальными функциями
У()(и] (х)) = у(г)(хоУг)(и; (х)), I = 0^, х е ЕХо, (2)
I
где Ех = ^Е]х , Е^ - множество точек, для которых соответствующие и ■ (х) < х при х > х0
0 ;=0 0 0
У = 1,1, а Е.°0 = [а, х0 ] и с линейными билокальными краевыми условиями
XКу0'^) + Рггy(i)(xj)] = ГТ, х=0,n-1, a < x0 < x, < b . (3)
