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

Вычислительные технологии

Том 7, № 2, 2002

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

Н. Г. Блндурин Волгоградская государственная архитектурно-строительная

академия, Россия e-mail: bandurin_n@bayard.ru

The new arbitrary high order numerical method for the solution of general integro-differential equations based on the interpolation procedure is suggested. This procedure is not adapted to the function but to it derivative.

Введение

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

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

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

© Н.Г. Бандурин, 2002.

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

1. Вывод основных соотношений

Для получения формулы интегрирования предположим, что на отрезке [а, Ь] вещественной оси х задана сетка, в N узлах которой известны значения производной к-го порядка функции у = f (х), определенной на [а,Ь], и значения производных низших порядков в начальном узле х = х^ Ставится задача вычисления в этих же узлах значений функции, ее производных и интегралов, выраженных через узловые значения производной порядка к и производные при х = XI.

Интерполируя на отрезке [х^хп] (п < N) производную у(к)(х) многочленом степени п — 1, после соответствующих выкладок можно получить соотношение [2]

п п

у(й)(х) = ЕЕ ^ 1у?)х<-1- (1)

г=1 j=1

Здесь определителем матрицы Z является определитель Вандермонда.

Значения производной порядка (к — 1) в узлах интерполяции получаются в результате последовательного интегрирования (1) в пределах отрезков [х1,х^] (г = 1, 2, ...,п):

у,('"') = уГ> + зу у<к) (2)

(в — квадратная матрица интегрирования порядка п). Если N > п, то, применяя последовательно формулу (2) на отрезках [ху(п_1)+1, ху(п_1)+п] = 0,1, 2,...), можно получить компоненты вектора производной Y(k_1) во всех узлах отрезка [а, Ь]:

у(к_^ = 1у(й-1) + £У(к), (3)

где I — единичный вектор, Б — квадратная матрица интегрирования порядка N.

Выполнив интегрирование в соответствии с (3) р раз, получим формулы для вычисления производных в узлах интерполяции, выраженных через значения производных порядка к в этих же узлах и производные низших порядков в начальном узле х = х1 :

р_1

у(к_р) = ^ Бт1у(к_р+т) + БрУ(к) при р < к. (4)

т=0

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

Если р > к, то компонентами вектора Y(k_p) в (4) будут значения интегралов, а формула примет вид

р_1

Y(к_р) у ^ БР_й+т1у(т) + Б(5)

т=р_к

В (5) учтено, что у[к р) = 0 при р > к, поскольку в соответствии с принятыми обозначениями эти величины являются значениями интегралов в узле х = Х\.

Формулу дифференцирования можно получить, полагая к = 1 в равенстве (3) и решая его относительно вектора производных первого порядка V. В результате можно записать

У = (Е - Б Б )1у'1 + ОТ, (6)

где Е и Б — единичная матрица и матрица дифференцирования порядка N.

Выполнив дифференцирование в соответствии с (6) р раз, получим формулу для вычисления производных р-го порядка в узлах интерполяции, выраженных через производные в начальной точке х = Х1 и значения функции в этих узлах:

р- 1

у(р) = ^ Бт(Е - ББ)1 у(р-т) + БРУ. (7)

т=О

Формула (7) применима для решения как обыкновенных ИДУ, так и уравнений с частными производными.

Для любого целого р > 0 имеют место равенства:

N

= (х - Х1)р/р!, (8)

3 = 1

N

= 0 (^ =1, 2,..., N), (9)

3=1

которые выражают результаты интегрирования и дифференцирования функции f (х) = 1.

Очевидно, что предлагаемая интерполяционная процедура применима не только для функций, но и для их приращений. Использование формул (4) и (7) при численном решении краевых задач для систем интегродифференциальных уравнений дает возможность:

1) с высокой точностью аппроксимировать системы ИДУ, имеющие достаточно гладкие решения, их дискретными аналогами, так как формулы являются точными на отрезке [а, Ь] для всех многочленов степени не выше п;

2) решать уравнения с заданными значениями производных не только на краях, но и внутри области интегрирования, поскольку при выводе формул имеется возможность за начальную точку на отрезке [а, Ь] принимать любой узел Хг (г = 1, 2,..., N);

3) при решении ИДУ высокого порядка исключить необходимость выполнения алгебраических преобразований на краях области, так как полученные формулы не зависят от номеров узлов, расположенных вне области интегрирования, и содержат в явном виде необходимые для постановки граничных условий производные.

Используя формулу Тейлора, можно получить оценку погрешности аппроксимации функции и ее производных 0(Ка+1), где к — шаг сетки. При вычислении смешанных производных при решении т-мерных задач дифференцирование по формуле (5) выполняется

последовательно по всем аргументам, а оценка погрешности (при постоянном п) будет т

равна 0( = кп+1\.

2. Алгоритм численной реализации

Ниже приводится вывод формул итерационного процесса для вычисления узловых значений функции у = f (х), являющейся решением нелинейного ИДУ г-го порядка на отрезке [а,Ь], который имеет в своем составе интеграл I(х) с постоянными пределами интегрирования х1 = а и х^ = Ь и интеграл 3(х) с переменным верхним пределом:

f [х,у(х),у'(х), ...,у(г)(х), I(х), 3(х)] = 0, (10)

где

1 (х) = I Ф,у(х),у'(х),---,у{г)(х),з,у(з>),у')з>),---,у{'г)(з>)¥з;

3(х) = J ф[х,у(х),у'(х),--,у(г)(х),в,у(в),у'(в),--,у(г)(в)]^в.

а

Краевые условия, зависящие только от а и Ь, могут быть записаны в виде

/¿[а, Ь, у(а), у'(а),..., у(г_1)(а), у(Ь), у'(Ь),..., у(г_1)(Ь)] = 0, г = 1, 2,...,г. (11)

Последовательно подставляя в (10) х1 = а, х2,..., х^ = Ь, можно получить N уравнений относительно узловых значений функции у и ее производных. Эти уравнения и краевые условия (11) можно записать в виде одной системы ^ + г) уравнений с N (г + 1) неизвестными узловыми значениями функции и ее производных:

Е(У, У',..., У(г)) = 0, (12)

где (У(у))т = (у(у),у(у),...,у}У)), ^ = 0,1, 2,...,г.

В результате дифференцирования (12) можно получить приближенное равенство

г

ДЕ ДУ(у), (13)

у=о

где элементами матрицы являются производные

дУ = + + т =12 N (14)

д*т = дуУ + д1дуУ) + д3дуУ), т =1,(14)

При решении нелинейных уравнений производные в (14) вычислялись по простейшей формуле

^ = 0.5№ + Да) — ^ (а — Да)]/Да. (15)

Ниже предполагается, что после т-й итерации известны узловые значения интегралов 1т и 3т, вектор Ут и вектор невязки Ет. Интегралы вычисляются по формуле (3) при к = 1. Если с использованием (7) выразить в (13) векторы ДУУ) через вектор ДУ, элементы которого должны быть линейно независимыми величинами приращений узловых значений функции и ее производных в фиксированных точках, то для вычисления ДУ на (т + 1)-й итерации можно записать систему линейных алгебраических уравнений

СтДУт+1 + Ет = 0. (16)

Ь

а

х

При определении набора основных неизвестных в ДУ и числа уравнений в (16) должны быть выполнены соответствующие преобразования для того, чтобы матрица О была квадратной и невырожденной. Например, если в рассматриваемом случае одного ИДУ предположить, что на краях задано значение функции и все (г — 1) производных, то корректно сформированный вектор ДУ должен иметь ^ + г — 1) компонентов: N значений функции в узлах и (г — 1) производных на краях. В то же время определенное выше число уравнений равно ^ + г). Так как на одном из краев значение функции задано, то уравнение, получаемое в результате подстановки в (10) соответствующего этому краю значения аргумента, оказывается лишним и должно быть удалено из системы уравнений (16). Такие преобразования имеют формальный характер и могут быть запрограммированы.

Вектор узловых значений функции и ее производных в фиксированных точках после (т + 1)-й итерации вычисляется суммированием

Ут+1 Ут + ДУт+1-

Остальные операции итерационного процесса являются стандартными и здесь не приводятся. При решении систем уравнений и многомерных задач изложенный алгоритм естественным образом обобщается. Заметим, что решение систем линейных ИДУ получается в результате выполнения одной итерации.

3. Тестовые примеры

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

Пример 1. Краевая задача для системы нелинейных интегродифференциальных уравнений с неизвестными функциями V(x,y,z) и W(x,y,z) и параметрами a и b. Решение следует найти внутри параллелепипеда R = {0 < x,y < п/2,0 < z < 1}:

п/2 п/2 1 x y z

Vxxx+Vyy+V J jy aVxx(a, e, y) Wz(a, prf)dadpdj+J J J Vxxx(a, Prf)Wx(a, prf)dadpdj+ 0 0 0 0 0 0 + exp(z) cos y(cos x — sinx(1 — e)) + exp(a) sinx sin y(exp(z) — 1) = 0,

n/2 1

Wx — Wz + J Jv(x,e,Y)dpdf — sinx(e — 1) = 0, 00

b2 + b + exp(a — 1) — 7 = 0.

Краевые условия:

1) на краях x = 0, x = n/2: V¡(0, y, z) — exp(z) cos y = 0;

n/2 1

V(0,y,z) + y У Vx(0, в, y) Wz(0,e,Y)dpdf + a — e = 0, Vxx(n/2, y, z) + exp(z) cos y = 0, 00

W(0, y, z) — z = 0;

2) на краях y = 0, y = п/2: V(x, n/2, z) = 0,

x z

V (*, 0, г)+Л V (c,n/2,7ÍWx(a,n/2, f)dad7 — exp(z) sin x — (exp(z) — 1«cos x — 4 = 0;

00

3) на краю z =1: W(x, y, 1) + VXx(x, y, 1) = 0.

Точное решение V = sin x cos y exp(z), W = x + z, a = 1, b =2.

Начальное приближение VO = x, Wo = 2, ao = 3, bo = 4.

В табл. 1 приведены результаты решения на ПЭВМ Pentium (433 МГц) задачи (8), (9) при различных размерах сетки узлов. Степень интерполяционного полинома n принималась равной числу узлов сетки N для каждой координатной оси. Относительная е и максимальная абсолютная А погрешности вычислялись по формулам

£

(Fijk - Fjfc)2/N

1/2

А = max |Fjfc - Fjfc| , i, j, k = 1, N,

где Fjfc и Fjk — вычисленные и точные значения функции в узле (i,j, k). Итерационный процесс завершался при уменьшении невязки в уравнениях (8), (9) до значения Iter = 1.0e — 11. Установлено, что число итераций главным образом зависит от степени близости начального приближения к точному решению, значительно увеличивается при уменьшении Iter и практически не зависит от размеров сетки. Отметим, что приведенное точное решение задачи (8), (9) не является единственным и при некоторых других вариантах начального приближения итерационный процесс может сходиться к другому решению.

Таблица 1

Размеры сетки £y £W Ay Aw Число итераций Время счета, с

3 x 3 x 3 1.4e-2 4.5e-3 5.8e-2 8.6e-3 8 8

5 x 5 x 5 6.4e-4 7.7e-5 4.8e-3 2.1e-4 8 193

7 x 7 x 7 9.0e-6 2.7e-6 7.6e-5 8.2e-6 8 3193

Пример 2. Рассматривается начально-краевая задача для существенно нелинейного интегродифференциального уравнения. Его решение необходимо найти в области R = {0 < ж < 1, 0 < t < tk}:

ж t ж t

Ухх - ytt + y2 - У + J j Ужж(а, в+ J yxx(a,t)ytt(«,t)da+y yxt(x,e-0 0 0 0

-ж3(ж - 2/3) sin21 - 2(x + 1)sint - 0.5x2(1 - cos2x) = 0. Краевые условия при x = 0, x =1

y(0,t) + у(1,в)yt(1,в- 0.25(1 - cos 2t) = 0,

y(1,t) - sint = 0.

Начальные условия при t = 0

y(x, 0) + yxt(x, 0) - 2x = 0, yt(x, 0) - x2 = 0.

3

е

t

o

Точное решение y = x2 sin t.

Для получения решения отрезок [0, tk] разбивается на малые отрезки длиной Ht, на которые в свою очередь наносится сетка узлов с шагом ht. На ось x наносится сетка с шагом hx. Решение задачи (10) - (12) на отрезке [0, tfc] получается в результате последовательного решения уравнений (10), (11) на каждом малом отрезке оси t, причем начальные условия для решения на последующем отрезке определяются из решения на предыдущем. Для tfc = 20 результаты решения задачи представлены в табл. 2. Степень интерполяционного многочлена n принималась равной числу узлов сетки при решении системы (10), (11). Так как число итераций, необходимое для уменьшения невязки до значения Iter, на каждом этапе решения различно, то в таблице представлены максимальные значения этого числа. Число узлов на оси x принято минимальным, равным 3, и его увеличение не повышает точность решения, так как точное решение зависит от координаты x лишь во второй степени. Из сравнения первых двух строк таблицы видно, что для экономии машинного времени значение Iter должно приниматься в зависимости от требуемой точности решения, которая в свою очередь определяется степенью интерполяционного многочлена n.

Таблица 2

Размеры сетки hx ht Iter A Число итераций Время счета, с

3 х 3 0.5 0.2 1.0e—5 1.5e—3 26 23

3 х 3 0.5 0.2 1.0e—13 1.5e—3 68 55

3 х 5 0.5 0.1 1.0e—11 2.0e—5 67 102

3 х 7 0.5 1/240 1.0e—13 9.6e—13 71 1520

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

Увв

у1 + = 1

Здесь аргументом является длина дуги упругой линии в; второе уравнение определяет неизменность полной длины стержня в процессе деформирования. Если ось стержня в недеформированном состоянии совпадает с осью х, то будем иметь краевые условия при в = 0

х(0) = 0, хв(0) = 1, у(0) = 0, ув(0) = 0.

Так как в соответствии с (13) кривизна упругой линии стержня не зависит от в, то для превращения стержня в кольцо кратности к необходимо приложить момент М = 2кп. В табл. 3 представлены результаты решения этой задачи, где пх — число узлов сетки; Хк,Ук — координаты конца стержня после деформирования ("дефекты" кольца); п — число выполненных итераций; ¿с — время решения.

4. Программы для решения ИДУ

Для решения корректных нежестких задач на основе интерполяционной формулы (7) разработана и отлажена содержательная часть следующих программ.

Таблица 3

Пх к Хк Ук П ¿с

11 1 —5.6е-4 7.1е—5 23 1

19 1 8.1е-7 6.6е—8 23 5

29 1 5.5е—8 1.3е—9 23 15

19 2 9.7е—3 7.0е—3 42 10

29 2 2.2е—4 3.2е—5 47 32

1. Программа для решения краевых задач для систем обыкновенных ИДУ, которые могут содержать подлежащие вычислению числовые параметры и соответствующее число дополнительных уравнений.

2. Программа для решения задачи Коши для систем обыкновенных ИДУ.

3. Программа для решения двумерных краевых задач для систем ИДУ, которые могут содержать подлежащие вычислению числовые параметры.

4. Программа для решения начально-краевых задач для ИДУ, зависящих от аргументов X и ¿.

5. Программа для решения четырехмерных краевых задач для систем ИДУ порядка не выше второго.

6. Программа для решения пятимерных краевых задач для систем ИДУ порядка не выше второго.

7. Программа для решения двумерных краевых задач для систем ИДУ порядка не выше второго при произвольной области интегрирования, на краях которой заданы значения искомых функций.

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

Список литературы

[1] Бахвалов Н. С. Численные методы. М.: Наука, 1973.

[2] КОРН Г., КОРН Т. Справочник по математике. М.: Наука, 1968.

[3] Попов Е. П. Теория и расчет гибких упругих стержней. М.: Наука, 1986.

Поступила в редакцию 26 марта 2001 г., в переработанном виде —16 октября 2001 г.