CC BY
49
Вычислительные технологии
Том 19, № 5, 2014
Метод и пакет программ для численного решения систем существенно нелинейных интегро-дифференциально-алгебраических уравнений (корректные по Адамару двумерные и трёхмерные краевые задачи)
Н.Г. Блндурин1, С.Ю. Калашников1
Кратко описываются метод, алгоритм и пакет программ для решения корректных по Адамару двумерных и трёхмерных задач для систем существенно нелинейных (т. е. не разрешённых относительно старших производных) интегро-дифферен-циально-алгебраических уравнений. Для приведения последних к системе обычных нелинейных алгебраических уравнений в алгоритме используются полученные ранее формулы интегрирования и дифференцирования. В качестве неизвестных система уравнений может также содержать числовые параметры. Приводятся результаты решения четырёх тестовых примеров.
Ключевые слова: существенно нелинейные интегро-дифференциальные уравнения, интегро-дифференциально-алгебраические уравнения, численные методы, компьютерные программы.
1. Введение и постановка задачи
Необходимость численного решения краевых задач для систем нелинейных дифференциальных, интегральных и интегро-дифференциально-алгебраических уравнений (ИДАУ), зависящих от двух и более аргументов, часто возникает при математическом моделировании различных физических процессов в науке и технике [1-7]. Наиболее используемым в настоящее время методом для численного решения таких задач является метод конечных разностей [8-11], на основе которого решаются прикладные задачи математической физики [12, 13]. В связи с развитием вычислительной техники этот метод успешно используется на многопроцессорных системах [14, 15]. Несмотря на то что фундаментальные основы для решения уравнений с частными производными хорошо известны [16-24], их программная реализация связана со значительными трудностями, поскольку не существует алгоритмов для решения таких уравнений в общем виде, что заставляет исследователей при решении каждой конкретной прикладной задачи в качестве начального шага составлять соответствующую программу и выполнять её отладку. Широко используемые программные комплексы (Maple, Macsyma, Mathematica, Reduce, Matlab и др.), предлагаемые на рынке программных продуктов известными зарубежными корпорациями, не могут быть применены для решения систем существенно
1 Волгоградский государственный архитектурно-строительный университет, Россия Контактный e-mail: bandurin_ng@mail.ru
нелинейных интегро-дифференциально-алгебраических уравнений (СНИДАУ) с частными производными. Ниже кратко описывается предложенный ранее численный метод [25-28] и приводятся сведения о программах для автоматического решения на ПК корректных по Адамару двумерных и трёхмерных систем, малых по размеру, в общем случае существенно нелинейных интегро-дифференциально-алгебраических уравнений [29], содержащих неизвестные числовые параметры, которые могут входить не только в состав основных уравнений, но и в краевые условия. Области интегрирования рассматриваемых уравнений должны быть ограничены координатными линиями или поверхностями в общем случае криволинейных неортогональных систем координат. Решаются следующие задачи:
— краевые задачи для двумерных систем ИДАУ;
— краевые задачи для трёхмерных систем ИДАУ.
2. Метод решения и программная реализация
Для преобразования интегро-дифференциально-алгебраической задачи к системе нелинейных алгебраических уравнений относительно узловых значений функций, которые затем должны быть решены известными методами, используется формула дифференцирования для функций W, которая в одномерном случае имеет вид [25, 28]
г-1
W(r) = ^ Бг(Е - ББ)№— + ВгW,
г=0
где W(r), W, I — вектор узловых значений производных функций порядка г, вектор узловых значений функций и единичный вектор соответственно; Б, Б, Е — матрицы интегрирования, дифференцирования и единичная матрица соответственно; и\-1 — значения функций и их производных в начальном узле х = х1.
Векторы узловых значений смешанных производных получаются в результате применения формулы дифференцирования к приведенному равенству в другом направлении. Поскольку эта операция является формальной, то повышение порядка производных в системах корректных ИДАУ не влечет за собой появления дополнительных трудностей в программной реализации алгоритма. Следует заметить, что элементы матриц интегрирования и дифференцирования не могут быть записаны в виде формул, используемых для формального программирования, так как представляют собой интегралы и программируются на основе их вычисления [25, 28].
Разработанные для решения систем ИДАУ программы работают на основе минимума исходных данных, т. е. необходимость писать код программы для решения каждой задачи и выполнять её отладку исключается. Для старта итерационного процесса в память ПК следует ввести следующие данные: число уравнений, тексты решаемых уравнений в обычной математической форме, параметры сетки узлов, максимальные значения порядка производных, точное решение задачи (если оно известно) для проверки адекватности вычислений, начальное приближение и краевые условия. Перед заполнением очередной строки исходных данных на дисплей выводится подсказка о содержимом этой строки. После ввода строки выполняется проверка правильности данного содержимого, а в случае обнаружения ошибки ввод прекращается и на экран выводится соответствующее сообщение.
В отличие от всех известных программных комплексов в нашем случае программы вычисляют и выводят на печать узловые значения не только функций, но также и их производных, в том числе, возможно, смешанных.
3. Пакет программ
Пакет состоит из следующих двух отдельных программ.
1. Программа БОиКВБ2 предназначена для численного решения корректных двумерных краевых задач в системе координат (ж, у) для существенно нелинейных систем ИДАУ, которые в своём составе кроме подлежащих вычислению функций могут иметь также неизвестные числовые параметры.
Для успешного численного решения подобных задач принимаются обычные условия и ограничения:
— система уравнений должна быть корректной и не жёсткой;
— число основных ИДАУ должно равняться числу неизвестных функций;
— если система уравнений содержит некоторое число подлежащих вычислению числовых параметров, то к основным уравнениям необходимо добавить соответствующее число дополнительных уравнений, содержащих эти параметры;
— каждое основное и каждое дополнительное уравнение может содержать в своём составе производные любой функции, интегралы и искомые параметры, входящие в подынтегральные выражения интегралов. Предусмотрено использование 10 типов интегралов: два типа с постоянными и переменными верхними пределами и областью интегрирования, совпадающей с областью, на которой ищется решение задачи, остальные восемь типов с постоянными и переменными верхними пределами при интегрировании по четырём координатным линиям, ограничивающим область (эти последние могут быть включены в соответствующие краевые условия);
— для каждой неизвестной функции основные уравнения должны иметь в своём составе по крайней мере одно уравнение с несмешанными производными этой функции наивысшего порядка;
— для каждого неизвестного параметра дополнительные уравнения должны иметь в своём составе уравнение, содержащее данный параметр;
— в области интегрирования решение задачи не должно быть сильно осциллирующим;
— размеры решаемых задач ограничиваются ресурсами ПК.
В области интегрирования предусмотрены два способа формирования сетки узлов, при первом из которых число узлов следует задавать для каждой из осей отдельно, при втором — шаг сетки принимается одинаковым для обеих осей.
2. Программа БОиКВБЗ используется для численного решения существенно нелинейных ИДАУ в системе координат (х,у, г). Условия и ограничения для успешного решения этих уравнений в основном те же, что и для программы БОиКВБ2, но с соответствующими поправками. Например, предусмотрено использование 14 типов интегралов: два типа с постоянными и переменными верхними пределами и областью интегрирования, совпадающей с областью, на которой ищется решение задачи, остальные 12 типов интегралов с постоянными и переменными верхними пределами при интегрировании по шести координатным поверхностям, ограничивающим область.
4. Тестовые примеры
4.1. Краевая задача для трёх нелинейных ИДАУ, содержащая два числовых параметра (решение выполняется по программе BOUNDS2)
Uxxxx + u3 + uxxx sin[2(x + y)] + u cos[a(x + y)] - sin(x + y) = 0, vx + vy + uw — sin(x + y) cos(xy3) — x — y = 0,
n n/2 x y
w + J j vy (s,t)dsdt + J j xsv(p, — cos(xy3) — x3y3/6 — n3/16 = 0,
0 0 0 0
x y
a + 4b — b2 + J j xsv(p, — x3y3 — 3 = 0,
00
b — a + w — cos(xy3) = 0.
Границы области интегрирования 0 < x, y < n/2. Краевые условия
Uxxx (0, y) + u(0, y) + v(0, y) + [v(0, y)]2 — sin y + cos y = 0,
Ux(0,y) + v(n/2,y) — cos y — ny/2 = 0,
uxx(0,y) + sin y = 0; n/2
v(0, y) = 0, uxx(n/2, y) + У aw(0, s)ds/3 + sin(n/a + y) — y = 0,
0
n/2
v(x, 0) + v(x, 0) + x J [a/3 + u(s,n/2)]ds — (1 + n/2)x = 0. 0
Точное решение u(x, y) = sin(x + y), v(x, y) = xy, w(x, y) = cos(xy3), a = 3, b = 4. Начальное приближение u0(x,y) =v0(x,y) =w0(x,y) = 0, a0 = 2,1, b0 = 3, 2. В результате решения на сетках nxm получены следующие значения максимального отклонения для функций приближённого решения от точного А:
Д5Х5 = 0.021, А13Х13 = 1.6 ■ 10-8, А7Х13 = 1.5 ■ 10-4, Ai3x7 = 2.5 ■ 10-8.
4.2. Краевая задача для трёх существенно нелинейных ИДАУ, содержащая два числовых параметра (BOUNDS2)
uxx + uyy + w3 + 2 sin(x + y) — (sin x cos y)3 = 0,
x y
vx + vy + J j u(s,t)dsdt + u — sin x — sin y — 1 = 0, 00
n/2 n/2 x y x y
w + J j v(s,í)dsdí^ J u(s,t)dsdt + J j vx(s,t)dsdt — sinx cos y — xy = 0, 0 0 0 0 0 0
п/2 р/2
a + b3 + J У v(p, s)dpds - 64 = 0, о о
b + sin a — sin 3 — 4 = 0.
Границы области интегрирования 0 < x, y < п/2. Краевые условия
u(0,y) — a sin y/3 = 0, u(n/2,y) — cos y = 0, v(0, y) + y + a — b +1 = 0, u(x, 0) — sin x = 0,
n/2
м(х,п/2)^У uy(s,n/2)ds — cosx + sinb — sin4+1 = 0, о
v(x, 0) — x = 0.
Точное решение u(x, y) = sin(x + y), v(x, y) = x — y, w(x, y) = sin x cos y, a = 3, b = 4. Начальное приближение u0(x,y) =v0(x,y) =w0(x,y) = 0, a0 = 2.1, b0 = 3.2. В результате решения на сетках nxm получены следующие значения максимального отклонения приближённого решения от точного А:
Д5Х5 = 8.1 ■ 10-4, А15Х15 = 2.0 ■ 10-10, А5Х15 = 4.6 ■ 10-4, А15Х5 = 7.1 ■ 10-4. 4.3. Краевая задача для трёх существенно нелинейных ИДАУ (BOUNDS3)
x y z п/2 п/2 п/2
uxx + uyy + uzz + u + J J J u(p, s, t)vy(p, s, t)dpdsdt + J J J u(p, s, t)dpdsdp—
0 0 0 0 0 0 — cos(x + z) — cos(y + z) + cos(x + y + z) — cos x — cos y — cos z — cos(x + y) — 1 = 0,
n/2 n/2 n/2
vx + vy + v + j J J u(p, s, t)dpdsdt — x — y + z — 4 = 0, 000
x y z
w + vx + J J J[u(p, s,t) — uxy (p, s, t)]dpdsdt — sin x cos y cos z — 1 = 0. 000
Границы области интегрирования 0 < x, y,z < n/2. Краевые условия
тг/2тг/2
u(0,y,z) + ^ J ux(n/2, s, t)dsdt — sin(y + z) + 2 = 0, 00
u(n/2, y, z) — cos(y + z) = 0, v(0, y, z) — y + z = 0,
u(x, 0, z) — sin(x + z) = 0, u(x, n/2, z) — cos(x + z) = 0, v(x, 0, z) — x + z = 0, u(x, y, 0) — sin(x + y) = 0,
тг/2тг/2
u(x,y,n/2) + y J ux(p, s, 0)dpds — cos(x + y) = 0. 00
Точное решение u(x, y, z) = sin(x + y + z), v(x, y, z) = x + y — z, w(x, y, z) = sinx cos yx cos z.
Начальное приближение u0(x,y,z) = 0, v0(x,y,z) = 0, w0(x,y,z) = 0. В результате решения на сетках 5 x 5 x 5, 7 x 7 x 7 и 9 x 9 x 9 получены соответственно следующие значения максимального отклонения приближённого решения от точного А:
3.1 ■ 10-2, 4.7 ■ 10-4, 10-6.
4.4. Краевая задача для двух существенно нелинейных ИДАУ (BOUNDS3)
uxx + uyy + uzz + u + sin(uyy) + v + 2 sin(x + y — z) — sin(sin(z — x — y)) — x + y — z = 0,
v3 + v — exp(ux) + sin
x y z
///u(M,r)<MS<fr
000
+ exp(cos(x + y — z)) — sin(xyz(x — y + z)/2) —
—(x — y + z)3 — x + y + z = 0.
Границы области интегрирования 0 < x,y, z < n/2. Краевые условия
y z
u(0,y,z) + y J u(0,s,t)dsdt + sin y — sin z — 2sin(y — z) = 0, 00
u(n/2, y, z) — cos(y — z) = 0,
u(x, 0, z) — sin(x — z) = 0, u(x, n/2, z) — cos(x — z) = 0,
(тг/2 п/2 \
/ / u(s, t, 0)dsdt I — sin(x + y) — 1 = 0, 00
x y
u(x,y,n/2) + y J u(s, t, n/2)dsdt + cos x + cos y — 1 = 0. 00
Точное решение u(x, y, z) = sin(x + y — z), v(x, y, z) = x — y + z. Начальное приближение u0(x, y, z) = v0(x, y, z) = 0.
В результате решения максимальное отклонение приближённого решения от точного А получено для второй функции v. Для разных сеток значения этого отклонения были равны:
А(7х7Х7) = 1.7■ 10-4, А(11х11х11) = 1.0■ 10-8, А(7х7х11) = 2.0■ 10-4, А(11х7х7) = 5.6■ 10-5.
В заключение заметим, что, по мнению авторов, в математической литературе не содержится сведений о методах численного решения многомерных систем существенно нелинейных интегро-дифференциально-алгебраических уравнений.
Список литературы
[1] Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: Наука, 1977. 535 с.
[2] Пирумов У.Г., Гидаспов В.Ю., Даниелян А.А. и др. Численный анализ двухфазного течения в газодинамическом фильтре // Матем. моделирование. 1998. Т. 10, № 11. С. 19-28.
[3] Головизнин В.М., Коршия Т.К., Самарский А.А. и др. Вариационные схемы магнитной гидродинамики в произвольной системе координат // Вычисл. математика и матем. физика. 1981. Т. 21, № 1. С. 54-68.
[4] Иванова К.Н., Четверушкин Б.Н., Чурбанова Н.Г. Квазигазодинамическая система уравнений и уравнения Навье — Стокса // Матем. моделирование. 2004. Т. 16, № 4. С. 98-104.
[5] Ковеня В.М. Об одном алгоритме решения уравнений Навье — Стокса вязкой несжимаемой жидкости // Вычисл. технологии. 2006. Т. 11, № 2. С. 39-51.
[6] Ичетовкин Д.А., Паасонен В.И. Численное исследование независимой аппроксимации граничных условий на решениях с разрывами производных // Там же. 2010. Т. 15, № 1. С. 77-84.
[7] Шайдуров А.А., Щепановская Г.И., Якубович М.В. Численное моделирование течений вязкого теплопроводного газа в канале // Там же. 2013. Т. 18, № 4. С. 77-90.
[8] Марков А.А. Исчисление конечных разностей. 2-е изд. Одесса: Mathess, 1910.
[9] Гельфонд А.О. Исчисление конечных разностей. М.: Наука, 1967.
[10] Годунов С.К., Рябенький В.С. Введение в теорию разностных схем. М.: Наука, 1973. 400 с.
[11] Волчинская М.И., Гольдин В.Я., Калиткин Н.Н. Сравнительное исследование разностных схем для уравнений акустики // Вычисл. математика и матем. физика. 1974. Т. 14, № 4. С. 919-927.
[12] Шокин Ю.И., Федотова З.И. О достижениях в теории разностных схем // Вычисл. технологии. 1999. Т. 4, № 5. С. 56-69.
[13] Ладонкина М.Е., Милюкова О.Ю., Тишкин В.Ф. Численный метод решения уравнений диффузионного типа на основе использования многосеточного метода // Вычисл. математика и матем. физика. 2010. Т. 50, № 8. С. 1438-1461.
[14] Четверушкин Б.Н. Кинетические схемы и высокопроизводительные многопроцессорные вычисления в газовой динамике // Вычисл. технологии. 2002. Т. 7, № 6. С. 65-89.
[15] Чеванин В.С. Численное моделирование развития гидродинамических неустойчивостей на многопроцессорных системах // Матем. моделирование. 2012. Т. 24, № 2. С. 17-32.
[16] Конторович Л.В., Крылов В.И. Приближённые методы высшего анализа. М.: Гостех-издат, 1949.
[17] Yanenko N.N. Methode a pas Fractionnaires: Resolutions de Problemes Polydimensionnels de Physique Mathematique. Paris: Librairie Armand Colin, 1968. 205 p.
[18] Марчук Г.И. Методы вычислительной математики. М.: Наука, 1980.
[19] Морс Ф.М., Фешбах Г. Методы теоретической физики: Пер. с англ. М.: Мир, 1958.
[20] Березин И.С., Жидков Н.П. Методы вычислений. Т. 2. М.: Физматгиз, 1962. 639 с.
[21] Курант Р. Уравнения с частными производными: Пер. с англ. М.: Мир, 1964. 830 с.
[22] Бахвалов Н.С. Численные методы. Т. 1. М.: Наука, 1973. 631 с.
[23] Калиткин Н.Н. Численные методы. М.: Наука, 1978. 512 с.
[24] На Ц. Вычислительные методы решения прикладных граничных задач. М.: Мир, 1982. 294 с.
[25] Бандурин Н.Г. Новый численный метод порядка n для решения интегро-дифферен-циальных уравнений общего вида // Вычисл. технологии. 2002. Т. 7, № 2. С. 3-10.
[26] Бандурин Н.Г. Большие динамические перемещения и деформации массивной растяжимой нити в движущей среде // Прикладная математика и механика. 2003. Т. 67, № 2. С. 266-271.
[27] Бандурин Н.Г. Численное решение существенно нелинейных интегро-дифференциаль-ных уравнений с запаздывающим аргументом // Вычисл. технологии. 2010. Т. 15, № 3. С. 31-38.
[28] Бандурин Н.Г., Гуреева Н.А. Метод и пакет программ для численного решения систем существенно нелинейных обыкновенных интегро-дифференциально-алгебраических уравнений // Матем. моделирование. 2012. Т. 24, № 2. С. 3-16.
[29] Математическая энциклопедия. Т. 3. М.: Советская энциклопедия, 1982. 951 с.
Поступила в 'редакцию 18 февраля 2014 г.
Method and program package for solution of systems of essentially nonlinear integro-differential-algebraic equations (two- and three-dimensional boundary value problems)
Bandurin Nikolay G.1, Kalashnikov Sergey Yu.1
Goal: the development of method, algorithm and software package for numerical solution of the well-posed (in the sense of Hadamard) two- and three-dimensional problems for the systems of highly nonlinear (i. e., in general, not solved with respect to higher derivatives) integro-differential algebraic equations.
Method: contrary to the known methods, based on the procedure of interpolation of the unknown function by algebraic polynomials, an absolutely original approach is applied, which is based on the procedure of interpolation of the derivative of the unknown function. As a result, the necessity is eliminated to use the grid nodes, which do not belong to the domain, where the integration is performed. Consequently, the process of solving highly nonlinear integro-differential algebraic equations can be maximally formalized: it is sufficient to enter the equations in the usual mathematical form and necessary numerical
1 Volgograd State University of Architecture and Civil Engineering, 400074, Volgograd, Russia Corresponding author: Bandurin Nikolay G., e-mail: bandurin_ng@mail.ru
parameters into the computer memory. When solving the partial differential equations, the described procedure is used in the first direction first and then the interpolation procedure in the second direction is performed using the obtained results, etc. The results of the solution of test problems for highly nonlinear integro-differential algebraic equations, which are obtained using the software in automatic mode, i. e. using the minimum of initial data without additional programming, confirm the high accuracy of the suggested approach.
Conclusion: the first time in the history of computational mathematics the method is suggested and the first test programs are developed for the automatic solution of the well-posed (in the sense of Hadamard) systems of highly nonlinear integro-differential algebraic equations. Contrary to the all known software packages, the developed software computes and outtypes the node values not only of functions, but also their derivatives and, among them, probably, mixed ones.
Keywords: highly nonlinear integro-differential algebraic equations, numerical methods, computer programs.
Received 18 February 2014
