Алгоритм решения обыкновенных дифференциальных уравнений, комбинирующий а- и L-устойчивые методы различных порядков точности Текст научной статьи по специальности «Математика»

УДК 519.622

АЛГОРИТМ РЕШЕНИЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ, КОМБИНИРУЮЩИЙ А- И ¿-УСТОЙЧИВЫЕ МЕТОДЫ РАЗЛИЧНЫХ ПОРЯДКОВ ТОЧНОСТИ

© 2008 г. В.Н. Бирюков

Таганрогский технологический институт Южного федерального университета, 347928, г. Таганрог, ГСП-17А, пер. Некрасовский, 44, biryukov@users.tsure.ru

Taganrog Technological Institute of Southern Federal University, 347928, Taganrog, GSP-17A, Nekrasovskiy Lane, 44, biryukov@users.tsure.ru

Предложен алгоритм объединения двух А-устойчивых известных методов (более точного, но не обладающего L-устойчивостью, и менее точного, но обладающего ею) таким образом, что при малых шагах достигается точность первого метода, а при больших — устойчивость второго. Для анализа свойств алгоритма определялась интервальная погрешность численного решения тестовой задачи dy/dt = — y, t е[0, T], y(t0) = 1, T = tN = 216.

Ключевые слова: неявные методы Рунге-Кутты, А-устойчивость, L-устойчивость.

One-step A- and L-stable methods similar in design to the class of Runge-Kutta methods are developed for the efficient numerical integration of stiff systems of first-order nonlinear differential equations. The algorithms developed combine asymptotic 4-th and 6-th accuracy order with a region of absolute stability, which does not intersect the imaginary and the negative real axis.

Keywords: Runge-Kutta method, implicit method, A-stability, L-stability.

Основным свойством, затрудняющим анализ электронных цепей во временной области, обычно считается их жесткость, которая в простейшем случае определяется как отношение максимальной и минимальной постоянных времени цепи. Существующие ¿-устойчивые методы численного решения систем обыкновенных дифференциальных уравнений (СОДУ) позволяют эффективно анализировать подобные цепи за исключением случая высокой добротности цепи [1]. Матрица Якоби СОДУ такой цепи обладает почти мнимыми собственными значениями, и численный анализ в этом случае требует применения А-устойчивых методов решения СОДУ. К сожалению, фактический порядок точности известных А-устойчивых методов не превышает двух, так как методы 3-го порядка предназначены для решения только линейных СОДУ [2, 3].

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

На интервалах быстрого изменения движений жесткой системы ¿-устойчивость разностной схемы не очень важна, так как шаг решения при этом выбирается самый мелкий. При переходе к интервалу медленного изменения переменных именно ¿-устой-чивость определяет возможность быстрого увеличения шага решения. Очевидно, что сингулярную (сингулярно-возмущенную) задачу на разных участках интервала наблюдения можно решать разными методами: на участках медленного изменения движений использовать А- и ¿-устойчивый метод низкого порядка (с крупным шагом), а на участках быстрого изменения - А-устойчивый метод высокого порядка (с мелким шагом). При решении конкретной СОДУ изменение метода решения возможно «вручную» путем решения задачи с низкой точностью на всем интервале наблюдения ¿-устойчивым методом, выделения

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

Метод переменного порядка

В А-устойчивом методе экспоненциальной подгонки для задачи йх / й = /(х) , х(0) = х0 , 0 < t < Т разностная схема имеет вид

хп+1 = хп + И [(0,5-п)/„ + (0,5 + ч)/п +!], (1)

где хп - сеточная функция; 0<т<0,5 - константа [2]. При ] — 0,5 (1) совпадает с разностной схемой ¿-устойчивого метода Эйлера 1-го порядка точности, при ] — 0 -с разностной схемой метода трапеций 2-го порядка, ¿-устойчивостью не обладающего. Если константу п в (1) заменить функцией шага, монотонно возрастающей от 0 при И = 0 до 0,5 при И — Итах(Итах < Т, обычно Ищах ~ Т/50), то, как показано в [4], глобальная (интервальная) ошибка метода (1) близка к аналогичной ошибке метода трапеций при малых шагах и к ошибке метода Эйлера - при больших. Формально метод (1) при введении зависимости ] (И) имеет первый порядок точности, но фактически при шаге И < МШ(Яе )-1 - второй.

В качестве функции г* = /(к) с целью экономии вычислительных ресурсов целесообразно использовать рациональную функцию с целочисленными коэффициентами, например, степенной многочлен

* "

Г „(к) = 2 а (к / ктах)' , коэффициенты а которого опре-

i = 1

деляются условиями: i|*(hmax) = 0,5, dki / dhk\

= 0,

к = 1,2,..., п — 1. Решение тестовых задач показало, что порядок полинома можно выбрать из соотношения Г> ^Отах/Ттт) ■

Область устойчивости предлагаемого метода при к ^ 0 стремится к области устойчивости метода трапеций, а при к ^ - к области устойчивости неявного

метода Эйлера. Следовательно, метод -устойчив. Граница области устойчивости в верхней комплексной

полуплоскости приведена на рис. 1 для ряда текущих

*

значений п ■

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

-5 -10

- 15

- 20

- 25 30 Ь

- 35

Im(hA-) -л

■ 1 0 1 Fte(hA)

Рис. 1. Граница области устойчивости метода (1) при различных значениях п

Для подтверждения высокой устойчивости метода (1) при переменном коэффициенте г* (к) в области больших шагов на рис. 2 приведены зависимости от шага глобальной погрешности Д = | хм — х(* м)| решения задачи / & = —у , * е [0, Т], у(*0 ) = 1, Т = *ы = 216 для трех методов.

При больших шагах погрешность (1) с п (Л) близка к погрешности неявного метода Эйлера, а при малых - к погрешности неявного метода трапеций, что позволяет называть рассматриваемый метод методом переменного порядка. Отметим, что, во-первых, вид зависимости Д(к) не критичен к выбору порядка

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

Рассмотренный метод переменного порядка точности не имеет принципиальных преимуществ перед другими А -устойчивыми методами 2-го порядка, но

Рис. 2. Зависимость глобальной погрешности от шага: + - для неявного метода трапеций; □ - для неявного

метода Эйлера; о - для метода (1) при использовании

„ „

Гз (Л); ◊ - того же метода при г„ (Л)

Рассмотрим метод, состоящий из двух частичных шагов, на первом из которых используется неявная схема Эйлера, а на втором - схема (30) 4-го порядка точности из [5]:

хп + а = хп + к/(хп + а> *п + ак), (2)

хп + 1 = хп + а + [(1 — а)к/2]{/п + а + /п +1 +

+4/[(хп + а + хп +1 )/2 — (1 — а)к(/п +1 — /п + а)/8, *п + к/2]} ■

Весовой коэффициент а здесь изменяется от 0 до 1, т.е. можно принять а* (к) = 2/*(к). Результаты расчета глобальной погрешности рассмотренной выше задачи показывают (рис. 3), что метод (2) имеет 3-й порядок точности, более низкий, чем 2-й вложенный в (2) метод.

-15 - 20

- 35

- 40

0,0 01 0,01

0,1

1,0

10,0

Рис. 3. Зависимость глобальной погрешности от шага: х - для метода Эйлера; + - для метода (30) из [5]; ◊ - для

метода (2) при а *(к) = 2/4 (Л); □ - для НМРК2; о - для метода (3) при а*(Л) = 2 г„ (Л)

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

Заменим неявную схему Эйлера в (2) полностью неявной схемой Рунге-Кутты 2-го порядка точности (НМРК2), обладающей ¿2-устойчивостью, максимальной для методов 2-го порядка точности. Разностная схема такого метода переменного порядка имеет вид

хп + а = хп + ИЛхп + а - (аИ/2)/п tn +аИ/2] хп + 1 = хп + а

+ [(1 -а)И/2]{/п

+ а ' J п +1 1

(3)

+4/[(хп + а + хп +1)/2-(1 -а)И(/п +1 -/п +а)/8,п + И/2]}.

Характер зависимости А(И) метода (3) (рис. 3) не изменяется по сравнению с методом (1): при малых шагах зависимость близка к соответствующей зависимости метода высокого порядка (т.е. метод (3) фактически имеет 4-й порядок точности), а при больших шагах - к методу низкого порядка (т.е. метод (3) обладает ¿2-устойчивостью).

Представляет интерес возможность получения метода более высокого порядка. Если в (3) на втором частичном шаге вместо А -устойчивого метода 4-го порядка использовать аналогичный метод (31) из [5] 6-го порядка, то новый метод, так же как и (2), не позволяет получить точность более высокого порядка (рис. 4).

1в А

- 5■

- 10■

- 15■

- 20■

- 25■

- 30■

- 35■

- 40

0,01 0,1 1,0 10,0 к

Рис. 4. Зависимость глобальной погрешности от шага: + - для метода 6-го порядка точности (31) [5];

□ - для НМРК2; ◊ - для метода переменного порядка при использовании а (к) = 2п 4 (к); х - при а 5 (к), о - при а 6 (к)

Тем не менее 6-й порядок точности может быть достигнут и в этом случае, если с уменьшением шага константы, определяющие разностные производные, будут стремиться к точным значениям быстрее, чем по линейному закону. Если потребовать квадратичную сходимость, то функция а* (И) может быть записана в виде

* п

а„(И) = ^Ь(И/Ищах) , где коэффициенты bi опреде-

I = 2

ляются из условий: а*(Ащах) = 1 , йка*/йИк\ = 0,

Итах

Поступила в редакцию_

к = 1,2,..., п - 2. Результаты, представленные на рис. 4, показывают, что применение последнего выражения действительно позволяет получить 6-й (эффективный) порядок точности, причем порядок функции а* (И) слабо влияет на устойчивость.

Обсуждение результатов и выводы

Комбинация двух разных методов с помощью весовых функций шага ]*п (к) и а*п (к) выбрана выше

для удобства анализа свойств полученного метода. При его практической реализации комбинация может быть любой удовлетворяющей сформулированным ограничениям. Предложенный алгоритм позволяет осуществлять выбор текущего (переменного) шага, используя только оценку текущей локальной погрешности. Следует отметить, что решение жестких задач требует переменного шага решения. Его применение неявно налагает ограничение на общее число шагов, так как алгоритм автоматического выбора шага обусловлен оценкой локальной погрешности решения, что при большом числе шагов делает оценку интервальной (глобальной) погрешности труднопредсказуемой. Из рассмотренных методов оптимальным соотношением между точностью, устойчивостью и арифметической сложностью обладает метод (3), поскольку состоит из двух методов одинаковой сложности. Сложность метода, полученного из (3) заменой 2-го вложенного метода методом 6-го порядка, определяется сложностью последнего.

Строго говоря, предложен не новый метод решения СОДУ, а алгоритм объединения двух известных методов (более точного, но не обладающего ¿-устой-чивостью, и менее точного, но обладающего ею) таким образом, что при малых шагах достигается точность первого метода, а при больших - устойчивость второго. С другой стороны, поскольку предлагаемая комбинация имеет свойства, которыми не обладает ни один метод в отдельности, то предложенный алгоритм можно считать и методом («метод Гира» разными авторами именуется и «алгоритмом Гира», и даже «программой Гира»).

Литература

1. Трохименко Я.К. Предисловие к тематическому выпуску «Автоматизация схемотехнического проектирования» // Изв. вузов. Радиоэлектроника. 1991. Т. 34. № 6. С. 3 - 4.

2. Современные численные методы решения обыкновенных дифференциальных уравнений / Под ред. Дж. Холла и Дж. Уатта: Пер. с англ. М., 1979.

3. Булатов М.В. Построение неклассических многошаговых схем для линейных ОДУ // Докл. АН 2005. Т. 404. № 1. С. 11-13

4. Бирюков В.Н., ПилипенкоА.М. Численный анализ жестких узкополосных систем // Радиосистемы. Радиоэлектронные устройства и системы управления, локации и связи. 2002. № 54. С. 36-41.

5. Ракитский Ю.В., Устинов С.М., Черноруцкий И.Г. Численные методы решения жестких систем. М., 1979.

24 марта 2008 г.