УДК 517.9
И. В. Бойков, Ю. Ф. Захарова
ОБ ОДНОМ ПРИБЛИЖЕННОМ МЕТОДЕ ИДЕНТИФИКАЦИИ СИСТЕМ С РАСПРЕДЕЛЕННЫМИ ПАРАМЕТРАМИ
Предложен приближенный метод идентификации параметров динамических систем, описываемых линейными параболическими и гиперболическими уравнениями. Идентификация заключается в определении функции Грина или коэффициентов уравнения. Метод основан на сведении задачи к уравнению в свертках. Показано, что этот метод применим для восстановления начальных условий при известном выходном сигнале и функции Грина. Приведены численные примеры.
Введение
Идентификация динамических систем с распределенными параметрами является некорректной задачей и относится к классу обратных задач. Различные методы нахождения коэффициентов систем параболических и гиперболических уравнений, а также обширная библиография приведены в [1, 2].
Однако методы, предложенные во многих из этих работ, трудно применить при решении практических задач. Поэтому возникает необходимость в разработке и программной реализации более простых, устойчивых и быстро сходящихся алгоритмов.
В работе [3] предложен метод сведения уравнений в частных производных к интегральным уравнениям и их последующая аппроксимация системами алгебраических уравнений высокой размерности. В связи с некорректностью задачи, данные системы требуют специальных методов решения.
Ниже приводится достаточно простой и эффективный метод идентификации динамических систем, описываемых линейными параболическими и гиперболическими уравнениями. Для определенности рассуждения проводятся на примере уравнений параболического типа. При этом данный метод является общим для всех линейных уравнений в частных производных, которые на основе интегральных преобразований или в результате использования функции Грина можно представить в виде интегральных уравнений в свертках.
1 Идентификация параметров динамических систем с распределенными параметрами
Известно [4], что задача Коши для параболического уравнения
Т , , du \ . . du ^^7 / \ du . . . .
L(u) = -— У akj(t) + У bk(0— + c(t)u = 0, (1)
dt j dxk dxj k=1 dxk
lim u (t, x) = uo ( x) (2)
t ^0
и уравнение
L(u ) = f (t, x) (3)
при нулевых начальных и граничных условиях интегральным преобразованием Фурье сводятся к интегральным уравнениям в свертках.
В частности задача Коши (1), (2) сводится к интегральному уравнению
и ((’ х) =-[[°((’х -^)“0(^)ё(4)
(2п)п/2 О
где интеграл берется по области О, которая является носителем функции щ (х).
Полагая функцию ^ (х) = 0 при х е Еп \ О, уравнение (4) можно рассматривать в пространстве Е2:
и((’ х) =-^пиДО°((’ х -^)“0(^)4^ (5)
Еп
где I е [0,~), х е Еп .
В случае, если рассматривается уравнение (3) при нулевых начальных и граничных значениях, то его решение представимо в виде
и ((’ х) =-(П+1)/2 | 4т | ОЦ - т, х - £)/(т, 1)ёI,
(2П) — Еп
где
ГО(I, х), I > 0
0 ('-х) = ' 0, ,< 0.
Ниже исследуется несколько задач идентификации параметров динамических систем, описываемых параболическими уравнениями.
Задача 1. Требуется, располагая решением задачи Коши (1), (2) или решением уравнения (3) при заданной функции /(^, х), найти функцию О(¿, х).
При приближенном вычислении функции О(¿, х) по решению задачи Коши (1), (2) будем считать известными значения ^(х) на сетке, состоящей из Ып узлов, а также значения функции и (¿, х) на прямоугольной сетке, состоящей из Ып узлов по переменной х, - Аг- < х1 <... < хN < А^ 1 = 1, п на каждом временном слое 0 < ¿1 < ¿2 < ■■■ < М < Т .
При вычислении функции О(¿, х) по решению уравнения (3) будем считать известными значения функций и (¿, х) и /(¿, х) в Ып узлах по переменной х на М временных слоях 0 < ¿1 < ¿2 < ■■■< М < Т . Значения А^,
1 = 1, п и Т определяются физической постановкой задачи. Ниже для определенности будем считать, что Аг- = А при всех значениях 1 и для простоты
обозначений положим п = 2 .
Построим алгоритм восстановления функции О (^, х) в задаче Коши (1), (2).
Зафиксируем ^ = ¿1 и применим к уравнению (5) преобразование Фурье по пространственным переменным. В результате приходим к уравнению в частотной области:
и (¿1,о>1,Ю2) = О (¿1,юью2)и0(0^2). (6)
Отметим, что прямое и обратное преобразование Фурье при этом будем вычислять по формулам
Р(и(х1,х2)) = и(01,Ю2) = | | и(х1,х2)е1(Ю1х'+Ю2х2)
—^ ^
<( х1; х2) = р _1(и (ю1,ю2)) = — Г Г и (ю1,ю2)е“1(ю'х1+°2 х2) й ю14ю2
2п
Для вычисления функций и(¿1,о>1,Ю2) и и0 (о>1, О) воспользуемся ку-батурными формулами преобразований Фурье.
Пусть функция щ (х1, х2) задана на сетке узлов [Мк }, к = 1, N , распо-
2
ложенных в области Б = [ - А, А] . Если задана прямоугольная сетка [Ук, V/},
к, / = 0, N -1, Ук =-А + 2АN, к = 0, N -1, то каждому узлу {Ук/}, Ук1 = V, V/) ставится в соответствие прямоугольник А к/ = У, Ук+1; V/, V/+1].
к, / = 0, N -1, и прямое преобразование Фурье вычисляется по формуле
N-1 N-1
и0( '
, N-1 N-1
,(01,02) = £ £ и(Ук,у1 )Це1(юЛ +°2х2)йх1йх2 . (7)
2п - -
^ к=0 /=0 Ак/
2
Если сетка [Мк}, к = 1, N неравномерна, то область Б покрывается
квадратами Ак/, к, / = 0, N -1, введенными выше. Каждому такому квадрату
ставится в соответствие число и* (Ак/), которое определяется следующим образом:
1) если в квадрате А к/ содержится узел сетки [МУ}, то полагаем и* (Ак/) = и (Му );
2) если таких узлов несколько, то можно использовать любой из них (формула будет более точна, если взять узел, ближайший к точке ^к, V),
_ (Ук+У к+1) V ^ =------2----);
3) если в квадрате Ак/ нет узлов сетки [МУ}, то полагаем
и * (Ак/) = и (Мм,), где М^ - ближайший к А к/ узел сетки [Мк } .
Тогда преобразование Фурье вычисляется по формуле
, N-1 N-1
и0(01,Ю2) = — £ £ и*(Ак/)Цег(°Л +Ю2*2)йх1йх2 . (8)
к=0 /=0 аи
Значение функции и0(Ю1,Ю2) вычисляется на сетке узлов (юк ,ю2),
к, / = 0, N -1, где ок = - В + , к = 0, N, 1 = 1, 2. Здесь В - достаточно
большое положительное число, В > А .
Аналогичным образом вычисляются и значения функции и(^а^,^)
на сетке узлов (¿1 ,юк ,ю2 ), к, / = 0, N -1. 10
Приравнивая левую и правую части уравнения (6) на указанной сетке, 2
приходим к системе из N алгебраических уравнений вида
и (¿1,»к ,о2) = 0(^1,юк ,»2)и0К ,02), к, / = 0, N-1. (9)
К каждому из уравнений системы (9) применим итерационный метод
х
Оп+1 ('1 >4 >ю2) = апО ('1 >0 >»2) + (1 -аИ) х ^ ('1,Юк,»2) -Ук,1 {Оп ('1>™к,ю2)и0(юк,ю2) - и(Л’О,ю2)}
(10)
где 0 <а* <ап <а <1, а параметры у к / подбираются таким образом, чтобы выполнялись условия
у к ,іио(юк >ю2)=-•
(11)
Сходимость метода следует из теоремы Обломской [5]. Предельные значения последовательных приближений обозначим через О* (¿1,0 ,“2) . Если в
точке (юк,о2) значения и0(юк,“2) = 0 , то уравнение (9) не рассматривается, а
* к /
значение О (¿1,00 ,“2) в квадрате А к/ полагается равным среднему значению
* к /
О (¿1,00 ,02) по всем квадратам, имеющим общие грани с А к/.
Так как все коэффициенты Ук / удовлетворяют условиям (11), то последовательность приближений (10) сходится со скоростью геометрической прогрессии со знаменателем, равным 1/2:
О* (¿1,0 ,о2) - Оп+1(^1,шк ,о2)
<1 -2
адЧ ,ю2)ио(®к ,ю2) - и ('1,юк ,»2)
следовательно, критерием остановки итерационного процесса (10) может являться точное число итераций т, которое может быть определено из уравнения
,ю2)и0(ю1Г >ю2) -и ('1,ю^ ,ю2)
<£,
где е - заданная точность.
Определив От (¿1,юк ,“2) с точностью е, вычислим значения
От (¿1, Ук, у/), к, / = 0, N -1, используя кубатурные формулы вычисления обратного преобразования Фурье:
1 N-1 N-1
От (¿1, х1, х2) = — £ £ От (¿1,Ю1?1 ,Ю^2) Ц е^^+^йю^ .
Л =0 /2 =0
Л ./2
Выбрав подходящую сетку узлов по переменным х1 и х2 и используя стандартные методы теории приближения: интерполяционные полиномы,
сплайны и т.д., - восстанавливаем функцию Грина От(¿1,х1,х2), после чего, располагая, таким образом, всеми необходимыми значениями, можно легко восстановить функцию и(¿1, х1, х2) по формулам (4).
Замечание 1. Если коэффициенты аи Ьк в уравнении (1) не зависят от времени, то приведенный выше алгоритм позволяет восстановить функцию и(¿, х1, х2) при ^ е [¿0, Т] по значениям функций и(¿0, х1, х2) и и(¿1, х1, х2), причем Т > ¿1. Аналогично, если от времени не зависят коэффициенты а^1, Ьк и функция /(¿, х) уравнения (3), то приведенный выше алгоритм позволяет восстановить функцию и(¿, х1, х2) при I е [¿0, Т] по значениям функций и (¿1, х1, х2) и /(х).
Это следует из того, что, согласно уравнению (6),
где п - целое положительное число.
Это позволяет получить значения функции и(^,“1,02) на сетке
Замечание 2. Выше отмечалось, что восстановление функции G(t, x) является некорректной задачей. Ее регуляризация проводится за счет того, что, во-первых, прямое и обратное преобразования Фурье вычисляются в конечных пределах; во-вторых, за счет того, что для решения систем уравнений (9) используется итерационный процесс, обладающий свойствами стабилизации и фильтрации.
Рассмотрим алгоритм восстановления функции G (t, x) в динамической системе, описываемой уравнением (3) при нулевых начальных и граничных условиях.
В случае, если в уравнении (3) функция f (t, x) зависит от времени, то для восстановления функции G (t, x ) требуется информация о функциях u (t, x) и f (t, x) на сетке, составленной из узлов по временным и пространственным переменным. Для простоты обозначений рассмотрим прямоугольную сетку (tv,yki,...,ykn }, v = 0,M , kt = 0,N, i = 1,n, tv = vt, где т - шаг по временной переменной; yk. = —A + hki, h = 2A/N - шаг по каждой из пространственных переменных.
Представим уравнение (6) в виде
U (nti ,Ю1 ,Ю2) = (G(ti ,o>i ,ro2))nUo(œ>i ,о>2),
(^1,01,02), а затем восстановить функцию и(¿, х1, х2) на каждом временном слое ^ , v = 1, т , т = Т +1.
¿1
'П
K (t — t)G (t — t, x -£) f (T,£)d t, (12)
Применим к уравнению (12) преобразование Фурье по пространственным переменным. В результате получим
и (',ю) =—I" К ІЇ-т)О(1 -т,ю)^(т,ю)^ т. (13)
(2п)п/2 і
Введем для краткости записи обозначение
К(' - т)О(' -т, ю) = О1 (' - т, ю) и применим к уравнению (13) преобразование Фурье по переменной '. В результате получим
и(Х,ю) = О^Х, ю)Е(Х, ю). (14)
Введем сетки узлов по переменным Х и ю :
^ = {^1,,.ДМ}, юг- = {ю1 }, і = 1,п . (15)
Воспользовавшись кубатурными формулами вычисления преобразований Фурье, находим значения функций и (Х,ю) и Е (Х,ю) на сетке узлов (15).
Приравнивая левые и правые части уравнения (14) на сетке (15), получаем систему из МЫп линейных алгебраических уравнений, каждое из которых решается в отдельности. Зафиксировав некоторый узел (Хк, ю^,..., юПп), рассмотрим уравнение
и(Хк, ю1,..., ю1пп) = О1(Хк, ю^1, юПп)Е(Хк, ю^,..., юПп).
Значение ^1(Хк, ю^, .,юПп) определим итерационным методом, аналогичным процессу (10):
ОГ+1 (Хк,ю^1, юПп) = аГОГ (Хк,ю^1,...,юПп) + (1 -ат)х
ОГ (Хк X1 ,.,юПп)-Ук, А>іп {От (Хк ,ю1 ,.,юПп) Е (Хк ,ю1 ,.,юПп) -
-и(Хк,ю}, ...,ю1пп)}
(16)
где 0 <а*<ат <а*<1, а параметры у к,/ ,/ подбираются таким образом, чтобы выполнялись условия
1 <Ук,..., іпЕ (Хк X1 >.,юПп) <4.
Сходимость итерационного метода (16) также следует из теоремы Об-ломской и исследуется аналогично сходимости итерационного процесса (10).
Обозначим через О*(кк, 011,..., ) предельное значение итераций (16)
или, что больше соответствует вычислительной практике, значение, на котором останавливается данный итерационный процесс. Вычислив значения
функции О*(кк,01, ..., 0)п) на сетке узлов (15), по кубатурным формулам вы-
числения обратного преобразования Фурье во временной области находим функцию G\(t,ю) и тем самым функцию G(¿,ю). Окончательно, применяя кубатурные формулы вычисления обратного преобразования Фурье в пространственной области, находим функцию G(t, х).
Задача 2. Остановимся на еще одном аспекте применения предложенного выше алгоритма. Пусть динамическая система описывается задачей Коши (1), (2) с неизвестными коэффициентами ak j (t), k, j = 1, n и bk (t),
k = 1, n . Требуется, располагая начальными значениями Hq (x) и решением задачи Коши u (t), t > 0, найти приближенные значения функций ak j (t) и
bk (t) при 0 < t < T , T = const.
Представим искомое решение задачи Коши (1), (2) в виде обратного преобразования от функции U (^ю), ю = (ю>1,..., юп):
u(t, x) =--1—~ ITe~i(х’ю)и(t,o>)dю, (17)
(2n)n/2 E
En
n
где (х,ю) = Z Xi Юi . i=1
Подставив (17) в уравнение (1) и предполагая, что можно провести все необходимые действия под знаком интеграла, приходим к обыкновенному дифференциальному уравнению
dt
при начальном условии
dU (t, ю)
_ Z ak,j(t)юkюj +1Z bk(t)юk+c(t)
k,j=1 k =1
(18)
U(0, ю) = —ife-i(x,ю)Uo(x)dx, (19)
(2n)n/2 E
En
здесь ю = (ю^, ..., юп) является параметром. Задача Коши (18), (19) имеет решение
U (^ю) = U (0,ю) exp
Z ak,j(t)(%(aj+i Z bk (t),%+c(t)
k ,j =1 k=1
(20)
Введем сетку узлов ¿к = N, к = 0,п . Из уравнения (20) найдем приближенные функции ак у ((), к, у = 1, п, Ък ((), к = 1, п и с ((), аппроксимирующие функции ак,у ((), Ък (() и с(() соответственно. Здесь подразумевается, что ак у (() - кусочно-постоянная функция, равная константе а\ у на сегменте [¿о, ^1 ] и равная константе ау,у на интервале (¿у_1, ], V = 2, N . Анало-
гичным образом определяются функции Ък ((), к = 1, п и с ((). Таким образом, при каждом фиксированном значении (к , к = 1, N из уравнения (20)
2
нужно найти (п + п + 1) неизвестных чисел.
Положив в уравнении (20) ( = (к и ю = (ю1,юп) = (0,.., 0), получаем
и (¿к, 0) = и(0,0)ес((к). (21)
Тогда уравнение (20) можно записать в виде
U (t, ю) = U * (0, ю) exp
£ ak, j(t)юкю] + і £ bk(t)юк
k,j=l к=1
(22)
где U (0,ю) = U(0,0)exp{c)}.
Для нахождения n2 значений ak j, к, j = 1, n при каждом фиксированном значении v , v = 1, N, поступим следующим образом. В области й = [-A, A]n выберем n2 узлов Мк = (<4-,<n), к = 1,n2 . Представим значения U(fy,<) и U * (0,юк) в виде
U(ti,<) = exp{ln U(ti,<)} = exp{ln | U(ti,<) | +i arg U(ti,<)};
U * (0,тк) = exp{ln U * (0,тк)} = exp{ln | U * (0,тк) | +i arg U * (0,тк)},
I = 0, N, к = 1, п2 .
Замечание 3. При этом предполагается, что функции и((/ ,юк) и и * (0,юк) не обращаются в ноль на рассматриваемой сетке узлов. Это ограничение является несущественным, т.к. функции и ((I ,юк) и и * (0,юк) являются аналитическими по переменной ю вне координатных плоскостей и, следовательно, может обращаться в ноль вне этих плоскостей только на счетном множестве точек, не включающем в себя предельную точку.
Из системы уравнений
ln | U(ti,Mk) | - ln | U* (0,Mk) |= £ af> j (ti К тк, к = 1, n2
i,j=1
определяем значения a^j (ti), i, j = 1, n .
Для определения функций Ьк (t), к = 1, n будем учитывать, что разность аргументов
n
arg U(ti ,юк) - arg U* (0,юк) = £ Ск (t)<% .
к=1
Нетрудно видеть, что функции Ck (t) являются приближениями к функциям bk (t), k = 1, n .
Замечание 4. В случае, если функция (argU(t/ ,ff>k) - argU* (0,rok)) известна только на сетке rok , k = 1, N, то аналогичным образом определяются приближенные значения b/ (tk), l = 0, n .
Пример. В качестве примера, иллюстрирующего изложенный выше метод восстановления функции Грина, рассмотрим процесс, описываемый уравнением
г A+_d2. ^
dt + дх2
u(t, х) = f (t, х), хє [-0.5,0.5], t є [0,1] (23)
с нулевыми начальными и граничными условиями.
Будем считать известными значения функций u (t, x) и f (t, x). Методами линейной алгебры данный процесс был исследован в работе [3]. Там же было показано, что при
( П
f (x, t) = n(cos nt - 4n sin nt) cos I 2nx - —
точным решением уравнения (23) при нулевых начальных и граничных значениях является функция
u(x, t) = sin nt cos 12nx - П
Полагая u(t, x) = 0 при x í [-0,5; 0,5], уравнение (23) можно рассматривать при -^ < x < ^ , t > 0 .
Поскольку рассматривается дифференциальный оператор
(_д+_д^ ^ dt + дх2
постоянными коэффициентами, то решение уравнения (23) можно представить в виде
и(I,х) = — Г йт Г — т,х-£)/(т,£)й£, (24)
2п
2п
где
G (t, х) = <
— expj-~_ і = G(t, х), при t > 0,
0, при I < 0.
Для определения функции О х) выберем сетки узлов
х^ = —А + 2АкN, к = 0,N и = ут, V = 0,1,..., где А - достаточно большое положительное число, а т - достаточно малое положительное число, являющееся временным шагом. Будем считать, что на данной сетке известны значения и х) и / х).
Применяя к уравнению (24) преобразования Фурье по переменным ^ и х, приходим к уравнению
где X - переменная, соответствующая времени ^; ю - переменная, соответствующая пространственной переменной х.
Вычислив значения функций и (X, ю) и Е (X, ю) на указанной сетке узлов, решаем итерационным методом (15) систему уравнений
после чего, используя кубатурные формулы вычисления обратного преобразования Фурье, находим функцию О(ґ, х), а затем, используя полученные данные, находим функцию и (ґ, х).
Погрешности восстановления функций О(ґ, х) и и (ґ, х) представлены на рис. 1, 2.
2 Приближенный метод восстановления начальных значений
Во многих областях физики и техники возникает задача определения начальных значений в динамических процессах, описываемых дифференциальными уравнениями в частных производных. Как известно, эта задача некорректна и для ее решения для каждого класса уравнений предлагались специфические методы регуляризации [6].
Приведем один из возможных методов определения начальных значений на примере уравнения теплопроводности.
Рассмотрим уравнение
и (X, ю) = О (X, ю) ^ (X, ю),
(25)
и (Х к ,юу ) 0 (Х к ,юу ) ^ (Х к ,юу ),
х
Рис. 1
Рис. 2
ди (ґ’х) -Ди (ґ,х) = 0, х єО, ґ > 0,
дґ
(26)
где
с начальным условием
Известия высших учебных заведений. Поволжский регион
и(0, х) = Мд( х) (27)
и граничным условием
и($, х) = 0, при хе Г = дй, I > 0, (28)
где щ (х) - функция, определенная в области й; Г = дй - граница области й .
Требуется, располагая значениями функции и(х, ^) при ^ = Т , восстановить функцию щ($). Предположим вначале, что й = Яп . Хорошо известно [4], что если ) - непрерывная функция, то решение уравнения (26) при ^ > 0 определяется интегралом Пуассона:
u(t,x) =------— Г ... Г exp
{2у[Лл)п J J
(Х1 — ^1)2 + ■■■ + (xn — ^и )2
4t
X
xu0(^i, ...,£и)dh-dIn. (29)
Таким образом, при Q, = Rn определение функции uo(t) сводится к решению интегрального уравнения в свертках (29). Так как размерность данного уравнения не имеет принципиального значения при построении вычислительных схем, подробно приведем случай, когда n = 1, т.е. ограничимся рассмотрением уравнения
^ [ 2 "I
u(t’Х) = 27ПТ ^ ехр | (Х 4t) I Uo(^)d^ . (30)
Учитывая, что функция u (t, x) известна при t = T, представим это уравнение в виде
u(T’x)^V2ni )1техР{(ХT ju0())d) . (31)
Применим к данному уравнению преобразование Фурье по пространственной переменной. В результате получаем уравнение
U (T ,ю) = H (ю)и0(ю), (32)
где U(T,ю), H(ю), Uо(ю) - преобразование Фурье функций u(x,T),
f(x) = ехр {и uo(t) соответственно.
4Т
Уравнение (32) будем решать итерационным методом, описанным в предыдущем разделе. В результате получим значения функции и0(ю^) на
сетке Ю£ , к = 1, п . Используя эти значения и квадратурные формулы вычисления обратного преобразования Фурье, найдем значения Нд(хі), і = 1, п на заранее заданной сетке {хі}, і = 1, п .
—оо
Приведем конкретный пример.
Пример. Рассмотрим уравнение (30) с начальной функцией
Гяп(2ях), х є [-1,1]
и (ґ, х) = <
[ 0, х £[-1,1].
Полагаем, что функция и (ґ, х) известна на равномерной сетке узлов
2к -------------
Хк =-1 + "^ ’ к = 0, N, N = 103, при ґ = 0,2. Требуется, располагая значениями функции и (ґ, х), на данной сетке восстановить функцию и0(х). Результат восстановления указанной функции и0(х) по предложенному в данном разделе алгоритму приведен на рис. 3.
/ 0.004- •. 8
-о!м ^<¿72 ^■е-.Зз -о!з4 -0.004- Г-—-__0л Л- 0.43 0*72 \ 0І96
Рис. 3
На рис. 3 е - абсолютная погрешность восстановления функции мд(х), а 5 - абсолютная погрешность восстановления функции и (х) по значениям новой функции {¡о (х), полученной в ходе реализации алгоритма.
Выводы. В работе предложен достаточно общий алгоритм идентификации динамических систем и восстановления входных сигналов динамических систем, обладающий следующими достоинствами:
1) высокая точность;
2) возможность распараллеливания;
3) быстродействие;
4) небольшой объем памяти.
Список литературы
1. Лаврентьев, М. М. Некорректные задачи математической физики / М. М. Лаврентьев, В. Г. Романов, С. П. Шишатских. - М. : Наука : 1980. - 288 с.
2. Романов, В. Г. Обратные задачи математической физики / В. Г. Романов. - М. : Наука, 1984. - 263 с.
3. Скопецький, В. В. Про задачу ідентифікації та керування для дискретно спостережуваного просторово-часового процессу / В. В. Скопецький, В. Н. Стоян // Reports of the Akedemy of Sciences of Ukraine. - 2006. - № 8. - Р. 102-108.
4. Эйдельман, С. Д. Параболические системы / С. Д. Эйдельман. - М. : Наука, 1964. - 444 с.
5. Обломская, Л. Я. О методах последовательных приближений для линейных уравнений в банаховых пространствах / Л. Я. Обломская // Журнал вычислительной математики и математической физики. - 1968. - Т. 8. - № 2. - С. 417-426.
6. Латтес, Р. Метод квазиобращения и его приложения / Р. Латтес, Ж.-Л. Лионс. -М. : Мир, 1970. - 336 с.