CC BY
219
УДК 519.714.3,517.972.6
А.В. Жариков, Е.В. Матюнин
Численное решение системы интегральных уравнений Фредгольма 2-го рода в среде Maple
A.V. Zharikov, E.V. Matyunin
Numerical Solution of Fredholm Integral Equations System of the Second Kind by using Maple Software
В работе рассматривается численное решение системы интегральных уравнений Фредгольма 2-го рода. Схема решения реализована в среде Мар1е. Проведен анализ предложенной численной схемы для различных начальных параметров и функций ядра на основе зависимости невязки от количества разбиений исходного интервала интегрирования. Ключевые слова: системы интегральных уравнений, интегральные уравнения Фредгольма, квадратурные методы.
In the paper authors consider numerical solution of Fredholm integral equations system of the second kind. Numerical solution was realized by using Maple software. The analysis of the offered numerical scheme was carried out for various initial parameters and kernel functions on the basis of dependence discrepancy for the different numbers of partitions on initial interval of integration.
Key words: integral equations systems, Fredgolm integral equations, quadrature methods.
1. Постановка задачи. В работе рассматривается система интегральных уравнений Фредгольма 2-го рода:
X (t) - ^} K (t, s) X (s)ds = f (t),
(1)
где X ( t ) =
Xi( t)
V X 2 (t )
0 Ki(t, s )Л
- искомые функции;
VK2 (t, S) 0
интегральных уравнений; f (t)
J
C1 (^ X2 ) = J C1 (^ X2 )P ( W )dW :
W
c(xi (() + ax7 (w))2 . .
^ ( ^ }} P(w)dw;
W
2r1 (w)
C2 (^ X2 ) = J C2 (^ X2 )P (w)dw ■
W
c(x? (w) + ax, (w))2 . .
Jv n \ P(w)dw,
(2)
- неко-
ядро системы
№ У V ЛО) у
торые функции.
Область О = [а, Ь] - отрезок. Решение уравнения рассматривается в пространстве С (О).
Данный вид уравнений возникает при решении задач стимулирования второго рода [1] с двумя активными элементами (АЭ) при асимметрии информированности [2]. Предполагается, что информационный вектор Н распределен на множестве ^ = Ж1Х Ж, С Я2 с плотностью Р (н) и функции затрат АЭ имеют вид:
2Г2 (Н)
где Г (н),Г2 (н) - функции квалификации активных элементов; а - некоторый параметр.
Функция дохода центра:
Н (х1з х2) = 1 ( (н ) + Х2 (н )) () Н)$н . (3)
W
При использовании центром компенсаторной системы стимулирования задача центра сводится к поиску оптимальных реализуемых действий:
[Ф(х1, x2 ) = H (xj, x2)- с, (x,, x2)- c2 (x,, x2) max,
(4)
[с ( х2) + с2 ( х2 )< Я.
В выражении (4) Я - ограничение фонда заработной платы.
При разной информированности АЭ:
І5- = 0, ^ = 0.
Эн Эн2
Как показано в работе [3], на одном из этапов решения возникает система интегральных уравнений вида:
xi (t) - j x2(s)K1 (t> S )ds = f1(tX
D
x2 (t) - j xi(s)K2 (t> S)ds = f2(tX
D
<
D
где правая часть системы зависит от положительного параметра X. Система интегральных уравнений (5) является частным случаем системы (1).
В данной работе рассматривается реализация численного решения системы интегральных уравнений Фредгольма 2-го рода (1) в программной среде Maple с помощью метода квадратур.
2. Численное интегрирование. Для аппроксимации интегралов в системе уравнений (1) воспользуемся следующей квадратурной формулой [4]:
i-1
), <6>
X1 ('l )- E— x2 (Si )K1 (j .si )= f1 ()+ e j = 1~m.
- b
(7)
A E Л2 Л ЛЛ E
(8)
(
A1 =
V
/
“ K1 ('n, Sn/
\
A2 =
K 2 ('^ Sn/2+1
/2! J \
K2 (t1,Sn )
- K 2 \'y2 , Sn
(
E Л Д
x' x2 ' -W 4
xn W
x1 /2(t1)
, < , , f2(tn),
(9)
3. Интерполирование вектора решений. Решением матричного уравнения (9) является набор значений функций Х1 (), Х2 () в узлах сетки. Обозначим Х1 =( Х1,..., Х1” ), Х2 = ( х2,..., Х2” ). Интерполирование искомых функций Х1(ґ) и Х2(ї) в виде
полиномов некоторой степени осуществляется с помощью многочлена Лагранжа:
< (t) = Ё Xl (t);
i=0
n
(t)=ё x2/| (t),
(10)
где П - количество отрезков разбиения интервала интегрирования [а, Ъ].
Учитывая формулу (6), систему интегральных уравнений (1) можно представить в следующем виде [5]:
i=0
где ха (?), ха (?) - приближенные решения системы интегральных уравнений (1) и базисные полиномы определяются как:
1 _
Х2 ()- X-------------х1 (5<- К (, ^)=Л ()+е з=1-т
1=1 п
где е - ошибка, связанная с заменой интеграла конечной суммой.
Для решения системы уравнений (7) составим матрицу коэффициентов ее левой части:
1=0,*jX -
x1 - x10 x1 - x1i-1 x1 - x1+1
x1 - x
x - xj-1 xi - xi+1
x1 - x1 x1 - x1n ’
l2 (t) = П
(11)
1=0,ІФ1 x2 x2
x2 - x2 x2- x20
x2 - x2-1 x2 - x2+1 x2 - x2-1 x2 - x2+1
x2 - x2 x2 - x2n
Vі ^ J
где E - единичная матрица соответствующей размерности,
Таким образом, система (7) запишется в матричном виде АХ=В:
Для полученного приближенного решения интегрального уравнения (1) в виде полиномов некоторой степени вычисляется невязка.
Для аппроксимации X1 (t), x2^ (t) решения Xj (t), X2(t) уравнения (1) вычисляется функция невязки [6]:
R(t) = Xa (t) - ¿JK(t, s)Xa (s)ds - f (t). (12)
D
Абсолютную величину невязки находим как:
Ra = max | Xa (t) - ¿J K (t, s)Xa (s)ds - f (t) |.
D
4. Реализация алгоритма в программной среде Maple. В среде Мар1е 8 реализован алгоритм численного решения (5), который включает следующие этапы:
1. Инициализация входных параметров. В алгоритме задается число разбиений исходного интервала, границы интегрирования, функции ядра, правая часть системы (5), интервалы для вывода графиков, параметры интерполяции.
2. Построение расчетной сетки.
3. Вычисление элементов матрицы (8). Заполняется матрица коэффициентов системы уравнений (9) на основе выражений (7).
4. Решение уравнения (9).
5. Интерполяция решений. Для решений X1 (t) и X2 (t) уравнения (9) строится интерполяционный полином заданной степени.
a
i=1
6. Вычисление невязки. По формуле (12) вычисляется функция невязки для полученных полиномов, определяющая порядок точности проводимых вычислений.
Листинг 1. Решение уравнения (5).
hfiateOfEquation := number „• (число разбиений интервала
hRateOfSystem := 2'hRateOfEqua tion; # размерность матрицы
integral Interval I := 0;
integralnterval2 := 1; if границы интервала
K1 :=K1 Iwl, w2) ; K2 :=K2 (wl ,w2) ; (ядра уравнения
right Part: =1/1 -A; if правая часть уравнения
borderOfPolynom := Array ( [-0. 01, 0.6J);
borderOfDiscrepancy := Array([-0.0002, 0.0002]);
intervalLenght: =abs (integralnterval2-abs (integrallnterval 1)) ;
h := l/hRateOfEquation;hl := 1/hRateOfSystem; (разбиения сетки
polynomStep:=l/6’hRa teOfEquation;
polynomPointl ;= Array ([polynomStep, 2 'polynomStep,
3,polynomStep, 4*polynomStep, 5*poIynomStep!) ; poIynomPoint2 := Array ([hRateOfEquation+poIynomStep, hRateCfEquation+2'polynomStep, hRateOfEquation+3'polynomStep, hfiateOfEguatj on+4’,polynomStep,
hRateOfEquation+5’polynomStep]) ;(интерполяционный многочлен nString;=hRateOfSystem;
nColumn :=hRateOfSystem; #размерность матрицы x:=array(l .. hRateOfSystem, []); x[l; :=hl; x[hRateQfEquation+l] : =h 1 ;
for i from 2 to (1/2) *nString do x[i] ;= x[i-l]+h end do;
for j from (1/2) *nColimn+2 to nColumn
do x ¡j j : =x [ j -1 j + Ji end do ;
( with (linalg); # Заполнение матрицы A
A ;= matrix(nString, nCalmn, 1);
for i from 1 io (1/2) * nStri ng do
for j from 1 to (1/2) *nColumn do if i=j then A[i, j]:=l
else A[i, j] := 0 end if end do end do;
for i from 1 to (1/2) *nString do
for j from (1/2) *nCoIШПП+1 to nColumn
do a[i, j] := simplify(subs(wl = x[i], v2 = x[j], Kl));
A[i, j] ;= -a[i, j]/ ((1/2) *nString) end do end do; for i from (1/2) *nString->-1 to r.String do
for j from 1 to (1/2)'"nColumn
do a[i, jj:=simplify(subs(wl=x[ij,w2 = x [j j, K2 '/) ;
A[i, j j := -a [ i , j ¡/( (1/2 j *n5tr ing) end do end do;
for i from (1/2) *nStzing+1 to ritring' do
for j from (1/2) *nColimn->-l to nColumn do
if i = j then A, j ] := I
A[i, jj := 0 end if end do end do;
7. Вывод результатов. Для реализации алгоритма подготовлена программа в среде вычислений Maple 8. Фрагмент решения системы интегральных уравнений представлен на листинге 1.
= proc (i, j) options operator, arrow;
= matrix (nString, 1, f);
= linsoIve( А, В); ( Решение матричного уравнение AX=B
5. Расчетный пример использования описанного алгоритма. Рассмотрим пример использования данного алгоритма для системы уравнений (5) с ядрами:
К^,м2) =------- —-----+ С08(^ + м2); (13)
К2( ^ »2) =
---( М?1 + ^2 +^1-^2)
2 - е 2
1
-------V-------------7'
--(»і+»2 +»1-»2)
4 - е 2
1
( »1+ »2 + Му»2)
-• (14)
Уравнение определено на интервале [0,1], а правая часть /(») = —1—, где X = 5 .
1 — X
На рисунках 1, 2 представлены графики численных решений Х1(») и Х2(») системы уравнений (7) для различного числа разбиений - т исходного интервала.
Рис. 1. Численное решение Х[(») и Х2(№) при т = 30
Рис. 2. Численное решение Х[(№) и х2(») при т = 60
Х?(4
■#»
Рис. 3. График интерполяционных многочленов а а
х1 (»), Х2 (») , при т = 30
Для восстановления искомых функций X (м2) , Х2 (м1) - по полученным точечным решениям - используется интерполяционный многочлен 4-й степени.
Получаем следующие интерполяционные полиномы х? (/), х2? (/):
х^ (»2)=0.0094»24 + 0.0241»/ + 0.0952»22 + 0.2166» + 0.0057 Х2 (»)=0.0103»! + 0.0262»3 + 0.1029»2 + 0.2327» + 0.0224. ,
(15)
На рисунке 3 представлен график интерполяционных полиномов для 30 разбиений исходного интервала.
Интерполяционные многочлены для 60 и 120 разбиений исходного интервала, формулы (15) и (16)-(17) соответственно:
XX») = 0.0095»/ + 0.0240»3 + 0.0953»2 + 0.2166» + 0.0056, х£(») = 0.0104»4 + 0.0260» + 0.1030»2 + 0.2326» + 0.0223,
Xа (») = 0.0096м!,4 + 0.0239»3 + 0.0953»2 + 0.2166» + 0.0056, X (») = 0.0104»4 + 0.0259»3 + 0.1031»2 + 0.2326» + 0.0223.
(16)
(17)
Далее вычисляется функция невязки решений. На рисунках 4-6 представлены графики функций невязок для различного числа разбиений исходного интервала (т = 30, т = 60, т = 120).
б
Рис. 4. а - график функции невязки для Х1 (м) при т = 30; б - график функции невязки для Х2 (м) при т = 30
а
Рис. 5. а - график функции невязки для Х1 (w) при т = 60; б - график функции невязки для Х2 (V) при т = 60
б
а
Рис. 6. а - график функции невязки для х1 (V) при т = 120; б - график функции невязки для х2а (V) при т = 120
Анализ приведенных зависимостей показывает, что при изменении числа разбиений исходного интервала от 30 до 120 абсолютная величина невязки уменьшается для х^ (м>) от 0.000071 до 0.000016,
для х2 (•»>) - от 0.00014 до 0.000041. Таким образом,
полученные результаты предлагаемого метода решения систем интегральных уравнений Фредгольма 2-го рода демонстрируют, что представленный метод применим для решения изложенных в работе задач.
б
а
Библиографический список
1. Новиков Д. А. Теория управления организационными системами. - М., 2005.
2. Жариков А.В. Модели стимулирования агентов промышленной корпорации в условиях асимметрии информированности // Известия АлтГУ. - 2010. - №1.
3. Векуа Н.П. Системы сингулярных интегральных уравнений. - М., 1989.
4. Волков Е.А. Численные методы. - М., 2004.
5. Васильева А.Б., Тихонов Н.А. Интегральные уравнения. - М., 2002.
6. Манжиров А.В., Полянин А.Д. Методы решения интегральных уравнений. - М., 2008.
