Контроль устойчивости метода Фельберга седьмого порядка точности Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 11, № 4, 2006

КОНТРОЛЬ УСТОЙЧИВОСТИ МЕТОДА ФЕЛЬБЕРГА СЕДЬМОГО ПОРЯДКА ТОЧНОСТИ*

E. A. Новиков Институт вычислительного моделирования СО РАН

Красноярск, Россия e-mail: [email protected]

Ю. В. ШОРНИКОВ Новосибирский государственный технический университет, Россия

An inequality for the stability control of Felberg method of 7th order of accuracy is formulated, which is based on the evaluation of maximum eigenvalue of Jacoby matrix by the power-law method. This evaluation does not increase the number of calculations of the right-hand side of ODE. The results of computations are provided, which confirm the almost double increase of efficiency of the method due to the additional stability control.

Введение

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 05-01-00579-а).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

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

В настоящее время можно выделить два подхода к контролю устойчивости [1, 2]. Первый связан с оценкой максимального собственного числа матрицы Якоби через ее норму с последующим контролем (наряду с точностью) неравенства Н||/у|| < Д [1], где положительная постоянная Д связана с размером области устойчивости. Ясно, что для явных методов, где матрица Якоби не участвует в вычислительном процессе, это приводит дополнительно к ее нахождению и, следовательно, к значительному увеличению вычислительных затрат. Второй подход основан на оценке максимального собственного числа матрицы Якоби степенным методом через приращения правой части системы дифференциальных уравнений с последующим контролем неравенства Н|Атах| < Д [2]. Во всех рассмотренных авторами ситуациях такая оценка фактически не приводит к увеличению вычислительных затрат [3].

В настоящей статье с применением предложенного в [2] способа оценки максимального собственного числа матрицы Якоби построено неравенство для контроля устойчивости метода Фельберга седьмого порядка точности. Приведены результаты расчетов, подтверждающие повышение эффективности алгоритма интегрирования почти в два раза за счет контроля устойчивости.

1. Метод Фельберга

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

у' = /(¿,у), у(*о) = Уо, ¿0 < г < гк, (1)

где у и / — Ж-мерные вещественные вектор-функции; г — независимая переменная. Для решения (1) будем использовать явные формулы типа Рунге — Кутты следующего вида:

13 / г—1

Уп+1 = Уп + ^Ртгкг, кг = Н/ I Ьп + агк, Уп + ^ А,к,

V ,=1

г=1

где Н — шаг интегрирования а1

2 115 15

0, а2 = —, а3 = -, а4 = -, а5 = —, а6 = -, а7 = -, ' 2 27 3 9' 4 6' 5 12' 6 2' 7 6'

аз

1 6;

ад

2

3'

2

а 10 1

3, ап = 1, а12 = 0, а1з = 1,

1

1

1

^21 = 77, ^31 = 777. , ^32 = Т7Т, ^41 — 77, ^42 = 0, ^43 = X

8'

0,

65 27,

27 36' ' 12' ' 24'

5 25 25 1

в51 =12 , ^52 = 0 в53 = —16, ^54 = 16, в61 = 20, вб2 = вб3 11 25 125

в64 = 7, в65 = 7, в71 = — ^77, в72 = в73 = ° в74 = 777, /$75

4 5 108 108

125 31 61

в76 = , в81 = , в82 = в83 = в84 = 0, ^85 = ^^ , в86 = — 9 ,

в87 = 900, в91 = 2, в92 = в93 = 0, в94 = 108, в95 = 704, вд6 = 9", 67 91 23

в97 = 77Т, в98 = 3, ^10,1 = — ^¡77, ^10,2 = ^10,3 = 0, ^10,4 = 777 ,

2

90

108

108

(2)

^10,5 = - 935, Ао.е = 3541, Ао,7 = -, ^ю.з = "6", Ао,9 = - ,

в11>1 = Ш' в11'2 = в11'3 = 0' Ам = -И4' в11'5 = 1Ш' в11'6 = -: =2133 в = 45 в = 45 в = 18

в11'7 = 4100 ' в11'8 = 82 ' в11'9 = 164 ' в11'10 = 41 '

3 6 3

в12,1 = 205 ' в12,2 = в12,3 = в12,4 = в12,5 = 0 , /$12,6 = — 41 , /$12,7 = — 205 '

„ _ 3 _ 3 _6 а 1777

в12,8 = — 77) в12,9 — 77) в12,10 — 77) в12,11 = 0 ) Р13,1

41 41 41 4100

341 4496 289

в13,2 = /#13,3 = 0 , в13,4 = —164 ' в13,5 = 1025' в13'6 = '

в13,7 = — 4190 ' в13,8 = 52 ' в13,9 = 164 ' в13,10 = 44 ' /#13,11 = 0 ' /#13,12 = 1- (3) При значениях коэффициентов

41 34 9

Р71 = 840 ' Р72 = Р73 = Р74 = Р75 = 0 , Р76 = 4^ , Р77 = Р78 = 35 ,

9 41

Р79 = Р7,10 = 280 ' Р7,11 = 840 ' Р7,12 = Р7,13 = 0 (4)

схема (2), (3) имеет седьмой порядок точности [4].

2. Контроль точности и устойчивости

Согласно [4] численная формула (2), (3) с коэффициентами

34 9

Р81 = Р82 = Р83 = Р84 = Р85 = 0 ' Р86 = —' Р87 = Р88 = — '

105 35

9 41

Р89 = Р8,10 = 280 ' Р8,11 = 0 ' Р8,12 = Р8,13 = 840 (5)

имеет восьмой порядок точности. Тогда локальную ошибку 8п метода седьмого порядка можно оценить по формуле

13

¿п = — Р7г)кг- (6)

->п

г=1

В результате для контроля точности вычислений применяется неравенство

1М< ^ ' (7)

где || ■ || — некоторая норма в Ки; е — требуемая точность расчетов. Учитывая, что = 0(й8), шаг йас по точности выбирается по формуле

йас = дй ' (8)

где д находится из уравнения

?8|к11 = е.

Если q < 1, то происходит повторное вычисление решения (возврат) с шагом h, равным qh. В противном случае вычисляется приближенное решение, а прогнозируемый шаг находится по формуле (8). Неравенство (7) хорошо зарекомендовало себя при решении многих практических задач, оно используется и в данной работе. В дальнейшем алгоритм переменного шага на основе схемы (2) с коэффициентами (3)-(5) и неравенством для контроля точности (7) будем называть FEL78.

Чтобы построить неравенство для контроля устойчивости, применим (2) для решения линейной задачи y' = Ay с постоянной матрицей A. Первые три стадии ki, k2 и k3 схемы (2) применительно к данной задаче имеют вид

ki = Xy„, k2 = (X + 27X2) yn, кз = (X + 1X2 + ^X^ y„,

где X = hA. Нетрудно видеть, что имеют место соотношения

12кз - 18к2 + 6ki = 27X3y„, к2 - ki = 27X2y„.

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

|(12кз - 18к2 + 6ki)j| vn = max —-—-——-—. (9)

П 1<j<N |(k2 - ki)j |

Тогда согласно [3] для контроля устойчивости метода Фельберга можно применять неравенство

vn < D, (10)

где постоянная D ограничивает интервал устойчивости. Устойчивость методов типа Рун-ге — Кутты обычно исследуется на скалярном тестовом уравнении

y' = Ay, Re (А) < 0. (11)

Применяя (2) с коэффициентами (4) для решения (11), получим, что функция устойчивости Q7(x) метода седьмого порядка точности имеет вид

ii

Q7(x) = 1 + ^^ c7)ixl, x = hA,

i=1

где

с7>1 = 1, с7,2 = 0.5, с7,3 = 0.16666666666667,

с7,4 = 0.41666666666667 ■ 10-1, с7,5 = 0.83333333333333 ■ 10-2,

С7,6 = 0.13888888888889 ■ 10-2, 07,7 = 0.19841269841270 ■ 10-3,

с7,8 = 0.23165371472663 ■ 10-4, с7,9 = 0.23671439526314 ■ 10-5,

07,10 = 0.51829448771964 ■ 10-7, 07,11 = -0.43191207309970 ■ 10-7.

С использованием [3] нетрудно убедиться, что интервал устойчивости метода седьмого порядка приблизительно равен 5. Функция устойчивости ^8(ж) метода восьмого порядка точности имеет вид

i2

Qs (x) = 1 + 5^ cg.ix", x = hA, i=1

где

с8)1 = 1 ' с8,2 = 0.5 ' с8,3 = 0.16666666666667'

С8,4 = 0.41666666666667 ■ 10-1 ' С8,5 = 0.83333333333333 ■ 10-2 '

с8,6 = 0.13888888888889 ■ 10-2 ' с8,7 = 0.19841269841270 ■ 10-3 '

с8,8 = 0.24801587301587 ■ 10-4 ' с8,9 = 0.23490700935724 ■ 10-5 '

С8,10 = 0.23620053064283 ■ 10-6 ' с&>п = —0.25914724385982 ■ 10-7 ' с8,12 = —0.14397069103323 ■ 10-7.

Интервал устойчивости метода восьмого порядка также приблизительно равен 5. Поэтому в неравенстве (10) положим О = 5. Учитывая, что уп = О(й), шаг й5* по устойчивости можно выбирать по формуле й5* = гй, где г вычисляется из равенства ггп = О.

Оценка (9) является грубой, потому что: 1) вовсе не обязательно, чтобы максимальное собственное число было сильно отделено от остальных; 2) в степенном методе применяется мало итераций; 3) дополнительные искажения вносит нелинейность задачи (1). Поэтому здесь контроль устойчивости используется как ограничитель на размер шага интегрирования. В результате прогнозируемый шаг вычисляется по формуле

й := тах[й ' тт(йас ' Л^)].

В дальнейшем алгоритм переменного шага с дополнительным контролем устойчивости численной схемы будем называть ЕЕЬ788Т.

3. Компьютерное моделирование

Новые алгоритмы анализа сложных динамических систем оказываются более эффективными для прикладного специалиста, если они окружены специализированными инструментальными средствами [6]. Индустрия создания таких средств моделирования становится сама по себе фундаментальной задачей исследования. Сейчас масштаб и объем трудностей при создании инструментальных средств настолько выросли, что можно говорить, что задача их преодоления сама стала задачей науки и представляет собой проблему фундаментального значения. [7].

Программная реализация алгоритмов с контролем точности и устойчивости встроена в библиотеку системы инструментальных средств машинного анализа ИСМА [8]. Класс исследуемых объектов в ИСМА ограничен системами обыкновенных дифференциальных уравнений с запаздыванием

х = /(х(г)' х(г — в) ' г)' г> 0 ' х(0) = х0 '

где х € Ям — вектор состояния; х(г) = ^(¿) при г € [—в' 0); г — независимая переменная; — т-мерная вектор-функция запаздывания, т < N; в = (т1 '... ' тт)т — вектор чистых запаздываний; / — нелинейная вектор-функция, удовлетворяющая условию Липшица; х0 = (х01 '... ' х0^)Т — вектор начальных условий.

Инструментальная среда ИСМА реализована с учетом простоты описания математических моделей и ориентирована на предметного пользователя. Переход от математической модели к ее компьютерному описанию происходит переписыванием правой части во встроенном редакторе ИСМА или в любом другом процессоре, допускающем текстовый

формат данных. Языковой процессор ИСМА интерпретирует введенные данные и придает им соответствующие типы, спецификации и т. д. В результате трудоемкость перехода к компьютерной модели существенно снижается. Графический интерпретатор GRIN предоставляет предметному пользователю широкомасштабную манипуляцию графическими данными (трассировку точек, мультиоконный режим отображения, конкатенацию графических объектов, редактирование текстовых полей, импорт данных и т.д.).

По сравнению с известными отечественными (MVS, AnyLogic) и зарубежными (Simulink, Modélica) системами моделирования [9] ИСМА отличается простотой и естественной формой представления обозначенного класса моделей и не требует от предметного пользователя дополнительных знаний в области программирования, а также современных парадигм (UML, ООП). Система ИСМА имеет разнообразную библиотеку традиционных и оригинальных численных методов анализа жестких и нежестких систем и поэтому в списке современных инструментов визуального моделирования предоставляет пользователю свои функциональные преимущества.

4. Результаты расчетов

Расчеты проводились на IBM PC Athlon(tm)XP 2000 с двойной точностью в среде ИСМА. В конкретных расчетах левая часть неравенства (7) вычислялась по формуле

I I

II ¿nil = max 1 ¿n|

i<j<N |уП | + r

где j — номер компоненты, r — положительный параметр. Если по j-й компоненте решения выполняется неравенство |yn | < r, то контролируется абсолютная ошибка er, в противном случае — относительная ошибка е. В расчетах параметр r выбирался таким образом, чтобы по всем компонентам решения фактическая точность была не хуже задаваемой.

Ниже через isa и iwo обозначены соответственно суммарное число шагов интегрирования и число повторных вычислений решения (возвратов) вследствие нарушения требуемой точности расчетов. Число вычислений правых частей ifu системы (1) можно определить по формуле ifu = 13 • isa +12 • iwo. Пример 1 [5, с. 145-146].

y1 = 2tyiy4, у2 = 10ty5 ш,

y3 = 2ty4, (12)

У4 = -2t(y3 - 1), t G [0,15п], yi (0)= У2(0) = Уз(0)= ш(0) = 1, ho = 10-2.

Точное решение (12) имеет вид

y1(t) = exp(sin t2), y2(t) = exp(5 sin t2), y3(t) = sin t2 + 1, y4(t) = cos t2.

Целью расчетов было определение влияния неравенства для контроля устойчивости на вычислительные затраты при решении нежестких задач. Из результатов расчетов с точностью е = 10-6 следует, что при вычислении решения (12) с помощью алгоритма FEL78 без контроля устойчивости затраты isa = 4055 и iwo = 1750, а для алгоритма FEL78ST — isa = 4094, iwo = 1554.

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

Пример 2 [10].

у' = -0.013yi - lOÜOyiys,

у2 = -2500у2Уз , (13)

у3 = -0.013yi - 1000yiy3 - 2500у2уз ,

t G [0, 50], yi(0) = 1, y2(0) = 1, Уз(0) = 0, ho = 2.9 ■ 10-4.

Данная задача является жесткой. Из результатов расчетов с точностью е = 10-6 следует, что при вычислении решения (13) алгоритмом FEL78 без контроля устойчивости затраты isa = 37 785, iwo = 37 752, ifu = 950 860, т.е. почти каждый шаг сопровождается повторным вычислением решения. Это естественно для всех алгоритмов без контроля устойчивости. Для алгоритма FEL78ST затраты таковы: isa = 37876, iwo = 454, ifu = 49 7836. Из сравнения результатов расчетов следует, что контроль устойчивости приводит к повышению эффективности расчетов почти в два раза. В конце интервала интегрирования фактическая точность вычислений для алгоритма FEL78 на порядок лучше задаваемой, а для алгоритма FEL78ST — на два порядка. При более грубой точности расчетов эффективность FEL78ST возрастает, а для FEL78 остается примерно такой же. Эта тенденция сохраняется при интегрировании всех рассмотренных задач и, в частности, десяти жестких тестовых примеров [10]. Заметим, что для одного из лучших методов — DVERK78 системы Maple 9.5 — на примере (13) затраты ifu = 598 890.

Заключение

Использование неравенства для контроля устойчивости фактически не приводит к увеличению вычислительных затрат, потому что оценка максимального собственного числа матрицы Якоби системы (1) осуществляется через ранее вычисленные стадии и не вызывает рост числа вычислений функции f. Как правило, для этих целей применяются три первые стадии. Такая оценка получается грубой. Однако контроль устойчивости в качестве ограничителя на рост шага позволяет избежать негативных последствий грубости оценки. Более того, в некоторых случаях это приводит к нестандартно высокому повышению эффективности алгоритма интегрирования. На участке установления за счет контроля устойчивости старые ошибки стремятся к нулю, а новые невелики из-за малости производных решения. В некоторых случаях вместо оценки максимального собственного числа оценивается следующее по порядку. Шаг интегрирования становится больше максимально допустимого, и с таким шагом интегрирование осуществляется до тех пор, пока не нарушается неравенство для контроля точности. Число таких шагов невелико, однако величина шага может на несколько порядков превышать максимальный шаг по устойчивости. После нарушения неравенства для контроля точности шаг уменьшается до максимально возможного. Такой эффект может повторяться многократно в зависимости от длины участка установления. В результате средний шаг интегрирования может превышать максимально допустимый.

Список литературы

[1] Shampine L.M. Implementation of Rosenbrock methods // ACM Transaction on Mathematical Software. 1982. Vol. 8, N 5. P. 93-113.

[2] Новиков Е.А., Новиков В.А. Контроль устойчивости явных одношаговых методов интегрирования обыкновенных дифференциальных уравнений // Докл. АН СССР. 1984. Т. 277, № 5. С. 1058-1062.

[3] Новиков Е.А. Явные методы для жестких систем. Новосибирск: Наука, 1997. 195 с.

[4] Fehlberg E. Classical fifth-, sixth-, seventh- and eighth order Runge — Kutta formulas with step size control // Computing. 1969. Vol. 4. P. 93-106.

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

[6] Моисеев Н.Н. Математические задачи системного анализа. М.: Наука, 1981.

[7] Яненко Н.Н., Карначук В.И., Коновалов А.Н. Проблемы математической технологии // Числ. методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние ВЦ; ИТПМ. 1977. Т. 8, № 3. C. 129-157.

[8] Шорников Ю.В. и др. Инструментальные средства машинного анализа: Свидетельство официальной регистрации программы для ЭВМ № 2005610126. М.: Роспатент, 2005.

[9] БЕнькович Е.А., Колесов Ю.Б., Сениченков Ю.Б. Практическое моделирование динамических систем. СПб.: БХВ-Петербург, 2002. 464 с.

[10] Enright W.H., Hull T.E. Comparing numerical methods for the solutions of systems of ODE's // BIT. 1975. N 15. P. 10-48.

Поступила в редакцию 22 марта 2006 г., в переработанном виде — 3 мая 2006 г.