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

СПИСОК ЛИТЕРАТУРЫ

1. Vafeiadis, V. A marriage of rely/guarantee and separation logic [TeKCT]/V. Vafeiadis, M. Parkinson//Proc. 18th Int'l Conf. on Concurrency Theory (CONCUR'07): Lecture Notes in Computer Science.-2007.-Vol. 4703. -P. 256-271.

2. Jones, C. Specification and design of (parallel) programs [TeKCT]/C. Jones//Proc. of IFIP'83.-North-Holland, 1983.-P. 321-332.

3. Feng, X. Local rely-guarantee reasoning [TeKCT]/X. Feng//ACM Symp. on Principles of Prog. Lang.-2009.-P. 315-327.

4. Feng, X. On the relationship between concurrent separation logic and assume-guarantee reasoning [TeKCT]/X. Feng, F. Rodrigo, S. Zhong//Proc. 16th European Symp. on Prog. (ESOP'07): Lecture Notes in Computer Science.- 2007.-Vol. 4421.-P. 173-188.

5. Reynolds, C. Separation logic: A logic for shared mutable data structures [TeKCT]/C. Reynolds//Proc. 17th Annual IEEE Symp. on Logic in Computer Science (LICS'02).-IEEE Computer Society, 2002.-P. 55-74.

6. Gotsman, A. Proving that non-blocking algorithms don't block [TeKCT]/A. Gotsman, B. Cook, M. Parkinson,

V. Vafeiadis//Proc. 36th ACM. Symp. on Principles of Prog. Lang. (POPL'09).-ACM Press, 2009.

7. Кларк, Э.М. Верификация моделей программ: Model Checking [Текст]/Э.М. Кларк, О. Грамберг, Д. Пелед; Под ред. Р. Смелянского.-М.: МЦНМО, 2002. -416 с.

8. Алешина, Н.А. Логика и компьютер. Моделирование рассуждений и проверка правильности программ [Текст]/Н.А. Алешина, А.М. Анисов, П.И. Быстров.-М.: Наука, 1990.-240 с.

9. Feng, X. Local rely-guarantee reasoning [Текст]/Х. Feng//ACM Symp. on Principles of Prog. Lang.-ACM, 2009.-P. 315-327.

10. Dodds, M. Deny-guarantee reasoning [Текст^/M. Dodds, X. Feng, M. Parkinson, V. Vafeiadis// Proc. ES0P.-2009.-P. 363-377.

11. Cliff, B.Jones. Wanted: a compositional approach to concurrency [Текст]/BJones, Cliff//In Prog. Methodology. -Springer-Verlag, 2003.-P 5-15.

12. Vafeiadis, V. Modular fine-grained concurrency verification: PhD thesis [Текст]/У. Vafeiadis.-Univ. of Cambridge. Computer Lab., 2008.

УДК 658.512:004.42,658.512:519.87

К.Г. Жуков, Д.Н. Бутусов

КОРРЕКЦИЯ ПОГРЕШНОСТИ РЕШЕНИЯ УРАВНЕНИЯ ГАРМОНИЧЕСКОГО ОСЦИЛЛЯТОРА МЕТОДОМ ПОСЛЕДОВАТЕЛЬНОГО ИНТЕГРИРОВАНИЯ

Широкий класс динамических систем описывается дифференциальными уравнениями второго порядка. Поэтому представляется целесообразным обеспечить на основе методики коррекции погрешности разомкнутых интеграторов [1] уменьшение величины погрешности интеграторов, охваченных отрицательной обратной связью.

Уравнение второго порядка в общем виде может быть записано следующим образом:

y + 2çray'+ra2y = га2x(t).

(1)

С нулевыми начальными условиями выражение (1) соответствует следующей передаточной функции:

Y (S ) га2

W (S ) =

X (S ) S2 + 2çraS +

(2)

В случае, когда х(?) = 0, уравнение (1) описывает свободные невынужденные колебания, а в случае, когда х(?) = 0 и с; = 0 - свободные незатухающие колебания.

Приведем уравнение (1) к виду:

y " =-2çray'-ra2 y + га2 x(t ), У "

(3)

(4)

— = -2qy -ay + ax(t). ю

Как правило, численные методы интегрирования требуют задания уравнений в нормальной форме Коши. Перейдем от уравнения (4) к системе уравнений:

[ у, = -2сасу1 - соу2 + шх(0; [у2 =®У1-

(5)

Рис. 1. Структурная схема решения системы (5)

Рис. 2. Структурная схема импульсной модели решения уравнения (4)

Для реализации модели (рис. 1) методом последовательного интегрирования (МПИ) необходимо на входах сумматора непрерывной модели установить ключевые элементы и фиксаторы нулевого уровня, тем самым осуществив переход к импульсной модели (рис. 2).

В работе [2] описана коррекция МПИ путем организации второго вычислительного процесса с определенным порядком вычислений. Второй процессор коммутируется в обратном порядке.

Установим окончательный порядок коммутации ключей для схемы решения уравнения (1) методом ПИ (рис. 3).

В соответствии с методикой, изложенной в работе [1], при реализации МПИ должна иметь место симметрия времен запаздывания ключевых элементов. При этом сумма времен запаздывания должна быть равна единице. Разделив промежуток интегрирования Т0 на восемь частей, получим временную диаграмму срабатывания импульсных элементов для каждого из процессоров (рис. 4).

х(с)

со

(0

гон-

гон-

0

0) 1 5 со

У\

(Я

гон

_1_

5

>'(0

Уг

Рис. 3. Порядок коммутации импульсных элементов

Рис. 4. Временная диаграмма срабатывания импульсных элементов

Частные случаи уравнения (1)

Случай свободных незатухающих колебаний. Уравнение гармонического осциллятора.

Когда х(^) = 0 и с; = 0 , уравнение (1) принимает вид:

у" + ш 2 у = 0. (6)

Перейдем к нормальной форме Коши:

У1 = -®У 2

У 2 =

(7)

при у (0) = 1 и у 2 (0) = 0.

Найдем корни характеристического полинома для (6):

г2 + ю2 = 0, г12 = 0 ± ш/.

Аналитическое решение уравнения (6) при заданных начальных условиях

(8)

y(t) = sin(®t).

(9)

Построим Simulink-модель согласно рис. 3. Для обеспечения максимальной точности применяемых «непрерывных» интеграторов выберем самый точный из имеющихся в инструментальной системе численных методов - метод ODE 8 (Dormand-Prince). Шаг интегрирования зададим равным 0,001 с при времени решения равном 5 с. Для удобства анализа результатов моделирования возьмем ю = 2п.

На рис. 5 показаны графики погрешностей вычисления переменной y2 в каждом из процессоров. На рис. 6 изображен график погрешности усредненного значения в сравнении с погрешностью метода Рунге-Кутта второго порядка. В табл. 1 приведены результаты компьютерного моделирования.

Скорректируем решение системы (7) согласно методике, изложенной в работе [1]. Для этого вычислим коэффициент коррекции по формуле:

2sin(--)

K

2

шГ0 cos(qT0 (0,5 - Л))

(10)

Таблица 1

Погрешности решения (7) методом ПИ без коррекции

Имя переменной Максимальная погрешность решения в первом процессоре Максимальная погрешность решения во втором процессоре Максимальная погрешность усредненного значения

У1 0,003139024288966 0,003190701642487 4,91208786803* 10"5

У 2 5,167761231231*10"5 5,167761231231*10~5 5,167761231231*10~5

Таблица 2 Погрешности решения (7) методом ПИ с одной коррекцией

Имя переменной Максимальная погрешность решения первого процессора Максимальная погрешность решения второго процессора Максимальная погрешность усредненного значения

У\ 0,003141602989062 0,0031416029891 1,11022302462* 10"14

У 2 4,934822503077641 * 10"6 4,9348225030776* 10"6 4,934822503078* 10"6

4

0,125_0,25 _ 0,375 _ 0,5 _ 0,625

0

□

п

гп

+

+ п

ш

□

о

* гп

(10

! ! ! ! -Еггог у1 РР?1 -Еггог у2 РР2

"7"\..... / \ \ / \ / "Т"Х.....: / \ ! / \ 1 /\ / \

/ \ У / :/ \

1 1 1 1 : | / 1 1

\ 1 ! \ Г ! / /

\ 1 / 1 ! \ 1 \ / / \ /

\/ / \/ ; / I \ / ;.......... \ / \ /

1 ........... 1 ......'•:"...... 1 1

Рис. 5. Блок-диаграмма импульсной модели уравнения (7) и результаты ее исполнения

Подставим значение Д = 0,5 и, вычислив значение коэффициента, проведем коррекцию решения (рис. 5) путем изменения коэффициентов компьютерной модели. Результаты компьютерного моделирования системы (7) с коррекцией представлены в табл. 2.

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

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

1

К„ =

(11)

Уш

Л \

/ г \........... /V Д \ .....\.......... / А

N \ / /V ' 1 \ \ \ \ \ / А 1 \ \ \ м \ \ п \ 1 1 1

--\ / / V___/ \ \ ^^ \ 1 / / / \ \ / \ \ / / / / \ / / \/

\У\ \ / / \ \У \ / / / 1 у

............... : ............. \ 1 V/ \

Рис. 6. Графики погрешности вычисления усредненных значений у и у2 методом ПИ без коррекции

Таблица 3

Погрешности решения (7) методом ПИ с двумя коррекциями

Имя переменной Максимальная Максимальная Максимальная

погрешность решения в первом процессоре погрешность решения во втором процессоре погрешность усредненного значения

0,0031416029891 0,003141602989056 1,110223024625* 10"14

З'г 4,9348225031 * 10"6 4,934822503078* 10"6 9,436895709314* 10"15

Таблица 4

Погрешности решения системы (7) методом ПИ с коррекцией

Метод последовательного интегррфования

Имя переменной Максимальная Максимальная Максимальная

погрешность решения в первом процессоре погрешность решения во втором процессоре погрешность усредненного значения

4,93683237146492* 10"4 4,9368323714627* 10"4 3,7192471324943*10"15

0,031426266043352 0,031426266043355 5,4400928206633* 10"15

Метод Рунге-Кутта второго порядка

У1 0,020639061806329

У2 0,019654989833941

Метод Рунге-Кутта четвертого порядка

У1 2,843382012408802* 10"5

У 2 2,991302502796882*10"5

Рис. 7. Графики погрешности вычисления усредненных значений переменныхy иy2

Умножив выходную величину y2 на Kin, получим результаты, показанные на рис. 8.

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

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

накоплением результата [3]. Блок-диаграмма модели решения системы (7) с цифровыми интеграторами приведена на рис. 9.

Рассчитаем корректирующие коэффициенты для шага 0,01 и изменим параметры блоков. В табл. 4 приведены значения погрешностей решения системы (7) методом ПИ и методами Рунге-Кутта второго и четвертого порядков на шаге 0,01 с.

Рис. 8. Погрешность вычисления усредненного значения переменной у2 после второй коррекции

Рис. 9. Блок-диаграмма цифровой модели

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

Исследование системы, заданной передаточной функцией. Зависимые гармонические колебания

При с; = 0, уравнение (1) принимает вид:

у" + ю2 У = ю2 х(1). (12)

4

Рис. 10. Блок-диаграмма 81тиНпк-модели решения системы (13)

Таблица 5

Погрешности решения (13) методом ПИ с двумя коррекциями

Имя переменной Максимальная погрешность решения в первом процессоре Максимальная погрешность решения во втором процессоре Максимальная погрешность усредненного значения

Уг 4,9348224997* 10_6 4,93482249819* 10"6 7,12624403931 * 10"15

У 2 0,003141602989058 0,003141602989057 9,82547376793* 10~15

0.125 0,25 0,375 0.5 0.625 0.75 0.875

Рис. 12. Погрешность вычисления у1 и у2 методами ПИ и РК-4

В нормальной форме Коши:

У

Из табл. 6 видно, что предлагаемый способ коррекции справедлив и для случая неоднородного уравнения. Погрешность решения уравнения методом ПИ при тех же параметрах моделирования существенно меньше, чем у классических численных методов, используемых в инструментальных системах моделирования [4].

Подтверждена возможность коррекции погрешности метода ПИ при решении задачи Коши на примере гармонического осциллятора до величины, определяемой длиной машинного слова.

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

При тех же параметрах моделирования погрешность методов РК2 и РК4 составляет 10-2 - 10-6, погрешность метода ПИ практически отсутствует.

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

Таблица 6

Погрешности решения (7) различными численными методами

--+ юу = юх(У)

ю

Г у1 =-юу2 + юх(); (13)

1У2 = юУ .

Аналитическое решение уравнения (12) выглядит следующим образом:

у^) = 1 - совОю). (14)

Построим БтиНпк-модель импульсной системы по выражению (13) с теми же параметрами моделирования, что и у модели на рис. 5.

В табл. 5 приведены результаты моделирования рис. 10.

Реализуем модель с цифровыми (рис. 11) интеграторами и проведем сравнительную оценку точности метода ПИ и методов Рунге-Кутта второго порядка и Рунге-Кутта четвертого порядка при шаге интегрирования 0,01 с. Результаты моделирования сведены в табл. 6. Графики погрешностей приведены на рис. 12.

Метод последовательного интегрирования

Имя переменной Максимальная погрешность решения в первом процессоре Максимальная погрешность решения во втором процессоре Максимальная погрешность усредненного значения

У1 0,031426266043355 0,0314262660433 6,4392935428* 10"15

У2 4,93683237*10~4 4,936832372*10^ 4,4408920985* 10~15

Метод Рунге-Кутта второго порядка

У1 0,020664988026269

У2 0,019679575866569

Метод Рунге-Кутта четвертого порядка

У1 4,074510124584747*10"6

У2 3,884044443291934* 10"6

СПИСОК ЛИТЕРАТУРЫ

1. Жуков, К.Г. Анализ и коррекция погрешности цифровых интеграторов встраиваемых систем управления [Текст]/К.Г. Жуков, Д.Н. Бутусов/Научно-технические ведомости СПбГПУ-2009.-№ 6(91). -С. 17-25.

2. Жуков, К.Г. Методы и средства реализации

последовательно-параллельных интегрирующих струк-

тур: Дис. ... канд. техн. наук[Текст]/К.Г. Жуков.-Л., 1988.

3. Ту, Ю.Т. Цифровые и импульсные системы автоматического управления [Текст]/Ю.Т. Ту; Под. ред. В.В. Солодовникова.-М.: Машиностроение, 1964.

4. Жуков, К.Г. Анализ погрешности цифровых интеграторов в 8тиНпк [Текст]/К.Г. Жуков/Проектирование научных и инженерных приложений в среде МАГЬАБ: Тр. Всерос. науч. конф.-М.: ИПУ РАН, 2004.