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

УДК 681.3

АЛГОРИТМ ПОНИЖЕНИЯ ПОРЯДКА РЕШАЕМОЙ СИСТЕМЫ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ С ПЕРЕМЕННЫМ ШАГОМ НА ПРИМЕРЕ

ПРЯМОЙ КИНЕТИЧЕСКОЙ ЗАДАЧИ

А.М. Белянин, С.Л. Подвальный

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

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

Введение

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

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

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

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

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

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

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

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

собственным числом матрицы Якоби исходной системы дифференциальных уравнений.

На других промежутках, называемых участками установления, производные от решения невелики, и поэтому шаг по точности может быть выбран достаточно большим. Длина отрезков установления значительно больше переходных участков. В случае расчетов с постоянным шагом его величина определяется переходным участком. В результате на интервале установления решения расчеты проводятся с маленьким шагом, хотя длина данного интервала может составлять более 90 % от всего промежутка

интегрирования. Это и является причиной низкой эффективности алгоритмов

интегрирования с постоянным шагом.

Белянин Алексей Михайлович - ВГТУ, аспирант, тел. (473) 243-77-18

Подвальный Семен Леонидович - ВГТУ, д-р техн. наук, профессор, тел. (473) 243-77-18

О 5 000 10 000 15 000 20 000 25 000 30 000 35 000 40 000 45 000 50 000 55 000 60 000 65 000 70 000 75 000 80 000 85 000 90 000 95 000

Пример нежёсткой задачи

На рисунке приведён пример нежёсткой системы уравнений (кинетический модуль со случайным обрывом) [4]. На рисунке отображены три монотонные кривые, которые являются моментами неактивных центров. Кривые с перегибом - это моменты активных центров. Для любых типов активных центров примерный вид кривых изменения моментов будет аналогичным.

Алгоритм понижения порядка решаемой системы ОДУ

Поскольку функции в системе обыкновенных дифференциальных уравнений изменяются с разной скоростью, имеет смысл более «медленные» функции считать постоянными на некотором отрезке. Таким образом понижается порядок решаемой системы ОДУ и уменьшаются вычислительные затраты. Очевидно, что если потратить «освобождённое» время на дополнительные вычисления, то решение системы будет более точным.

В явных одношаговых методах Рунге-Кутта разных порядков значение у 1+1 на

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

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

1. Вычисления векторов значений для

точки конца отрезка у 1+1 в количестве функций в системе, которые имеют различный шаг.

2. Хранение значений всех функций в точке начала отрезка, на котором приняли функцию за константу.

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

Пример использования алгоритма понижения порядка решаемой системы ОДУ

Описанный выше алгоритм был реализован и апробован на модели многоцентровой полимеризации [3].

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

скоростей реакций, протекающих в системе, и кинетической схемы процесса.

Система дифференциальных уравнений с моментами до 2-го порядка включительно имеет вид (кинетический модуль):

ёИ +,

+ -ш )С ац

— = -М I (к]р + кі )С/

-і і =1

-Л

■ = -Л Iк С

] Р ]

-і і =1

а^ ац

-Р

1 . . . . . _

—— = -крМР^ + (кім + ка Л) т 0 , і = 1, п,

А21 п і і і

—1 = і (кім + ка л )р/,

-і і =1

- т

аі

— = крМР-1 - (кім + ка Л )т0-

1

= к рм (Р^ +т — ) - (кім + кал )т(,

= к ірм (4Р1 + 2т 1 + т —) - (кім + ка л )т 2,

Аі

-1, п

—- = I (кім + ка Л )ті ,, = 0,1,2.

-І і =1

где кр , к і , ка — константы скоростей

реакций роста, передачи цепи на мономер и на АОС _)-го типа активных центров; М-мономер;

Р^ -активная частица, образованная на _)-м типе

активных центров; $ 1 - неактивная частица; А-

АОС; С 1ц - концентрация активных центров как сумма концентраций всех активных цепей в любой момент времени; ¡1 ,1 - моменты 8-го порядка активных (растущих) цепей полимера на _)-м типе активных центров, 1 - моменты 8-го

порядка неактивных цепей полимера; п- число типов активных центров.

Исходные данные для решения системы представляются в виде:

И (0) = И (0), А (0) = А (0),Р1(0) = Сац ;

е{0) = 0, т- (0) = 1- (0) = 0, - = 03

Рассмотрим процесс полимеризации бутадиена на каталитической системе К^13*3ТБФ-ТГА, параметры модели приведены в таблице 1, которая получена решением обратной кинетической задачи в предположении 4 активных центров [3]. В итоге кинетический модуль будет состоять из 22 уравнений.

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

Чем меньше количество одинаковых отрезков, на которые разбивается интервал, тем больше будет шаг интегрирования и тем меньше будет точность.В качестве эталона возьмём решение методом Рунге-Кутты четвёртого порядка при шаге Ь=0,00001.

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

1. Берём по 100 контрольных точек с эталона и решения исследуемого метода. Найдём значения относительных ошибок исследуемого метода в контрольных точка:

-71 от

ггил'1- _ уШшучешюе

^ ^| у« х ал иг..

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

Таблица 1

Кинетические параметры активных центров полимеризации бутадиена

Каталитическая система Тип АЦ л/моль* мин л/моль*мин ка, л/моль*мин Са, моль/л

ШСЬ*3ТБФ-ТГА I 32,5 0,11 0,9 2,2* 10-5

II 76 0,062 0,5 1,4* 10-5

III 211 0,043 0,22 5,3* 10-6

IV 625 0,031 0,07 1,4* 10-6

if-і

(У»

^эталэн _ полученное

ч

-yij

4-І

■^эталон

Ai

где п - общее число параметров системы, к - количество контрольных точек.

Результаты сравнение метода с переменным шагом и изменяющимся порядком явного метода Рунге-Кутты 4-го порядка с понижением и без понижения порядка решаемой системы обыкновенных дифференциальных уравнений представлены в таблице 2. В таблице 2 на пересечении строк и столбцов записано значение в конечной точке средней относительной ошибки по всем 22 функциям. Примерный вид кривых изменения параметров на каждом активном центре соответствует иллюстрированному графику рис.1.

Рассмотренный алгоритм составляет один из модулей предлагаемой программной системы [5] и следует из принципа многоальтерна-тивности [6], в соответствии с которым необходимо иметь набор методов и алгоритмов решения дифференциальных уравнений.

1. Новиков Е.А. Алгоритм интегрирования переменной структуры для решения жестких задач на основе явного и Ь-устойчивого методов. // Вестник СибГАУ, 1(2008). № 18. С. 75-78.

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

3. Усманов Т.С., Спивак С.И., Усманов С.М. Обратные задачи формирования молекулярно-массовых распределений. — М.: Химия, 2004. — 252 с.

4. Подвальный С.Л. Моделирование промышленных процессов полимеризации. М.: Химия. 1979.-256 с.

5. Подвальный С. Л. Методы решения прямой и обратной кинетических задач в зависимости от сложности химической системы / С.Л. Подвальный, А.М. Белянин, А. В. Плотников // Вестник Воронежского государственного технического университета. 2012. Т.8. № 5. С. 18-21. Подвальный С. Л. Многоальтернативные системы: обзор и классификация / С. Л. Подвальный // Системы управления и информационные технологии. 2012. Т. 48. № 2. С. 4-13.

Воронежский государственный технический университет

Таблица 2

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

Время, мс метод 250 400 1300 2400 4000

Без понижения порядка 0,007945 0,00099 5,474E-05 7,93E-06 1,07E-05

С понижением порядка 0,000347 3,72E-05 1,144E-05 5,176E-06 4E-06

Выводы

Таким образом, разработанный алгоритм решения систем обыкновенных дифференциальных уравнений, отличающийся возможностью понижения порядка решаемой системы за счёт принятия за константы на некоторых отрезках функций с наименьшими значениями по модулю производных, позволяет снизить машинное время вычислений. Из таблицы 2 видно, что «сэкономленное» время было потрачено на дополнительные вычисления, что дало прирост в точности.

Литература

ALGORITHM FOR REDUCTION OF

ORDER TO SOLVE SYSTEMS OF

ORDINARY DIFFERENTIAL EQUATIONS WITH VARIABLE STEP AN EXAMPLE DIRECT KINETIC PROBLEM

A.M. Belianin, S.L. Podvalny

In the article the order reduction algorithm solved a system of ordinary differential equations. Results of comparing the explicit Runge-Kutta method of order 4 with the Runge-Kutta 4th order using the order reduction solved a system of ordinary differential equations

Key words: ordinary differential equations, Runge-Kutta method, a variable pitch, non-rigid system