CC BY
36
ми. Это можно объяснить тем, что предложенный алгоритм учитывает свойства оптических частей микроскопа, поэтому сигнал является более “размытым”, чем те, которые были получены с помощью других подходов. Таким образом, результаты, полученные с использованием приведенного алгоритма, более напоминают регистрируемые экспериментальные данные. Сравнение результатов вычислений и экспериментальных данных будет приведено в следующих публикациях.
7. Выводы
Предложен алгоритм, позволяющий вычислять локальную интенсивность ЭХЛ свечения, регистрируемую конфокальным микроскопом. Данный подход является важным этапом в решении обрат -ной задачи, состоящей в восстановлении профиля концентрации синглет-частиц по регистрируемой локальной интенсивности ЭХЛ свечения, что позволит корректно сравнивать теоретические и экспериментальные данные, зная погрешность, допускаемую при измерении сигнала конфокальным микроскопом.
Признательность. С.н.с. Олейник А.И. благодарит мэрию г. Парижа за post-doc стипендию для работы в лаборатории академика Кристиана Аматора в Ecole Normale Superieure (Париж, Франция).
Литература: 1. FossetB., Amatore C.A., Bartelt J.E., Michael A. C, Wightman R.M. The use of Conformal Maps to Model the Voltammetric Response of Collector-Generator Double-Band Electrodes. // Analytical Chemistry. 1991.V. 63.P.306-314. 2. Amatore C, in:I. Rubinstein, (Ed.), Physical Electrochemistry: Principles, Methods and Applications (Chap. 4), Marsel Dekker, New York, 1995. 3. Amatore C., Bonhomme F, Bruneel J.L., Servant L, Thouin L. Mapping concentration profiles within the diffusion layer of an electrode. Part I: Confocal resonance Raman microscopy // Electrochemistry Communications. 2000. V.2. P.235239. 4. Amatore C, Bonhomme F, Bruneel J.L., Servant L, Thouin L. Mapping dynamic concentration profiles with micrometric resolution near an active microscopic surface by confocal resonance Raman microscopy. Application to diffusion near ultramicroelectrodes: First direct evidence for a conproportionation reaction // J. Electroanalytical Chemistry. 2000. V. 484. P. 1-17. 5. Bard A.J., FaulknerL.R. Electrochemical methods: fundamentals and applications.
УДК 517.958:541.14
АНАЛИЗ НОВЫХ МЕТОДОВ ЧИСЛЕННОГО РЕШЕНИЯ ЖЕСТКИХ ЗАДАЧ КОШИ
ШУЛЫК В.Н., КЛИМЕНКО А.В., СВИРЬ И.Б.
Рассматриваются новые методы решения жестких систем обыкновенных дифференциальных уравнений — метод квази-Рунге-Кутта и метод Алуффи-Пентини. Предлагается алгоритм адаптивного изменения шага интегрирования по времени. Методы сравниваются с известным методом Г ира на примере жесткой системы ОДУ, возникающей при математическом моделировании химической реакции Бриггса-Рошера.
John Wiley & Sons: N-Y. 2002. 834 p. 6. Свирь И.Б. Формулировка обобщенной краевой задачи для моделирования электрохемилюминесцентных процессов при различных режимах нестационарного электролиза // АСУ и приборы автоматики. 2003. № 122. С. 129134. 7. Лаврентьев М.А., Шабат Б.В. Методы теории функций комплексного переменного. М: Наука, 1973. 736с.8. http://www.physics.emory.edu/~weeks/confocal/
9. http://ru.ntmdt.ru/SPM-Techniques/Basics/3_SOM/
10. http://www.microscopyu.com/articles/confocal/ confocalintrobasics.html 11. Борн М., Вольф Э. Основы оптики. М: Наука, 1973. 720с. 12. БутиковН. С. Оптика: Учебное пособие для студентов физических специальностей вузов. М.: Высш. шк. 1986. 412с. 13. Бронштейн И.Н., Семендяев К.А. Справочник по математике для инженеров и учащихся втузов. М.: Физматгиз, 1962. 608 с. 14. Флетчер К. Вычислительные методы в динамике жидкостей: В 2-х томах: Т. 1: Пер. с англ. М.: Мир, 1991. 504 с.
Поступила в редколлегию 03.09.2005
Рецензент: д-р техн. наук, проф. Хаханов В.И.
Олейник Александр Игоревич, канд. техн. наук, с.н.с. лаборатории математического и компьютерного моделирования ХНУРЭ. Научные интересы: численное моделирование физико-химических процессов, уравнения математической физики. Адрес: Украина, 61166, Харьков, пр. Ленина 14. В настоящее время проходит научную стажировку (post-doctoral research) в Ecole Normale Superieure (ENS, Париж, Франция). Адрес: ENS, 24 rue Lhomond, 75231 Paris Cedex 05, France.
Дроговозов Алексей Геннадьевич, бакалавр по специальности прикладная математика, ХНУРЭ. Научные интересы: численные методы, программирование. Адрес: Украина, 61166, Харьков, пр. Ленина, 14.
Аматор Кристиан Андрэ, профессор Ecole Normale Superieure (ENS), действительный член Академии Наук Франции, декан факультета химии ENS. Научные интересы: аналитическая химия, прикладная математика, физика, физическая химия, биохимия, электроанализ, химия мозга. Address: ENS, 24 rue Lhomond, 75231 Paris Cedex 05, France.
Свирь Ирина Борисовна, д-р техн.наук, зав. лабораторией математического и компьютерного моделирования ХНУРЭ. Научные интересы: математическое моделирование в научных исследованиях: физической химии, биологии и медицине. Адрес: Украина, 61166, Харьков, пр. Ленина. 14 , е-mail: <svir@kture.kharkov.ua>.
1. Введение
Существуют системы обыкновенных дифференциальных уравнений, отдельные переменные которых имеют время ответной реакции значительно короче интервала, на котором рассматривается поведение решения этой системы. В таких системах быстро реагирующие переменные, как правило, «устойчивы» в том смысле, что если медленные переменные были в какой-то момент жестко зафиксированы, то быстро реагирующие переменные должны быстро стремиться к значениям равновесия или оставаться в окрестности этих значений.
Решение различных задач химической кинетики приводит к появлению означенных выше систем. Системы, содержащие быстро и медленно изменяющиеся компоненты, по определению являются
28
РИ, 2005, № 3
жесткими и часто решаются при помощи численных методов. Однако большинство известных численных методов решения задачи Коши для систем обыкновенных дифференциальных уравнений, таких как методы Эйлера, Рунге-Кутта, Адамса и др. [1-4], не подходят для решения жестких задач. Это обусловлено тем, что условия устойчивости этих методов накладывают слишком сильные ограничения на длину шага интегрирования. Поэтому численное решение жесткой задачи все-таки может быть получено этими методами, однако процесс решения потребовал бы неприемлемого количества шагов интегрирования и, соответственно, огромного расчетного времени. Для решения жестких задач были предложены специальные методы; это направление численного анализа активно развивается в настоящее время.
В данной работе исследуются два новых численных метода решения жестких задач Коши: квази-Рунге-Кутта (КРК) [5] и метод Алуффи-Пентини (АП) [6].
Общая идея метода КРК очень сходна с идеей метода Рунге-Кутта. Однако в отличие от обычного метода Рунге-Кутта в процессе решения задачи методом КРК также используется информация о первых и вторых производных искомого решения.
Метод Алуффи- Пентини является новым методом, разработанным профессором Римского университета (Universita di Roma “La sapienza”) Ф. Алуффи-Пентини с соавторами. Он основан на линеаризации исходной системы дифференциальных уравнений. В процессе решения линеаризованной задачи необходимо вычислять экспоненту матрицы Якоби исходной системы. В методе предложен новый подход к вычислению экспоненты матрицы, который включает автоматическую коррекцию ошибок округления.
Целью работы является анализ этих двух новых методов, сравнение результатов их работы с результатами известного метода Гира [7,8] и разработка алгоритмов для адаптивного изменения шага интегрирования для обеспечения необходимой точности численного решения жестких систем ОДУ.
2. Описание методов
2.1. Метод квази-Рунге-Кутта
Дана система обыкновенных дифференциальных уравнений: y'- f(x, у) = 0.
В начале (n- 1)-го шага численного интегрирования векторы решения аппроксимируются следующим образом:
y1n-1] ~ y(xn-1);
у2-1] - в hy,(xn-1);
y3n-1] - hV(xn-1Х
т.е., предшествующие расчеты дают приближенные значения как искомой функции у, так и ее производных первого и второго порядка.
В методе КРК, как и в классическом методе Рунге-Кутта, один шаг интегрирования разбивается на несколько промежуточных шагов. Через Yp Y2,... ,Ys обозначим промежуточные значения, аппроксимирующие y(yn_i + hq) i = 1,2,...,s , которые можно вычислить по формуле [5]:
s-1
Yi = S aijFj + Z uijy[n_1], i = 1,2
J=1
j=1
s
(1)
В конце шага вектор решения определяется по формуле
yin] = ZbijhFj + ]^У[,пЛ i = 1,2,3 , (2)
j=1 j=1
где aij, uij и bij, vij — элементы матрицы (3):
"A U“
B V _
(3)
fj = f(xn-1 + hcj,Yj) i = 1,2,...,s; cj - элемент s-размерного вектора c.
Из формул (1),(2) видно, что описанная группа методов является неявной, что является необходимым условием жесткой устойчивости метода. Элементы матрицы (3) определяются исходя из условий устойчивости метода [5]. Для метода третьего порядка коэффициенты из формул (1) и (2) составляют следующую матрицу:
2 4 2
5 0 0 1 15 "45
п 2 0 1 127 13
"144 5 720 "540
21 к 2 1 0
"20 5 5 1 20
21 к 2 1 0
"20 5 5 1 20
0 0 і 0 0 0
39 18 3 0 3 0
20 ’ 5 2 20
2.2. Метод Алуффи-Пентини
Запишем рассматриваемую систему уравнений в явной форме:
У' = f(y), (4)
с начальными условиями y(to) = Уо . Система (4) является автономной [9, 10], так как рассматриваемые гомогенные химические процессы приводят к математическим моделям в таком виде [8, 12]. Пусть y(t) — непрерывная, кусочно-гладкая аппроксимация решения, определенная ниже.
Пусть для k = 0,1,2,..., n -1 tk и 4+1 есть начальное и конечное время на шаге k +1, hk+1 = tk+1 - tk есть шаг по времени и пусть yk = y(tk). Будем интегрировать линеаризованную форму дифференциального уравнения (4) вдоль шага k +1 следующим образом.
РИ, 2005, № 3
29
Пусть J(y) — якобиан функции f(y). Для простоты примем fk = f(yk), Jk = J(yk). Тогда линеаризация f(y) в окрестности точки yk приводит к:
fk + Jk • (y -yk). (5)
Приближенное решение задачи Коши y на шаге k +1 является решением дифференциального уравнения
^ = fk + Jk • (y(t) - Yk) (6)
Исходя из определения р и ц0 , вектор р будет
представлен последним столбцом матрицы ex A , тогда как вектор \ — это первые n элементов вектора ^. Следовательно, в конце (к+1)-го шага
значение решения y(t) находится следующим образом:
yk+1 = yk +1 k(hk+iX (12)
с начальным условием y(tk) = yk . Аналитическое
решение данного уравнения на шаге k +1 имеет вид:
y(t) = т-ф(т-Jk) • fk + yk, (7)
где т = tk+i - tk и ф(в) = G_1 • (eG-I).
Начальной точкой следующего шага будет y(tk+i) = yk+i, т.е. конечная точка текущего шага.
Чтобы упростить обозначения, внутри шага k +1 положим ^(т) = J_1 • (ex J -1) • f , где I — единичная матрица размерности n х n .
Чтобы избежать аналитических и численных трудностей, когда J — особенная или плохо обусловленная матрица, мы рассмотрим расширенную задачу. Она дает возможность избежать операции обращения матрицы и, к тому же, проще и более устойчива.
Определим расширенные (n +1) — размерные векторы
Ї!
1 J
о
(8)
где n-мерный вектор lk(hk+1) состоит из первых n элементов вектора
p(hk+1) = ehk+1'A -Ро , (13)
здесь матрица а вычисляется в точке yk .
2.2.1. Вычисление экспоненты матрицы
Пусть z — квадратная матрица и I — соответствующая ей единичная матрица. Тогда по определению экспонентой матрицы z является:
eZ = lim (I +—1 .
N^-Д N J
На практике данная формула аппроксимируется более простым выражением:
e
Z
I+Z
N
N
(14)
которое в улучшенной форме принимает вид:
e
Z
1+#
b
2
(15)
и расширенную (n+1) х (n+1) матрицу:
J f
0T 0
(9)
где \ — n-мерный вектор; 0T — транспонированный n -мерный нулевой вектор. Отметим, что собственные значения матрицы а будут совпадать с собственными значениями матрицы J (за исключением добавленного нулевого значения).
В результате для расширенной задачи получим однородное дифференциальное уравнение
dp
dr
A-р
(10)
с начальными условиями р(0) = Р0 .
Точное решение данной задачи имеет вид:
р(т) = ex'A -Р0 , (11)
требует только вычисления матричной экспоненты и справедливо, даже если J — особенная матрица.
для любого натурального b . Используя выражение (15), можно вычислить экспоненту матрицы, при-
Z
меняя операцию возведения матрицы I н----в
2b
квадрат b раз. Это значительно сокращает вычислительные затраты по сравнению с прямым расчетом по формуле (14).
2.2.2. Контроль ошибок округления
Чтобы вычислить погрешности значащих цифр, которые возникают из-за округлений в процессе вычисления матричной экспоненты, применим следующий метод [6].
Идея его состоит в сохранении кумулятивного эффекта ошибок округления для введения их соответствующим образом на более поздних фазах вычислений.
Для вычисления экспоненты матрицы hA размера (n +1) х (n +1) мы определим новую матрицу размера (n + 2) х (n + 2):
M =
"h • A 0^
-t •
I s 0)
(16)
30
РИ, 2005, № 3
где о — нулевой вектор размера (n +1); s —
транспонированный вектор размера (n +1), компоненты которого задаются следующим образом:
П + 1
Sj =-h -Z alj5 j = 1,...,n +1, (17)
1=1
так что сумма элементов каждого столбца равна нулю.
Отметим, что благодаря структуре матрицы М, любая степень M содержит в себе подматрицу соответствующей степени hA, размещенную на той же позиции, где М содержит hA; следовательно,
искомую матрицу ehA можно получить из матрицы
Et
1.Z2 2 ‘ 2^ ■
Следовательно, для нормы матрицы:
Et
1 Z2 1 IIZI2
2 2b ' 2 2b
(19)
(20)
Необходимой оценкой величины b будет условие ||Et|| <є , где є — заданная точность для ||Et||. Таким образом, получим:
ґ л, 2
IIZ
b = int+ log2 2s
'v V jj
(21)
eM . Кроме того, для матрицы I л—— (которая будет
2b
возводиться в квадрат b раз для нахождения eM),
м
для всех ее квадратов, а также для e сумма элементов каждого столбца равна единице. Но из -за ошибок округления суммы вида (17) не будут точно равны единице.
Как меру ошибки округления, которая будет появляться на каждой стадии последующего процесса возведения в квадрат, определяют для каждого столбца «ошибку столбца», т.е. разность между единицей и фактическим значением суммы элементов столбца.
Чтобы вычислить кумулятивный эффект ошибок, мы распределим «ошибку столбца» между элементами столбца таким образом, чтобы их сумма снова равнялась единице.
Правило распределения «ошибки столбца» выберем следующим способом. Если мы сделаем предположение, что ошибки округления - независимые нормально распределенные случайные величины с нулевым математическим ожиданием и стандартные отклонения пропорциональны величине каждого элемента, то оптимальным правилом будет распределение ошибок между элементами столбца с весами, пропорциональными значениям квадратов данных элементов.
2.2.3. Выбор оптимального количества операций возведения в квадрат
где int+(x) — наименьшее положительное целое, большее или равное x.
Для простоты норму матрицы выберем следующим образом: ||z|| = maxZ | zkj |, и для надежности выбе-
k j *
ремb следующим образом: b = b + 3 .
3. Автоматическое управление длиной шага
Пусть, решая поставленную жесткую задачу Коши, мы исходим из некоторой точки (x0,y0). Применяя определенный метод КРК порядка p, проведем расчет для двух шагов при длине шага h. При этом получим следующие численные результаты y1 и y2 . Исходя из той же начальной точки (x0,y0), произведем расчет для одного шага длины 2h. Полученный результат обозначим w .
Представление главного члена погрешности метода на шаге найдем, исходя из следующей теоремы [1,11].
Пусть некоторым методом порядка p в результате выполнения двух шагов длины h найдено численное значение y2, а в результате выполнения одного большого шага длины 2h получено значение w . Тогда погрешность y2 может быть оценена по формуле
y(x0 + 2h) - y2 = ^1^ + Ohp+2),
2p -1 ' '
Грубую оценку подходящего значения числа b , необходимого для дальнейших вычислений, получим следующим образом [6].
Пренебрегая численной ошибкой, которая должна быть достаточно малой в нашем случае, будем рассматривать аналитическую ошибку, появляющуюся из-за ограничения значения b .
Можем показать, что матрицу (аналитической) ошибки Et:
2b
(1+Jrj = eZ • (1 + Et) (18)
можно аппроксимировать следующим образом:
y2 - w
а выражение y2 = y2 + ^p-аппроксимирует ве-
личину y(x0 + 2h) с порядком (p +1).
Следовательно, ошибку вычисляем по формуле:
error = —^ max | y2,i - wi | , (22)
2p -1 i=1,..,n
где n — количество неизвестных.
Сравнение полученной погрешности error с заданной величиной погрешности tol позволяет вычислить оптимальную длину шага:
1
h p+1 . (23)
^ error J
РИ, 2005, № 3
31
Введем гарантийный множитель (фактор) fac = 0,9 , чтобы в ближайшее время погрешность была приемлема с высокой вероятностью. Кроме того, необходимо следить, чтобы длина шага не возрастала и не убывала слишком быстро. Для этого можно выбрать новую длину шага следующим образом:
h
new
( (
h • min
facmax max
facmipfac •
і Yi
tol і p+1 error )
.(24)
V
V
JJ
Затем, если error < tol, два вычисленных шага принимаются и решение продолжается, исходя из у 2, причем в качестве новой длины шага выбирается hnew . В противном случае оба шага отбрасываются, а вычисления повторяются с длиной шага
hnew .
Максимальный коэффициент увеличения шага
facmax обычно выбирают между 1,5 и 5; facmin принимаем равным 0,01.
При работе с методом АП для выбора шага использовался тот же алгоритм, но со следующими особенностями.
По своей структуре метод АП для решения задачи Коши в векторном виде напоминает явный метод Эйлера [4] для одного обыкновенного дифференциального уравнения. Следовательно, можно предположить, что ошибка аппроксимации метода АП будет иметь такой же вид, что и ошибка аппроксимации явного метода Эйлера, т.е. s; = c;h2, где є; — ошибка i-й компоненты искомой вектор-функции, c; — константа, h — длина шага интегрирования. Таким образом, принимая порядок главного члена ошибки равным двум и используя формулы (22)-(24), получаем алгоритм контроля длины шага интегрирования для метода АП.
Знак обозначает, что реакция обратима, т.е. протекает как в прямом направлении, так и в обратном (при этом скорости протекания прямой и обратной реакции, как правило, различны). Символами А, B, C, ..., P обозначены химические вещества, участвующие в реакции Бриггса-Рошера.
Константы скоростей каждой из простейших реакций, приведенных в (25), помещены в табл. 1. Символом k| обозначена константа скорости прямой j-й реакции, k у — константа скорости обратной j-й реакции.
Концентрации Таблица 1
веществ имеют следующие начальные значения в момент времени t = 0 (для обозначения концентраций веществ используются их символические названия):
A = 0; J = 1,1;
B = 0,0001; I = 0;
C = 0,057;
K = 0,002;
D = 0,0001;
L = 0,004 E = 55; M = 0;
Номер k t k I
реакции j
1 3 1012 2,2
2 2109 -
3 1,4 • 103 -
4 45,3 -
5 1,5 • 104 1,6109
6 1,5 • 105 * * * -
7 31 -
8 104 -
9 3,2-104 -
10 410-3 91
11 9105 -
F = 0; N = 0,013;
G = 0,019; O = 0;
H = 0; P = 0. (26)
4. Математическая модель реакции Бриггса-Рошера
Химический механизм Бриггса-Рошера [13] представляется следующей схемой простейших 11 реакций:
1) A+B + C^D + E;
2) F + B + C ^ 2A;
3) G + B + 2C ^ F + A;
4) 2F ^ G + A + C;
5) G + F + С 2Н + Е;
6) I ^ J + K;
7) A + J ^ B + K + C + E;
8) H + L + E ^ F + M;
9) M + J ^ L + E + I;
10) N О;
11) D + O ^ P + E. (25)
Следует отметить, что реакция является гомогенной, т.е. вещества равномерно распределены во всем объеме и их концентрации не зависят от пространственных координат. Следовательно, будем рассматривать только изменение концентрации веществ во времени. Из начальных данных видно, что константы скоростей реакций чрезвычайно сильно различаются, поэтому соответствующая математическая модель с большой вероятностью является жесткой.
Используя алгоритмы численного моделирования гомогенных химических процессов [8] и схему реакции Бриггса-Рошера (25), получаем систему обыкновенных дифференциальных уравнений:
— = -k|ABC+ kj"DE+2k^FBC+ k^GBC2 - kyAJ+kjF2; dt
dB
— = -kf ABC + kj"DE - k£ FBC - k^ GBC2 + k^ AJ;
dt
— = -kf ABC + k 7DE - ktFBC - 2k7GBC2 + kjF2 -dt 1 12 3 4
- k^ GFC + ky H2E + ky AJ;
32
РИ, 2005, № 3
— = kf ABC - k IDE - k f, DO; dt 1 1 1 1
— = k f ABC - k Г DE + k< GFC - kj H2E + k} AJ -dt 1 15 5
- kg HLE + kc> MJ + k jb1 DO; dF
— = -ktFBC + k jGBC2 - 2k4F2 - k jGFC +
dt 2 3 45
+ kjH2E + kg HLE;
— = -kjGBC2 + kjF2 - ks GFC + kjH2E; dt 3 4 5 5
— = 2k< GFC - 2k7H2E - kg HLE; dt 5 5 g
— = -k£ I + ko MJ;
dt 6 9
— = kgI - k? AJ - ko MJ; dt 6 7
— = kg I + ky AJ; dt 6
— = -kg HLE + ko MJ;
dt g 9
— = kg HLE - kg MJ; dt g
—--kl0N + k!QO;
dt
= kfoN - kfoO - k^DO; dt
dP
= knDO. (27)
Начальными условиями для данной системы служат начальные концентрации веществ, участвующих в реакции.
Для решения приведенного уравнения применялись методы КРК и АП. Реакция исследовалась на протяжении 100 секунд от начального момента времени t = 0 .
На рис. 1,а-в показано изменение с течением времени концентраций различных веществ, участвующих в реакции Бриггса-Рошера. Результат получен с помощью метода КРК. Для того чтобы не загромождать рисунки, рассчитанные с помощью метода Алуффи-Пентини, данные наложены только на некоторые графики. Непоказанные зависимости, полученные методом АП, демонстрируют аналогичное сходство с решением, полученным методом КРК. Как видно из рис.1, два метода дают абсолютно одинаковый результат.
Время,
в
Рис.1. Рассчитанные концентрации веществ в реакции Бриггса-Рошера
5. Сравнительный анализ результатов
Графики, изображенные на рис.1, показывают, что решения системы дифференциальных уравнений (27), полученные методами квази-Рунге-Кутта и Алуффи-Пентини, визуально неразличимы. Это дает право утверждать, что оба метода сходятся к одному и тому же решению.
Поскольку система дифференциальных уравнений (27) содержит 16 неизвестных и является довольно громоздкой в плане вычислений задачей, то для её решения применялись алгоритмы адаптивного изменения шага, описанные в данной работе. Это позволяет значительно сократить время и количество вычислений (см. ниже). Применение алгоритма адаптивного изменения шага для метода КРК и метода АП приводит к изменению длины шага интегрирования, как показано на рис. 2 и 3 соответственно.
РИ, 2005, № 3
33
Рис. 2. Изменение длины шага при решении системы уравнений (27) методом квази-Рунге-Кутта
В общем случае для вычисления экспоненты матрицы требуется порядка 2bn3 + O(n2) операций,
где n — размерность матрицы, а величина b ищется по схеме, описанной в разделе 2.2.3. Для системы
(27) n = 16 , b и 69 . Следовательно, для решения задачи методом АП на каждой итерации необходимо выполнять приблизительно 549000 операций.
В случае применения метода КРК на каждом шаге интегрирования необходимо решать систему нелинейных алгебраических уравнений (1), содержащую n переменных, для чего применялся метод Ньютона. На каждой итерации метод Ньютона использует метод Гаусса для решения соответствующей линейной алгебраической системы. Как известно, для решения системы метод Гаусса требует
n3 n
— + n2 + з арифметических действий умножения
и деления [3]. На практике для решения задачи (27) на каждом шаге метод Ньютона совершал в среднем 3-4 итерации до сходимости, т.е. на каждом шаге
метод КРК выполняет 4
у
( 3 n 2 n
— + n + —
3 3
V У
операций.
Кроме того, при применении метода третьего порядка для вычисления вектора (1) требуется 30n операций. Значит, для решения (27) и на каждой итерации требуется выполнять приблизительно
4
у
( 3
n 2 n — + n2 + —
3 3
V У
+ 30n
операций. Для модели реак-
При этом необходимо заметить, что вследствие сильной жесткости рассматриваемой задачи Коши шаг интегрирования изменяется плавно не на всех участках рассмотренного отрезка временной оси. Некоторые шаги интегрирования отбрасываются в результате превышения заданного порога точности, как описано в разделе 3. В табл. 2 приведены время вычислений, количество шагов интегрирования, необходимых методам для решения системы ОДУ (27) с допустимой точностью є = 10 _1°, и количество отброшенных шагов интегрирования.
Как видно из табл. 2, время, затрачиваемое методом КРК на интегрирование, в несколько раз меньше времени вычислений методом АП. Такое различие обусловлено тем, что на каждой итерации метода АП выполняется большое количество операций, необходимых для вычисления экспоненты матрицы.
Таблица 2
Метод Время вычислений, с Количество итераций Количество отброшенных шагов
КРК 52,745 42587 42596
АП 1670,412 23734 12161
Гира 16,501 40000 163
34
ции Бриггса- Рошера это составляет прибл изитель -но 6984 операций.
Как видно из приведенных выше рассуждений, метод КРК значительно менее ресурсоемкий. Этим и обусловлена более высокая скорость его работы.
Как было указано выше, применение алгоритма адаптивного изменения шага по времени, описанного в разделе 3, позволяет существенно сократить вычислительные затраты. Однако, как видно из табл. 2, метод КРК требует большего количества шагов по времени, чем метод АП при одинаковом пороге точности методов. Это вызвано тем, что КРК является методом более высокого порядка и, соответственно, хуже работает при резких изменениях решения. В таких случаях лучшие результаты показывают методы более низкого порядка, одним из которых является метод АП (фактически явный метод Эйлера для векторных задач), которые также способны совершать более длинные шаги интегрирования. Тем не менее, в обоих случаях количество шагов, неудовлетворяющих заданному порогу точности, очень велико. Если в первом случае (метод КРК) причиной такого поведения является слишком высокий порядок метода на значительной части рассмотренного временного отрезка, то во втором случае большое количество отброшенных шагов обусловлено явным характером метода Алуф-фи-Пентини. Действительно, метод АП оказывается способным решать жесткие задачи вследствие
РИ, 2005, № 3
применения специального метода потенцирования матрицы с коррекцией ошибок округления. Однако данный метод оказывается слишком ресурсоемким на практике.
В табл. 2 также приведены данные, полученные с применением метода Гира [7], реализованного в пакете программ KinFitSim [14,15], предназначенном для моделирования любых гомогенных химических реакций и определения их параметров по экспериментальным данным. Сравнение этих данных с результатами методов КРК и АП показывает, что метод Г ира превосходит оба рассматриваемых метода по скорости вычислений, совершает меньшее количество шагов интегрирования, чем метод КРК, имеет наименьшее количество отброшенных шагов. Тем не менее, метод КРК имеет возможности для улучшения эффективности, например, за счет применения лучшей стратегии выбора длины шага, которая позволила бы уменьшить количество неудачных шагов и, возможно, использовать более длинные шаги на участках плавного изменения решения. В таком случае временные затраты метода КРК могут стать сравнимыми с методом Гира.
6. Выводы
Исследованные новые методы решения жестких систем ОДУ (метод квази-Рунге-Кутта и метод Алуффи-Пентини) были модифицированы с использованием алгоритмов для автоматического выбора длины шага интегрирования. Результаты работы данных методов сравнивались на примере модели реакции Бриггса-Рошера. Метод АП уникален по своим свойствам, так как, являясь явным, он в то же время является жестко устойчивым. Однако из -за больших вычислительных затрат его нецелесообразно применять для практических расчетов. Метод КРК достаточно прост в реализации и, хотя он превышает по временным затратам известный жестко устойчивый метод Гира, является перспективным методом решения жестких систем обыкновенных дифференциальных уравнений.
Признательность. Д-р Клименко А.В. и проф. Свирь И.Б. благодарят НАТО за грант (RIG 981488), поддерживающий проводимые исследования.
Литература: 1. Бахвалов Н.С. Численные методы. М.: Наука, 1973. 632 с. 2. Коллатц Л. Численные методы решения дифференциальных уравнений. М: Изд-во
иностр.лит., 1953. 459 с. 3. Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989. 432 с. 4. Волков E.A. Численные методы. М: Наука, 1982. 256 с. 5. Butcher J.C., Rattenbury N. Almost Runge-Kutta methods for stiff problems // J. Applied Numerical Mathematics. 2004. Vol. 53. P. 165-181. 6. Aluffi-Pentini С., De Fonzo V., Parisi V. A novel algorithm for the numerical integration of system of ordinary differential equations arising in chemical problems // Journal of Mathematical Chemistry. 2003. Vol. 33, No. 1. P. 1-15. 7. Gear C.W. Numerical Initial Value Problem in Ordinary Differential Equations. Prentice-Hall, Englewood Cliffs, NJ, 1971. 253 p. 8. КлименкоА.В., Свирь И.Б. Моделирование кинетических механизмов для фотохимического анализа // АСУ и приборы автоматики. 2002. № 121. С. 30-34. 9. Петровский ИГ. Лекции по теории обыкновенных дифференциальных уравнений. М.: Наука, 1970. 280 с. 10. Камке Э. Справочник по обыкновенным дифференциальным уравнениям. М.: Наука, 1976. 576 с. 11. Хайрер Э, Нёрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. М.: Мир, 1990. 512 с. 12. Клименко А.В. Методы определения жесткости и необходимой точности численного интегрирования моделей гомогенных химических процессов // Радиоэлектроника и информатика. 2004. № 3. С. 42-47. 13. Furrow S.D, Noyes RM. The oscillatory Briggs-Rauscher reaction. 3. A skeleton mechanism for oscillations // J. Am. Chem. Soc. 1982. Vol. 104, No. 1. P. 45-48. 14. Svirl.B, Klymenko A. V., PlatzM.S. KinFitSim — a software package to fit kinetic data // Радиоэлектроника и информатика. 2001. № 1. С. 132 — 136. 15. Svirl.B, Klymenko O. V, Platz M.S. ‘KINFITSIM’ — a software to fit kinetic data to a user selected mechanism. // Comput. Chem. 2002. Vol. 26. P. 379-386.
Поступила в редколлегию 03.09.2005
Рецензент: д-р техн. наук, проф. Хаханов В.И.
Шулык Вадим Николаевич, бакалавр по специальности прикладная математика, ХНУРЭ. Научные интересы: численные методы, программирование. Адрес: Украина, 61166, Харьков, пр. Ленина, 14.
Клименко Алексей Викторович, PhD, с.н.с. лаборатории математического и компьютерного моделирования ХНУРЭ. Научные интересы: численное моделирование физико-химических процессов, программирование, уравнения математической физики. Адрес: Украина, 61166, Харьков, пр. Ленина, 14.
Свирь Ирина Борисовна, д-р техн. наук, зав. лабораторией математического и компьютерного моделирования ХНУРЭ. Научные интересы: математическое моделирование в научных исследованиях: физической химии, биологии и медицине. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, е-mail: svir@kture.kharkov.ua.
РИ, 2005, № 3
35
