CC BY
9ISSN 2782-2001 Системы анализа и обработки данных том 88, № 4, 2022, с. 63-74
http://journals.nstu.ru/vestnik Analysis and data processing systems Vol. 88, No. 4, 2022, pp. 63-74
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ И ТЕЛЕКОММУНИКАЦИИ
INFORMATION
TECHNOLOGIES
AND TELECOMMUNICATIONS
УДК 519.876.5
DOI: 10.17212/2782-2001-2022-4-63-74
Анализ численных методов интегрирования жестких систем высокой размерности в 81ш1пТееЬ*
Ю.В. ШОРНИКОВ", К.А. ТИМОФЕЕВ4
630073, Новосибирск, пр. Карла Маркса, 20, Новосибирский государственный технический университет
а shornikov@corp.nstu.ru ь k.timofeev@corp.nstu.ru
При численном моделировании различных динамических систем часто возникает вопрос эффективности того или иного численного метода интегрирования применительно к определенному классу задач.
Важным классом задач Коши являются жесткие задачи высокой размерности. Одним из примеров такой задачи является модель проникновения помеченных радиоактивной меткой антител в пораженную опухолью ткань живого организма.
Методом прямых исходная начально-краевая задача для системы дифференциальных уравнений в частных производных аппроксимируется задачей Коши. Приведена компьютерная модель, выполненная в программном комплексе ВМиТесИ. Проведен сравнительный анализ эффективности явных адаптивных методов и диагонально-неявных методов при различных размерах шага сетки по пространственной переменной.
Показано, что наиболее эффективными при решении задач рассмотренного класса являются диагонально-неявные методы интегрирования типа Рунге-Кутты - Б1КК2 и БШК4 из пакета БМиТесИ. Метод БШКЗ уступает в связи с большим количеством вычислений матрицы Якоби. Предпочтительным в данном случаем скорее является метод БШК2, так как он имеет при практически той же производительности большее количество шагов по времени и не так сильно увеличивает шаг интегрирования при сравнительно низких настройках точности метода интегрирования. Из явных методов типа Рунге-Кутты с адаптивной численной схемой наиболее эффективным для решения задач подобного класса является метод «Адаптивный-5». Для задач подобного класса можно рекомендовать использование явных методов интегрирования «Адап-тивный-5», «Адаптивный-1» при небольшой размерности системы.
Традиционные неявные методы интегрирования Гира и Эйлера также эффективно решают данную задачу при условии эффективной реализации алгоритма вычисления матрицы Якоби.
Ключевые слова: дифференциальные уравнения, диффузия, радиоактивные метки, жесткие динамические системы, программная модель, численные методы, неявные методы интегрирования, диагонально-неявные методы интегрирования Рунге-Кутты, моделирование
* Статья получена 04 октября 2022 года.
ВВЕДЕНИЕ
При моделировании сложных динамических процессов возникает необходимость численного решения задачи Коши для систем обыкновенных дифференциальных уравнений (ОДУ). Эффективность применяемых численных методов зависит от степени жесткости и размерности задачи [1-3]. В статье проводится сравнительный анализ явных адаптивных и диагонально-неявных методов интегрирования, реализованных в отечественной инструментальной среде моделирования 8т1пТееЬ [4, 6, 12], которая предназначена для моделирования сложных динамических процессов. Программные модели в 81ш1пТееЬ специфицируются графическим языком структурных схем и символьным представлением в виде алгебро-дифференциальных уравнений средствами встроенного языка программирования РЬ. Библиотека численных методов в системе моделирования 81ш1пТееЬ состоит из традиционных и оригинальных численных схем интегрирования ОДУ.
1. ПОСТАНОВКА ЗАДАЧИ
Рассмотрим задачу проникновения помеченных радиоактивной меткой антител в пораженную опухолью ткань живого организма, сформулированную лабораторией Akzo Nobel Central Research [5, 11]. Исследование проводилось в диагностических и терапевтических целях. Процесс моделируется системой одномерных уравнений реакции-диффузии
dv ды д u ,
— = -kuv, — = —- - kuv, (1)
dt dt дх2
где к - постоянная скорости реакции; ы, v - концентрации помеченных радиоактивной меткой антител и пораженной опухолью ткани соответственно. Заданы начальные и граничные условия:
ы( х, 0) = 0, v( х, 0) = v0, х > 0; ы(0, t) = ф(0, 0 < t < tk .
После замены пространственной переменной х на С = —Х—, c > 0, аппрокси-
х + c
мации частных производных первого и второго порядков по новой пространственной переменной соответственно имеют вид:
дыj Uj+1 Uj—i
дС 2ДС
—
д ы ; ы ; —1 — 2ы ; + ы ; .I
—2- = j—j 2 j+1 , 1 < j < N.
дС2 (ДС)2
Значения Uo и UN+1 получаются из граничных условий: Щ = ф(0, uN+1 = Щ .
Т
Полагая у = (и1, VI,..., ч^, VN) и время окончания моделирования (к = 20, запишем задачу (1) в виде задачи Коши для обыкновенных дифференциальных уравнений:
% = Ж, у), (2)
у(0) = g, уе Я™, 0 < (< 20,
где N - изменяемый параметр, а компоненты функции Г((, у) определяются выражениями
У2 j+1 " У2 j -3 + р У2 j-3 - 2 У2 j-1 + У2 j+1
j (АС
/2 j =-ky2 7У2 j-1,
/2 j-1 = « j-2АС-+ р j-—2--ky2 j -1У2 j,
2АС (АС)2 (3)
где
«j = 2(jAC- 1)3/c2, Pj = (jAC-1)4/c2, 1 < j < N, j j N
У-1 (() = Ф(0, У2N+1 = У2^1, gе Ж2N, g = (0, 0, vo,..., 0, vo)T.
Функция ф(() является кусочно-непрерывной и имеет вид
[2, (е [0,5], ф(0 = \ ' 1 Ь [0, (е 5,20.
Подходящими значениями параметров являются к = 100, Vo = 1 и с = 4.
2. КОМПЬЮТЕРНАЯ МОДЕЛЬ
Структурная схема компьютерной модели (2) - (3) в SimlnTech приведена на рис. 1. Параметр phi соответствует функции ф(^), формируется блоком ступенчатого воздействия (0) и подается на вход модели (2) на встроенном языке программирования, в котором система уравнений описана в текстовом виде. Остальные блоки на схеме предназначены для выборки и вывода на график результатов моделирования.
Рис. 1. Компьютерная модель в SimInTech
Fig. 1. A computer model in SimInTech
На рис. 2 приведено текстовое описание системы уравнений в блоке PL (язык программирования) компьютерной модели SimInTech.
' = alpha[l]*(y[3] - ph^)*N/2 + beta[l]»(phi - 2»у[1] + y[3])*N*N - k»y[l]»y[2];
Рис. 2. Текстовое описание системы уравнений в PL Fig. 2. Text description of the system of equations in PL
3. ВЫЧИСЛИТЕЛЬНЫЕ ЭКСПЕРИМЕНТЫ
В настройках параметров расчета задана относительная точность решения £ = 10 3, начальный шаг интегрирования ho = 10-9. Допустимый диапазон
шага интегрирования задан от hmin = 10-20 до hmax = 1 с. Конечное время моделирования равно 20 с.
Для сравнительного анализа будем применять явные методы интегрирования типа Рунге-Кутты с адаптивным выбором численной схемы и переменным шагом «Адаптивный-1» - «Адаптивный-5» [7], а также диагонально-неявные методы типа Рунге-Кутты [10, 13-15] с переменным шагом DIRK2, DIRK3, DIRK4 [4, 8, 9], реализованные в программном комплексе SimlnTech [6].
На рис. 3-6 приведены компоненты решений У79, У172 и У199, а также компонента ф(^).
у[79]
1J3E-+0D
1.2Е+00
1.1Е+00 1Е-+00 гЕ-01 BE-01 7Е- 01 6Е-01 5Е-01 4Е-01 ЗЕ-01 2 Е-01 | 1Е-01
2.776Е-17 I_
1-
0 2 4 6 8 10 12 14 16 18 20
Время t, с
Рис. 3. Компонент решения y79 Fig. 3. Solution component У79
На рис. 7 показана зависимость времени расчета задачи от ее размерности для рассмотренных методов интегрирования.
Рис. 4. Компонент решения >>172 Fig. 4. Solution component >172
Рис. 5. Компоненты решения >19
Fig. 5. Solution components >199
Рис. 6. Функция ф(/) Fig. 6. Function ф(/)
Рис. 7. Зависимость времени расчета задачи от ее размерности Fig. 7. Dependence of the calculation time of the problem on its dimension
В таблице приведены численные характеристики методов при различных размерностях задачи, Ы - количество хороших шагов, - количество
ошибочных шагов - количество расчетов правой части ОДУ, Т - физическое время расчета в секундах.
Характеристики методов интегрирования Characteristics of integration methods
N = 50
Характеристики Адап-тив-ный-1 Адаптив-ный-2 Адаптив-ный-3 Адап-тив-ный-4 Адап-тив-ный-5 DIRK2 DIRK3 DIRK4
Nh 1187 871 2126 2565 579 321 126 108
Nbad 387 179 27 79 106 24 11 22
Nf 4334 8223 10740 7855 3321 1064 644 822
Т, с 0.26 0.44 0.61 0.51 0.18 0.11 0.23 0.09
N = 100
Nh 3994 3827 9516 10664 2151 395 144 125
Nbad 1740 1187 11 468 588 44 17 31
Nf 15467 38934 47626 32930 13109 1347 745 973
Т, с 1.71 4.05 5.06 3.81 1.37 0.29 0.91 0.26
N = 200
Nh 17608 15813 41459 45254 9255 442 165 132
Nbad 7691 5231 18 1902 2582 59 40 32
Nf 68241 163123 207369 139568 56605 1541 879 1003
Т, с 14.37 33.1 44.1 30.88 12.18 0.92 3.44 0.92
N = 400
Nh 74886 66508 172980 194507 39384 467 176 144
Nbad 33496 22126 19 7583 11243 59 32 39
Nf 291697 686948 864978 598689 241894 1611 927 1110
Т, с 122.43 281.75 360.85 269.08 99.32 2.89 14.4 3.33
N = 500
Nh 118957 105962 277235 304947 62500 486 168 145
Nbad 53552 35261 22 11798 17938 72 22 38
Nf 464050 1094525 1386265 938439 384254 1708 859 1111
Т, с 245.25 558.08 719.47 506.57 198.47 4.35 21.59 4.98
ЗАКЛЮЧЕНИЕ
Наиболее эффективными при решении задач рассмотренного класса являются диагонально-неявные методы интегрирования типа Рунге-Кутты -DIRK2 и DIRK4. Метод DIRK3 уступает в связи с большим количеством вычислений матрицы Якоби. Предпочтительным в данном случаем скорее является метод DIRK2, так как он имеет при практически той же производительности большее количество шагов по времени и не так сильно увеличивает шаг интегрирования при сравнительно низких настройках точности метода интегрирования.
Из явных методов типа Рунге-Кутты с адаптивной численной схемой наиболее эффективным для решения задач подобного класса является метод «Адаптивный-5». Для задач подобного класса можно рекомендовать использование явных методов интегрирования «Адаптивный-5», «Адаптивный-1» при небольшой размерности системы.
Полностью традиционные неявные методы интегрирования Гира и Эйлера также эффективно решают данную задачу; время расчета неявным методом Эйлера, реализованным в SimInTech, составляет 11 с против 4.98 с у DIRK4 на размерности N = 500. Таким образом, можно сделать вывод, что использование неявных методов интегрирования оправдано при условии использования алгоритмов решения СЛАУ с разреженными матрицами и эффективной реализации вычисления матрицы Якоби (в частности, выполнения LU-разложений не на каждом шаге интегрирования).
СПИСОК ЛИТЕРАТУРЫ
1. КолесовЮ.Б., СениченковЮ.Б. Моделирование систем. Динамические и гибридные системы. - СПб.: БХВ-Петербург, 2012. - 224 с.
2.Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений: жесткие и дифференциально-алгебраические задачи. - М.: Мир, 1999. - 685 с.
3. Cellier F.E., Kofman E. Continuous system simulation. - Springer US, 2006. - 644 p.
4. Скворцов Л.М. Численное решение обыкновенных дифференциальных и дифференциально-алгебраических уравнений. - М.: ДМК Пресс, 2018. - 230 с. - ISBN 978-5-97060-636-0.
5. Mazzia F., Iavernaro F. Test set for initial value problem solvers. Report 40/2003 / Department of Mathematics, University of Bari. - Bari, Italy, 2003.
6. Среда динамического моделирования технических систем SimInTech: практикум по моделированию систем автоматического регулирования / Б.А. Карташов, Е.А. Шабаев, О.С. Козлов, А.М. Щекатуров. - М.: ДМК Пресс, 2017. - 424 с. - ISBN 978-5-97060-482-3.
7. Скворцов Л.М. Явные адаптивные методы численного решения жестких систем // Математическое моделирование. - 2000. - Т. 12, № 12. - С. 97-107.
8. Скворцов Л.М. Диагонально неявные FSAL-методы Рунге-Кутты для жестких и дифференциально-алгебраических систем // Математическое моделирование. - 2002. - Т. 14, № 2. -С. 3-17.
9. Скворцов Л.М. Неявный метод пятого порядка для численного решения дифференциально-алгебраических уравнений // Журнал вычислительной математики и математической физики. - 2015. - Т. 55, № 6. - С. 978-984.
10. AlexanderR. Diagonally implicit Runge-Kutta methods for stiff O.D.E.'s // SIAM Journal on Numerical Analysis. - 1977. - Vol. l4, N 6. - P. l006-1021.
11. Lioen W.M., Swart J.J.B de. Test set for initial value problem solvers. Report MAS-R9832 / Centrum Wiskunde & Informatica (CWI). - Amsterdam, December 1998.
12. Программный комплекс для исследования динамики и проектирования технических систем / О.С. Козлов, Д.Е. Кондаков, Л.М. Скворцов, К.А. Тимофеев, В.В. Ходаковский // Информационные технологии. - 2005. - № 9. - С. 20-25.
13. Kv№rn0 A. Singly diagonally implicit Runge-Kutta methods with an explicit first stage // BIT Numerical Mathematics. - 2004. - Vol. 44, N 3. - P. 489-502.
14. HoseaM.E., Shampine L.F. Analysis and implementation of TR BDF2 // Applied Numerical Mathematics. - 1996. - Vol. 20, N 1-3. - P. 21-37.
15. Левыкин А.И., НовиковА.Е., НовиковЕ.А. (m, ¿)-схемы решения дифференциально-алгебраических и жестких систем // Сибирский журнал вычислительной математики. - 2020. -Т. 23, № 1. - С. 39-51.
Шорников Юрий Владимирович, доктор технических наук, профессор кафедры автоматизированных систем управления Новосибирского государственного технического университета, старший научный сотрудник Института автоматики и электрометрии Сибирского отделения Российской академии наук (ИАиЭ СО РАН). Основное направление научных исследований - математическое моделирование, исследования гибридных динамических систем. Имеет более 200 печатных работ. E-mail: shornikov@corp.nstu.ru
Тимофеев Константин Александрович, аспирант кафедры автоматизированных систем управления Новосибирского государственного технического университета. Основное направление научных исследований - математическое моделирование, исследования динамических систем. Имеет более 10 печатных работ. E-mail: k.timofeev@ corp.nstu.ru
Shornikov Yury V., Doctor of Technical Sciences, Professor at the department of automated control systems, Novosibirsk State Technical University, senior researcher in the Institute of Automation and Electrometry, Siberian Branch of the Russian Academy of Sciences. The main direction of his scientific research is mathematical modeling and research of hybrid dynamic systems. He has over 200 publications. E-mail: shornikov@corp.nstu.ru
Timofeev Konstantin A., a postgraduate at the department of automated control systems, Novosibirsk State Technical University. The main direction of his scientific research is mathematical modeling and studies of dynamic systems. He has more than 10 publications. E-mail: k.timofeev@corp.nstu.ru
DOI: 10.17212/2782-2001-2022-4-63-74
Analysis of numerical methods for integrating high-dimensional stiff systems in SimInTech*
Yu.V, Shornikova, K.A. Timofeevb
630073, Novosibirsk, 20 K. Marx Prospekt, Department of Automated Control Systems, Novosibirsk State Technical University
a shornikov@corp.nstu.ru b k.timofeev@corp.nstu.ru Abstract
When modeling complex dynamic processes, it becomes necessary to numerically solve the Cauchy problem for systems of ordinary differential equations (ODEs). The efficiency of the applied numerical methods depends on the degree of stiffness and dimension of the problem [1-3]. Depending on the task class, different methods behave differently.
* Received 04 October 2022.
This article provides a comparative analysis of explicit adaptive and diagonal-implicit integration methods implemented in the SimlnTech software [4]. The SimlnTech software package is designed to simulate complex dynamic processes in systems of various classes. The system supports the ability to develop models in the form of block diagrams, as well as describe systems of differential equations using the built-in programming language and simulate event-driven systems and finite automata.
It is shown that the most effective in solving problems of the considered class are the diagonal-implicit Runge-Kutta type integration methods - DIRK2 and DIRK4 from the SimlnTech package. The DIRK3 method is inferior due to the large number of calculations of the Jacobian matrix. The preferred method is rather The DIRK2 method is preferable in this case, because it has a greater number of time steps with almost the same performance and does not increase the integration step so much with relatively low accuracy settings of the integration method. Of the explicit methods of the Runge-Kutta type with an adaptive numerical scheme, the "Adaptive-5" method is the most effective for solving problems of this class. For problems of this class, we can recommend the use of explicit integration methods "Adaptive-5", "Adaptive-1" with a small system dimension.
The traditional implicit Gear and Euler integration methods also effectively solve this problem, provided that the algorithm for calculating the Jacobian matrix is effectively implemented.
Keywords: differential equations, diffusion, radioactive labels, stiff dynamic systems, software model, numerical methods, implicit integration methods, diagonal-implicit Runge-Kutta integration methods, simulation
REFERENCES
1. Kolesov Yu.B., Senichenkov Yu.B. Modelirovanie sistem. Dinamicheskie i gibridnye sistemy [Modeling systems. Dynamic and hybrid systems]. St. Petersburg, BHV-Peterburg Publ., 2012. 224 p.
2. Hairer E., Wanner G. Solving ordinary differential equations II. Stiff and differential-algebraic problems. Berlin, Springer-Verlag, 1996 ( Russ. ed.: Khairer E., Vanner G. Reshenie obyknovennykh differentsial'nykh uravnenii: zhestkie i dif-ferentsial'no-algebraicheskie zadachi. Moscow, Mir Publ., 1999. 685 p.).
3. Cellier F.E., Kofman E. Continuous system simulation. Springer US, 2006. 644 p.
4. Skvortsov L.M. Chislennoe reshenie obyknovennykh differentsial'nykh i differentsial'no-alge-braicheskikh uravnenii [Numerical solution of ordinary differential and differential-algebraic equations]. Moscow, DMK-Press, 2018. 230 p. ISBN 978-5-97060-636-0.
5. Mazzia F., Iavernaro F. Test set for initial value problem solvers. Report 40/2003. Department of Mathematics, University of Bari, Italy, 2003.
6. Kartashov B.A., Shabaev E.A., Kozlov O.S., Shchekaturov A.M. Sreda dinamicheskogo mod-elirovaniya tekhnicheskikh sistem SimInTech: praktikum po mo-delirovaniyu sistem avtomaticheskogo regulirovaniya [The environment for dynamic modeling of technical systems SimInTech: a workshop on modeling automatic control systems]. Moscow, DMK Press, 2017. 424 p. ISBN 978-5-97060-482-3.
7. Skvortsov L.M. Yavnye adaptivnye metody chislennogo resheniya zhestkikh sistem [Explicit adaptive methods for the numerical solution of stiff systems]. Matematicheskoe modelirovanie = Mathematical Modeling, 2000, vol. 12, no. 12, pp. 97-107. (In Russian).
8. Skvortsov L.M. Diagonal'no neyavnye FSAL-metody Runge-Kutty dlya zhestkikh i diffe-rentsial'no-algebraicheskikh sistem [Diagonally implicit Runge-Kutta FSAL methods for stiff and differential-algebraic systems]. Matematicheskoe modelirovanie = Mathematical Modeling, 2002, vol. 14, no. 2, pp. 3-17. (In Russian).
9. Skvortsov L.M. Neyavnyi metod pyatogo poryadka dlya chislennogo resheniya differentsi-al'no-algebraicheskikh uravnenii [An implicit fifth-order method for the numerical solution of differential-algebraic equations]. Zhurnal vychislitel'noi matematiki i matematicheskoi fiziki = Journal of Computational Mathematics and Mathematical Physics, 2015, vol. 55, no. 6, pp. 978-984. (In Russian).
10. Alexander R. Diagonally implicit Runge-Kutta methods for stiff O.D.E.'s. SIAM Journal on Numerical Analysis, 1977, vol. l4, no. 6, pp. l006-1021.
11. Lioen W.M., Swart J.J.B de. Test set for initial value problem solvers. Report MAS-R9832. Amsterdam, Centrum Wiskunde & Informatica (CWI), December 1998.
12. Kozlov O.S., Kondakov D.E., Skvortsov L.M., Timofeev K.A., Hodakovsky V.V. Programmnyi kompleks dlya issledovaniya dinamiki i proektirovaniya tekhnicheskikh sistem [Software for research dynamics and design of technical systems]. Informatsionnye tekhnologii = Information Technologies, 2005, no. 9, pp. 20-25.
13. Kvœrn0 A. Singly diagonally implicit Runge-Kutta methods with an explicit first stage. BIT Numerical Mathematics, 2004, vol. 44, no. 3, pp. 489-502.
14. Hosea M.E., Shampine L.F. Analysis and implementation of TR BDF2. Applied Numerical Mathematics, 1996, vol. 20, no. 1-3, pp. 21-37.
15. Levykin A.I., Novikov A.E., Novikov E.A. (m, k)-skhemy resheniya differentsial'no-alge-braicheskikh i zhestkikh sistem [Schemes of (m, k)-type for solving differential-algebraic and stiff systems]. Sibirskii zhurnal vychislitel'noi matematiki = Numerical Analysis and Applications, 2020, vol. 23, no.°1, pp. 39-51. (In Russian).
Для цитирования:
Шорников Ю.В., Тимофеев К.А. Анализ численных методов интегрирования жестких систем высокой размерности в SimInTech // Системы анализа и обработки данных. - 2022. -№ 4 (88). - С. 63-74. - DOI: 10.17212/2782-2001-2022-4-63-74.
For citation:
Shornikov Yu.V., Timofeev K.A. Analiz chislennykh metodov integrirovaniya zhestkikh si-stem vysokoi razmernosti v SimInTech [Analysis of numerical methods for integrating high-dimensional stiff systems in SimInTech]. Sistemy analiza i obrabotki dannykh = Analysis and Data Processing Systems, 2022, no. 4 (88), pp. 63-74. DOI: 10.17212/2782-2001-2022-4-63-74.
ISSN2782-2001, http://journals.nstu.ru/vestnik Analysis and data processing systems Vol. 88, No 4, 2022, pp. 63-74
