CC BY
23
УДК 519.62
МЕТОД И ПРОГРАММА ДЛЯ РЕШЕНИЯ СИСТЕМ ОБЫКНОВЕННЫХ НЕЛИНЕЙНЫХ ИНТЕГРО-АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ ОБЩЕГО ВИДА
© 2011 г. Н.Г. Бандурин
Волгоградский государственный архитектурно- Volgograd State University
строительный университет, of Architecture and Civil Engineering,
ул. Академическая, 1, г. Волгоград, 400074, Akademicheskaya St., 1, Volgograd, 400074,
info@vgasu.ru info@vgasu.ru
Кратко описывается численный метод, алгоритм и программа для автоматического решения систем корректных нелинейных интегро-алгебраическихуравнений (ИАУ) достаточно произвольной структуры. Под автоматическим решением понимается получение результата без выполнения этапов выбора метода, программирования и отладки программы. Приводятся результаты решения 3 тестовых примеров, которые представлены и решены в 2 видах: как системы ИАУ и дифференциальных уравнений, получаемые из ИАУ в результате их дифференцирования. При этом процедура подготовки исходных данных для старта программы и степень вычислительной сложности решения остаются практически прежними.
Ключевые слова: нелинейные интегро-алгебраические уравнения общего вида, численные методы, компьютерные программы.
The numerical method, algorithm and program for automatic solution of the systems of correct nonlinear integro -algebraic equations (IAE) which have a quite random structure are described briefly. Under the automatic solution the author understands getting result without executing the stages of selecting a method, of programming and program checking. Results of solution of 3 test examples are cited.
Keywords: nonlinear ordinary integro-algebraic equations, numerical methods, computer Programs.
Постановка задачи Метод решения
Для решения системы (1) принят метод шагов [1], в соответствии с которым весь промежуток интегрирования [0,b] должен быть представлен в виде Q больших шагов длиной т, т.е. принимается b = Qt . Для численной реализации метода каждый большой шаг разбивается на m узлов, и весь вычислительный процесс состоит из последовательного выполнения Q этапов. В результате выполнения каждого этапа получаются значения искомых функций в m узлах соответствующего большого шага. При выполнении последующего этапа в качестве начальных условий принимаются значения функций, вычисленных в конце предыдущего. Для приведения систем ИАУ и ДУ к системе обычных нелинейных алгебраических уравнений на каждом этапе решения в алгоритме используются формулы интегрирования и дифференцирования соответственно, опубликованные впервые в [2-5]. Ниже приводится в кратком виде вывод этих формул применительно к поставленной задаче.
Ниже кратко описывается метод, алгоритм и программа для решения системы нелинейных интегро-алгебраических уравнений (ИАУ) достаточно общего вида, т.е. начальной задачи Коши
РМ, МО, >2(0,..., Ум С)] + (1)
Я I
ач | [u, >11(и) У2(u),..., Ум (и)¥и = 1 0
Здесь Ы,Я - число уравнений и интегралов; а^ -числовые параметры; 1 < г < М .
Очевидно, что в результате дифференцирования (1) по / получается система обыкновенных в общем случае существенно нелинейных дифференциальных уравнений (ДУ).
Предполагается, что задан полный набор начальных условий при /=0.
Требуется с использованием разработанной программы найти решение системы (1) на равномерной сетке ¿1 = 0, ¿2,-.., ¿м = Ь отрезка [0,Ь] в 2 формах: как решение системы ИАУ и соответствующей системы ДУ
Для получения формулы интегрирования предположим, что на отрезке [с,^ вещественной оси t задана равномерная сетка, в P узлах которой известны значения производной ^го порядка функции у=/(/), определенной на [с^, и значения производных низших порядков в начальном узле ^ = с. Ставится задача вычисления в этих же узлах значений функции, ее производных и интегралов, выраженных через узловые значения производной порядка k и производные при ^ = с.
Производную у№)(0 на отрезке [/ь4] (п < P) интерполируем многочленом степени (п-1) [2, 5]
y(k)(t)=tTc, (2)
где ^(М,.../"-1); ст=(сьс2,...,с„).
Полагая в (2) / = /1, /2, ...,/„, можно получить матричное соотношение у®=2с, где y(k)- вектор, компонентами которого являются узловые значения производной порядка k функции /(/), а элементы матрицы Z - хц= 1.
Определителем матрицы Z является определитель Вандермонда, поэтому c=Z-1 y(k) и (2) примет вид у(\/)=^-1 y(k).
(3)
Интегрируя последнее равенство последовательно в пределах отрезков [/ь/2],..., [/ь/„], получим
"(Ы)(/)=! yl(k-1)+BZ-1 y(k). (4)
У
Здесь у(А-1) - вектор, компонентами которого являются узловые значения производной ^-1)-го порядка; I - вектор размером п, элементы которого равны единице; элементы матрицы В - Ь^ = (// - ^)/ ] .
Из (4) получаем формулу интегрирования на отрезке [/ь/п]
у (к-1) = 1у(к-1) + (к), (5)
где s=BZ-1- квадратная матрица интегрирования порядка п.
Если Р>„, то, применяя последовательно (5) на отрезках [//(„-О+ь/"-^"] (/=0, 1, 2,.), можно получить компоненты вектора производной у(к-1) во всех узлах отрезка [с,^
у(к-1) = 1у(к-1) + §у(к), (6)
где 8 - квадратная матрица интегрирования порядка Р.
Например, при с - d = 3, „ = 3, Р = 4 эта матрица имеет вид
S =
0
0,37500 0,33333 0,37500
0
0,79167 1,33333 1,12500
0
- 0,20833 0,33333 1,12500
0
0,41667 0
0,37500
Y
(k-r) _
= Z smiy[
m=0
my. .(k-r+m)
+ Sr Y(k) при r < k.
(7)
Если г'^, то компонентами вектора У№ г) в (7) будут значения интегралов, а формула примет вид
Y(k-r) = Z Sr -k+mIy[m) + Sr Y
m=0
(k)
Для получения формулы дифференцирования матрицу интегрирования представим в виде [5]
§ = 8п , где = о , элементы матрицы-строки
21 822J
я12 равны нулю, а квадратная матрица я22 порядка (Р-1) - не особенная.
Отсюда можно получить выражение вектора узловых значений производной функции у(/) через начальные величины производной и узловые значения собственно функции
Y' = (E - DS )Iy'+ DY,
(8)
где D =
0
s-11
s 2212
0
s -1
s
- матрица дифференцирования.
При c - d = 3, n = 3, P = 4 эта матрица имеет вид 0 0 0 0
D
- 0,9444 0,5000 - 0,5000 - 0,0555 0,7778 - 2,0000 1,0000 0,2222
- 2,1667 4,5000 - 4,5000 2,1667
Выполнив дифференцирование в соответствии с (8) г раз, получим формулу для вычисления производных г-го порядка в узлах интерполяции, выраженных через производные в начальной точке / = /1 и значения функции в этих узлах
r-1
Y(r) = Z DmE - DSIy[r-m) + DrY .
(r-m)
(9)
Выполнив интегрирование в соответствии с (6) г раз, получим формулы для вычисления производных в узлах интерполяции, выраженные через производные порядка k в этих же узлах и производные низших порядков в начальном узле х=х1
В последнем равенстве учтено, что у(к г) = 0 при
г>к, поскольку в соответствии с принятыми обозначениями эти величины - значения интегралов в узле / = /1.
Для любого целого p > 0 имеют место равенства
p p
PS? = (tt -t.)p /p!, ZDp+1 = 0, (i = 1,2,...,P), кото-
j=i j=i
рые выражают результаты интегрирования и дифференцирования функции fx) = 1.
Использование формул (7) и (9) [2, 5] при численном решении краевых и начально-краевых задач для систем интегро-дифференциальных алгебраических уравнений (ИДАУ) в общем случае дает возможность: 1) с высокой точностью аппроксимировать эти системы, имеющие достаточно гладкие решения, их дискретными аналогами, так как формулы являются точными на отрезке [c,d] для всех многочленов степени не выше n; 2) решать уравнения с заданными значениями производных не только на краях, но и внутри области интегрирования, поскольку при выводе формул имеется возможность за начальную точку на отрезке [c,d принимать любой узел ti (i = 1,2,.,P); 3) при решении ИДАУ высокого порядка с использованием интерполирования многочленами исключить необходимость выполнения алгебраических преобразований на краях области, как это рекомендуется, например, в [6], так как полученные формулы (6) и (9) не зависят от номеров узлов, расположенных вне области интегрирования, и содержат в явном виде необходимые для постановки граничных условий производные. Используя формулу Тейлора, можно получить оценку погрешности аппроксимации функции и ее производных O(hn+1), где h - шаг сетки.
Для старта программы, написанной на алгоритмическом языке Delphi, должны быть подготовлены исходные данные (тексты уравнений, начальные условия и числовые параметры задачи) в виде последова-
m=0
r-1
r-1
тельности строк, содержание которых определяется подсказками для каждой строки. По мере ввода исходных данных и в процессе вычислений выполняется проверка на наличие ошибок формального характера, и при их обнаружении выдается соответствующее сообщение. Одна из особенностей метода и программы - имеется возможность выводить на печать и сохранять графики y = ft), y\ = f (y.) yt = f (y j), в то
время как известные методы и программы для решения задачи Коши основаны на вычислении только узловых значений искомых функций [7-10].
Демонстрационные версии программ для решения ИАУ и других 1-, 2- и 3-мерных уравнений доступны для скачивания и использования в некоммерческих целях [11].
Тестовые примеры
Все приведенные ниже системы ИАУ и соответствующие ДУ решаются с помощью разработанной программы в течение нескольких минут, включая сюда и время для подготовки исходных данных. Итерационный процесс начинался из нулевого приближения.
Пример 1. Задача Коши для системы 3 ИАУ с неизвестными функциями u(t) v(t) и w(t) (b = 10):
t
u + Ju(s)ds + w3 - 2sinw + costsin21 -1 = 0,
0
t
v-Jv(s)ds + ev -esmt+cost = 0, (10)
0
t
w + J u(s)ds + sin w - sin(cost) -1 = 0,
0
u(0) = 0, v(0) -1 = 0, w(0) -1 = 0.
Точное аналитическое решение: u (t) = sin t, v(t) = sin t+cos t, w(t) = cos t.
Получены значения максимального отклонения выполненного численного решения от точного (Д) для различных малых шагов сетки h: Д^=002 = 1,4-10-3; Дй=0,01 = 9,110 4; Дh=0,05 = 8,7-10-5; Дh=o,o2 = 5,5-10-4; Дh=o,ol = 9,0-10-4.
В результате дифференцирования (10) по t получаем систему существенно нелинейных ДУ:
u + u + 3w2w - 2cost - sint + 3 cos21sint = 0,
V - v + evV - esin t+cos t (cos t - sin t) = 0,
w'+ u + cos(w)w' + cos(cos t) sin t = 0,
в результате решения которой максимальные отклонения Д^0,2 = 4,0-10-4; Д^од = 6,8-10-4; Дм,05 = 6,5^10-4; Д^0,02 = 5,5^10-4; Д^0,01 = 4,6^10-4. Точность решения одной и той же задачи в различных формах представления приблизительно одинакова, что является достаточно строгим доказательством надежности полученных результатов расчета.
Пример 2. Система 3 ИАУ (b = 10):
u + w3 - e' sint + (sint + cos t)3 = 0,
t
v + Jw(s)ds - e~' cos t - 0,5 sin2t = 0,
0
t
w + Jw(s)ds + sinv - 0,5(sin2t + sin21) - sine-t cos t = 0 ,
0
u(0) = 0, v(0) -1 = 0, w(0) = 0. Точное аналитическое решение: u(t) = elsin t, v(t) = e- sin t, w(t) = sin t cos t.
Максимальные отклонения: Д^0,2 = 1,8-10"7;
Ah=0,1 = 1,110 9; Дh=o,o5 = 1/M0-12; Дh=o,o2 = 1,M0-12; Дh=0,05 = 1,9^10-12.
Максимальные отклонения для соответствующей системы ДУ (которая здесь не приводится):
Дh=o,2 = 1,3 10 4; Ah=o,1 = 1,110 7; Дh=o,o5 = 6,540-10; Дh=o,o2 = 9,2-10-1°; Дh=o,ol = 7,640-10. Пример 3. Система 3 ИАУ (b = 10):
t
u + v - sin w + J (u(s)v(s))ds - sht - cht - e-t - sin(sin( t)) -1 = 0,
0
t
v - u + Ju(s)ds - 2sht + cht = 0, w + u - sint + cht = 0,
0
u(0)-1 = 0, v(0)-1 = 0 , w(0) = 0 . Точное аналитическое решение: u(t) = el sin t, v(t) = e-cos t, w(t) = sin t cos t.
Максимальные отклонения: Д^0,02 = 1,6-10-6; Дни = 7,440-10; Дh=o,o5 = 2,7-10-1°; Дhlo,o2 = 7,6^10-10; Дh=o,ol = 5,0-10-1°.
Итерационный процесс для соответствующей системы ДУ оказался крайне неустойчивым. Удалось получить по сути дела лишь одно решение для b = 4 и h = 0,05 с максимальным отклонением Д = 6,3 •Ю-12.
Все приведенные выше результаты решения уравнений были получены при шаге численного дифференцирования, равном 10-12, что приблизительно согласуется с точностью арифметических операций программного комплекса. При уменьшении или увеличении этого шага на два порядка точность решения изменяется незначительно.
В заключение следует подчеркнуть, что созданная компьютерная программа позволяет в автоматическом режиме получать решение систем нелинейных интегро-алгебраических уравнений, недоступных для известных зарубежных программных комплексов.
Литература
1. Эльсгольц Л.Э., Норкин С.Б. Введение в теорию дифференциальных уравнений с отклоняющимся аргументом. М., 1971. 296 с.
2. Бандурин Н.Г. Новый численный метод порядка n для решения интегро-дифференциальных уравнений общего вида // Вычислительные технологии. 2002. № 2. С. 3 - 10.
3. Бандурин Н.Г., Игнатьев В.А. Комплекс программ для решения систем нелинейных интегро-дифференциальных уравнений (ИДУ), содержащих производные и интегралы дробного порядка // Интеллектуальные системы и компьютерные науки : материалы IX междунар. конф. 23 - 27 октября 2006 г. Т. 2, ч. 1. М., 2006. С. 56-58.
4. Бандурин Н.Г., Игнатьев В.А. Пакет программ для решения систем нелинейных интегро-дифференциальных уравнений (одно-, двух- и трехмерные начально-краевые задачи) // Мат. моделирование. 2007. Т. 19, № 2. С. 105 - 111.
5. Бандурин Н.Г. Численное решение нелинейных ин-тегро-дифференциальных уравнений с запаздывающим аргументом // Вычисл. технологии. 2010. № 3. С. 31 - 38.
6. Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. М., 1968. 720 с.
7. Холл Д., Уатт Д. Современные численные методы решения обыкновенных дифференциальных уравнений. М., 1979. 312 с.
8. Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. М., 1990. 512 с.
Поступила в редакцию_
9. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М., 1999. 685 с.
10. Gear C. W. The automatic integration of stiff ordinary differential equations // Information Processing 68. Proceedings of the IFIR Congress. Edinburg, 1968. Р. 187 - 193.
11. URL: http://www.vgasu.ru/rus/science/inn/a-driff.php (дата обращения: 20.06.2010).
_29 ноября 2010 г.
