САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
ФАКУЛЬТЕТ ПРИКЛАДНОЙ МАТЕМАТИКИ — ПРОЦЕССОВ УПРАВЛЕНИЯ
КАФЕДРА ТЕОРИИ УПРАВЛЕНИЯ
Чернышева Любовь Андреевна
Выпускная квалификационная работа бакалавра
Матрица Ляпунова как решение уравнения
Фредгольма
Направление 010400
Прикладная математика, фундаментальная информатика
и основы программирования
z
Научный руководитель,
кандидат физ.-мат. наук,
доцент
Егоров А. В.
Санкт-Петербург
2016
Содержание
1. Введение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
2. Основные используемые понятия и теоремы . . . . . . . . . . . . . . .
6
3. Матрица Ляпунова как решение интегрального уравнения Фредгольма 2-го рода . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
4. Методы решения интегрального уравнения Фредгольма 2-го рода . . 16
4.1. Метод последовательных приближений . . . . . . . . . . . . . . . 16
4.2. Квадратурный метод . . . . . . . . . . . . . . . . . . . . . . . . . 18
5. Практическая реализация . . . . . . . . . . . . . . . . . . . . . . . . . 20
6. Заключение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
7. Список литературы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2
Введение
Предметом данной работы является построение матрицы Ляпунова для
систем дифференциально-разностных уравнений с запаздывающим аргументом. Подобные системы используются для описания многочисленных процессов, использующих передачу массы, энергии, информации и т. п. Существует множество факторов, таких как ограниченность скорости распространения взаимодействия, наличие инерционности отдельных элементов, ограниченность скорости протекания физических процессов, в связи с которыми появляется запаздывание в процессе. Иногда этим запаздыванием можно пренебречь и описать систему при помощи системы обыкновенных дифференциальных уравнений. Однако неучитывание запаздывания зачастую может
привести к результатам, несоответствующим реальности. Поэтому целесообразно использовать дифференциальные уравнения, в которых функция и ее
производные входят при разных значениях аргумента.
Известно, что для дифференциальных уравнений с запаздывающим аргументом при выполнении условия Ляпунова существует единственная матрица
Ляпунова [1]. В последние 10 лет проводились исследования матрицы Ляпунова, было дано новое определение [1], [3], а также найден ряд ее применений,
таких как исследование устойчивости [2], построение экспоненциальных оценок решений [1], исследование робастной устойчивости [3], вычисление значения функционала качества [1], вычисление H2 [6] нормы и другие.
Однако до сих пор не был найден универсальный метод построения этой
матрицы. Нахождение точного значения матрицы Ляпунова является непростой задачей, поэтому существуют различные методы вычисления матрицы,
например, полуаналитический [1], полиномиальная аппроксимация [5], [6],
3
кусочно-линейная аппроксимация [1].
Полуаналитический метод применим только для систем с кратными запаздываниями вида
ẋ(t) =
m
X
Aj x(t − kh),
k ∈ N.
j=0
Было показано, что вычисление матрицы Ляпунова может быть сведено к построению решения специальной краевой задачи для линейной системы обыкновенных дифференциальных уравнений [1]. При выполнении условия Ляпунова такая задача имеет единственное решение [1]. Недостатком этого метода является отсутствие оценки погрешности вычислений, а также тот факт,
что этот метод применим только для систем с кратными запздываниями. Как
правило, необходимо вводить дополнительные запаздывания с нулевыми матрицами, чтобы запаздывания в системе стали кратными, что в свою очередь
влечет значитльное увеличение времени работы метода, а также погрешности
вычислений.
Полиномиальная аппроксимация матрицы Ляпунова осуществима также
без оценки погрешности вычислений.
Численный метод построения кусочно-линейной аппроксимации матрицы
Ляпунова состоит из двух этапов. На первом этапе вычисляется кусочнолинейная аппроксимация начальных условий для матрицы Ляпунова. Затем
эти начальные условия используются при вычислении приближенного значения матрицы Ляпунова как решения уравнения
m
dU (τ ) X
=
U (τ − hj )Aj , τ > 0.
dτ
j=0
4
Целью данной работы является построение матрицы Ляпунова для систем
с некратными запаздываниями, а также получение оценки погрешности, применимой на практике.
Один из методов решения поставленной задачи это сведение задачи к системе, решение которой будет удовлетворять динамическому свойству, свойству симметрии и алгебраическому свойству из определения матрицы Ляпунова. Полученное интегральное уравнение будем решать численно при помощи метода последовательных приближений и квадратурного метода, используя квадратурную формулу трапеций. Для метода последовательных приближений будет найдена оценка погрешности.
5
Основные используемые понятия и теоремы
Определение 1 [11]. Дифференциальными уравнениями с отклоняющимся аргументом называются дифференциальные уравнения, в которых неизвестная функция и ее производные входят при различных значениях аргумента.
Рассмотрим систему дифференциально-разностных уравнений вида:
m
X
ẋ(t) =
Aj x(t − hj ), t > 0.
(1)
j=0
Здесь Aj , j = 0, . . . , m, — матрицы размерности n × n; hj , j = 0, . . . , m, —
запаздывания, упорядоченные по возрастанию 0 = h0 < h1 < h2 < . . . < hm =
H. Запаздывания в поставленной задаче являются некратными, не зависят от
времени.
Пусть ϕ : [−H, 0] → Rn — начальная функция, принадлежащая пространству P C([−H, 0], Rn ).
Определение 2 [1]. Пусть x(t, ϕ) — решение задачи Коши для системы
(5) с начальными условиями:
x(θ, ϕ) = ϕ(θ),
θ ∈ [−H, 0].
Требуется построить матрицу Ляпунова для рассматриваемой системы
дифференциально-разностных уравнений.
Определение 3 [1]. Будем говорить, что U (τ ) — матрица Ляпунова для
системы (5), если она удовлетворяет следующим свойствам:
1. Динамическое свойство:
m
dU (τ ) X
U (τ − hj )Aj , τ > 0;
=
dτ
j=0
6
(2)
2. Свойство симметрии:
U (−τ ) = U T (τ );
(3)
3. Алгебраическое свойство:
m
X
[U (−hj )Aj + ATj U (hj )] = −W.
(4)
j=0
Для построения матрицы Ляпунова необходимо убедиться в ее существовании и единственности. Для этого воспользуемся условием Ляпунова.
Определение 4
[1]. Будем говорить, что система (5) удовлетворяет
условию Ляпунова, если спектр системы,
!
)
(
m
X
Λ = s | det sI −
e−shj Aj = 0 ,
j=0
не содержит точки s0 , такой что −s0 также принадлежит спектру.
Теорема 1 (О существовании и единственности матрицы Ляпунова) [1]. Матрица Ляпунова для системы (5), ассоциированная с данной
симметрической матрицей W , существует и единственна тогда и только
тогда, когда система удовлетворяет условию Ляпунова.
Определение 5 [1]. Будем говорить, что матрица K(t) размерности n ×
n — фундаментальная матрица системы (5), если выполнено:
m
X
d
K(t) =
K(t − hj )Aj ,
dt
j=0
t > 0,
K(t) = 0n×n для t < 0, K(0) = I.
Определение 6 [1]. Фундаментальная матрица системы (5) также удовлетворяет равенству:
m
X
d
K(t) =
Aj K(t − hj ),
dt
j=0
7
t > 0.
Теорема 2 [1]. Дана начальная функция ϕ ∈ P C([−H, 0], Rn ). Выполнено
следующее равенство, называемое формулой Коши:
0
x(t, ϕ) = K(t)ϕ(0) +
m Z
X
K(t − θ − hj )Aj ϕ(θ)dθ,
t > 0.
j=1 −h
j
Здесь x(t, ϕ) — решение системы (5).
Поиск решения системы (5) осуществляется при помощи метода последовательного интегрирования [8].
Определение 7 [7]. Системой интегральных уравнений Фредгольма 2-го
рода называется равенство:
ZH
Y (s) = G(s) +
F (s, x)Y (x)dx.
0
Здесь Y (s) — искомая функция, а F (s, x) и G(s) — известные функции,
заданные на отрезке [0, H]. Функцию F (s, x) будем называть ядром интегрального уравнения.
Определение 8 [1]. Функцией Хевисайда будем называть функцию
1, t > 0;
χ(t) =
0, t < 0.
Будем использовать равномерную норму:
kϕk∞ = sup kϕ(t)k.
t∈[0,H]
8
(5)
Определение 9
[5]. Кронекеровым произведением матриц будем назы-
вать операцию:
b A
11
b12 A
A⊗B =
..
.
b1n
T
c11 d1
.
..
cq1 dT
1
.
T
.
C⊗D=
.
c d T
11 m
..
.
cq1 dTm
b21 A · · · bn1 A
b22 A · · ·
..
...
.
b2n A · · ·
c12 dT1 · · ·
..
...
.
cq2 dT1 · · ·
..
...
.
c12 dTm · · ·
..
...
.
cq2 dTm · · ·
bn2 A
;
..
.
bnn A
T
c1p d1
..
.
T
cqp d1
..
.
.
T
c1p dm
..
.
T
cqp dm
(6)
(7)
Определение 10 [1]. Будем говорить, что ũ = vec(U ) операция векторизации матрицы U, если ũ получается путем размещения последовательно
столбцов матрицы один под другим. Эта операция удовлетворяет равенству:
vec(AU B) = (A ⊗ B)ũ
и равенству
T
vec(CU T D) = (C ⊗ D)ũ.
9
Матрица Ляпунова как решение интегрального
уравнения Фредгольма 2-го рода
Рассмотрим дифференциально-разностное уравнение вида
ẋ(t) = A0 x(t) +
m
X
Aj x(t − hj ),
(8)
j=1
где Aj , j = 0, . . . , m, — заданные матрицы размерности n × n; hj , j =
0, . . . , m, — запаздывания, упорядоченные по возрастанию 0 = h0 < h1 <
h2 < . . . < hm = H.
Вычислим для ξ ∈ (0, t) частную производную
∂[U (ξ)K(t − ξ)]
.
∂ξ
Будем использовать динамическое свойство
m
dU (τ ) X
=
U (τ − hj )Aj , τ > 0.
dτ
j=0
Получим
∂[U (ξ)K(t − ξ)]
= U (ξ)A0 K(t − ξ)+
∂ξ
m
X
+
U (ξ − hj )Aj K(t − ξ) − U (ξ)A0 K(t − ξ)−
j=1
−
m
X
U (ξ)Aj K(t − ξ − hj ) =
j=1
=
m
X
[U (ξ − hj Aj K(t − ξ) − U (ξ)Aj K(t − ξ − hj )] .
j=1
Проинтегрируем это равенство по ξ на отрезке [0, t].
10
Получим, что интеграл левой части равен
Zt
∂[U (ξ)K(t − ξ)]
dξ = U (t)K(t − t) − U (0)K(t) =
∂ξ
0
= U (t) − U (0)K(t).
Интеграл правой части записывается в виде
J(t) =
Zt X
m
[K(t − ξ)U (t − hj )Aj − K(t − ξ − hj )Aj U (ξ)] dξ.
j=1
0
Или
J=
m Z
X
t
[U (ξ − hj )Aj K(t − ξ) − U (ξ)Aj K(t − ξ − hj )] dξ =
j=1 0
= hθ = ξ − hi =
t−h
m Z j
X
U (θ)Aj K(t − θ − hj )dθ −
j=1 −h
0
+
t
U (ξ)Aj K(t − ξ − hj )dξ+
j=1 0
j
m Z
X
m Z
X
U (θ)Aj K(t − θ − hj )dθ −
j=1 −h
m Zt
X
U (ξ)Aj K(t − ξ − hj )dξ.
j=1 t−h
j
j
Так как t − ξ − hj < 0 для ξ ∈ (t − hj , t], K(t − ξ − hj ) = 0n×n . Получаем
равенство
0
J=
m Z
X
U (θ)Aj K(t − θ − hj )dθ.
j=1 −h
j
В результате заключаем, что
0
U (t) − U (0)K(t) =
m Z
X
U (θ)Aj K(t − θ − hj )dθ.
j=1 −h
j
11
Отсюда
0
U (t) = U (0)K(t) +
m Z
X
U (θ)Aj K(t − θ − hj )dθ,
t > 0.
(9)
j=1 −h
j
Далее применим к полученному выражению свойство симметрии
U (−τ ) = U T (τ ).
Чтобы матрица Ляпунова вычислялась на положительной оси, необходимо
сделать замену
Z0
U T (−θ)K(t − θ − hj )dθ = hη = −θi =
−hj
Zhj
U T (η)Aj K(t + η − hj )dη.
0
В итоге получим
hj
U (t) = U (0)K(t) +
m Z
X
U T (η)Aj K(t + η − hj )dη,
t > 0.
(10)
j=1 0
Осталось применить третье свойство — алгебраическое:
m
X
[U (−hj )Aj + ATj U (hj )] = −W,
(11)
j=0
где W — произвольная положительно определенная матрица соответствующей размерности.
Чтобы учесть алгебраическое свойство необходимо сначала векторизовать
уравнения (10) и (11).
12
Равенство (10) примет вид:
hj
ũ(t) = (I ⊗ K(t))ũ(0) +
m Z
X
T
(I ⊗ Aj K(t + η − hj ))ũ(η)dη,
t > 0.
(12)
j=1 0
Алгебраическое свойство (11) запишется в виде:
m
X
T
[(I ⊗ Aj )ũ(hj ) + (ATj ⊗ I)ũ(hj )] = −w̃.
(13)
j=1
В уравнении (13) неизвестными являются ũ(hj ) и ũ(0), где ũ(hj ) можно
получить по формуле (12), подставив hj вместо t. Затем необходимо из алгебраического свойства выразить ũ(0) и подставить в исходную формулу Коши
(12).
В итоге получим систему интегральных уравнений Фредгольма 2-го рода:
ZH
ũ(τ ) = G(τ ) + F (τ, η)ũ(η)dη.
(14)
0
Опишем более подробно процесс получения ядра F (τ, η) и функции G(τ )
на примере скалярного уравнения.
Формула Коши в скалярном случае c учетом динамического свойства и
свойства симметрии будет выглядеть следующим образом:
u(t) = k(t)u(0) +
m
X
Zhj
k(t + η − hj )u(η)dη.
aj
j=1
0
Алгебраическое свойство в скалярном случае примет вид:
m
X
aj u(hi ) = −w/2.
j=0
Здесь w — любое положительное число.
Найдем u(hi ) по уже полученной формуле (15)
u(hi ) = k(hi )u(0) +
m
X
(15)
Zhj
k(hi + η − hj )u(η)dη.
aj
j=1
0
13
(16)
Подставим получившееся выражение в алгебраическое свойство (16)
h
j
Z
m
m
X
X
aj k(hi + η − hj )u(η)dη = −w/2.
k(hi )u(0) +
i=0
j=1
0
Или
a0 u(0) +
m
X
ai k(hi )u(0) +
i=1
Zhj
m X
m
X
k(hi + η − hj )u(η)dη = −w/2.
ai aj
i=1 j=1
0
Далее выразим u(0):
u(0)(a0 +
m
X
ai k(hi )) +
m
m X
X
Zhj
i=1 j=1
i=1
w/2 +
0
m P
m
P
ai aj
i=1 j=1
u(0) = −
k(hi + η − hj )u(η)dη = −w/2,
ai aj
a0 +
Rhj
k(hi + η − hj )u(η)dη
0
m
P
.
ai k(hi )
i=1
Полученное выражение для u(0) подставим в исходную формулу Коши
(15):
w/2 +
m
m P
P
ai aj
i=1 j=1
u(τ ) = −k(τ )
a0 +
Rhj
k(hi + η − hj )u(η)dη
0
m
P
+
ai k(hi )
i=1
+
m
X
j=1
Zhj
k(t + η − hj )u(η)dη.
aj
0
Или
u(τ ) = −
k(τ )w/2
+
m
P
a0 +
ai k(hi )
i=1
14
+
m
X
Zhj
aj
j=1
m
P
k(τ + η − hj ) −
i=1
ai k(τ )k(hi + η − hj )
a0 +
0
m
P
u(η)dη, τ > 0.
ai k(hi )
i=1
Воспользуемся функцией Хевисайда χ(t) так, чтобы получить один интеграл с пределами от 0 до наибольшего из запаздываний H.
u(τ ) = −
k(τ )w/2
+
m
P
a0 +
ai k(hi )
i=1
+
Zhj X
m
0
j=1
m
P
aj
k(τ + η − hj ) −
ai k(τ )k(hi + η − hj )
i=1
a0 +
m
P
ai k(hi )
χ(hj −η)u(η)dη, τ > 0.
i=1
Пусть
k(τ )w/2
g(τ ) = − P
;
m
ai k(hi )
m
X
aj k(τ + η − hj ) −
f (τ, η) =
i=0
m
P
ai k(τ )k(hi + η − hj )
i=1
m
P
j=1
ai k(hi )
χ(hj − η) .
i=0
Тогда
ZH
u(τ ) = g(τ ) +
f (τ, η)u(η)dη.
(17)
0
Полученное выражение является интегральным уравнением Фредгольма
2-го рода.
15
Методы решения интегрального уравнения Фредгольма
2-го рода
Метод последовательных приближений
Метод последовательных приближений применяется для решения интегральных уравнений Фредгольма 2-го рода. Каждое последующее приближение
находится по формуле:
ZH
f (τ, η)un (η)dη.
un+1 (τ ) = g(τ ) +
0
Первое приближение u0 выбирается произвольным образом.
Оценку погрешности можно получить при помощи принципа сжимающих
отображений. Введем оператор:
ZH
f (τ, η)u(η)dη.
A(u)(τ ) = g(τ ) +
0
Известно [10], что
kA(u) − A(v)k 6 qku − vk,
где 0 6 q < 1.
Необходимо оценить q:
H
Z
6
kA(u)(τ ) − A(v)(τ )k =
f
(τ,
η)(u(η)
−
v(η))dη
0
H
Z
max ku(η) − v(η)k.
6
f
(τ,
η)dη
η∈[0,H]
0
16
Отсюда
H
Z
.
f
(τ,
η)dη
q=
0
Если 0 6 q < 1, то A — сжимающий оператор. В таком случае можно
говорить о сходимости метода.
Пусть u — точное решение, un — приближенное решение, полученное на
n-ом шаге.
Тогда
kun − uk 6 kun+1 − un k + kun+1 − uk 6
6 kun+1 − un k + kA(un ) − A(u)k 6 kun+1 − un k + qkun − uk.
Отсюда
kun − uk 6
1
kun+1 − un k.
1−q
Также можно можно получить оценку погрешности для n+1 приближения:
kun+1 − uk 6 qkun − uk 6
q
kun=1 − un k.
1−q
Таким образом, методом последовательных приближений можно построить приближенное значение матрицы Ляпунова с любой наперед заданной
точностью. Существенным недостатком этого метода является жесткое условие сходимости, так как оно выполняется только при малых значения коэффициентов. Метод достаточно прост в практической реализации, но из-за
того, что на каждой итерации, а так же при вычислении оценки погрешности необходимо вычислять значения нескольких интегралов, метод сходится
довольно медленно.
17
Квадратурный метод
Численное решение интегрального уравнения Фредгольма 2-го рода можно
также построить методом квадратур. Для этого необходимо взять на отрезке
[0, H] сетку с узлами τ1 , τ2 , . . . , τn . Сетка равномерная с шагом h. Уравнение
(17) в узлах сетки запишется следующим образом:
ZH
u(τi ) = g(τi ) +
f (τi , η)u(η)dη,
i = 1, 2, . . . , n.
(18)
0
Аппроксимируя интегралы в равенствах (18) конечными суммами при помощи квадратурной формулы, получим:
ui = gi +
n
X
Bk fik uk ,
i = 1, 2, . . . , n.
(19)
k=1
Здесь ui = ũ(τi ), ũ — приближение к искомой функции u; Bk — веса квадратурной формулы; fik = f (τi , τk ).
При аппроксимации интегралов конечными суммами будет использоваться
квадратурная формула трапеций:
ui = gi + h
n
X
wk fik uk ,
i = 1, 2, . . . , n.
(20)
k=1
w1 = wn = 1/2, wk = 1, k = 2, 3, . . . , n − 1.
Пусть ∆ — определитель системы (19), а ∆ki — алгебраическое дополнение
элемента определителя ∆ с индексами k, i.
Обозначим
ZH
f (τ, η)g(η)dη −
pf (τ ) =
n
X
k=1
0
18
Bk g(τk ),
ZH
f (τ, η)f (η, v)dη −
pk (τ, v) =
n
X
Bk K(τk , v).
k=1
0
pf = max |pf (τ )|,
τ
ZH
|pk (τ, v)|.
pk = max
τ
0
Погрешность ε вычислений квадратурного метода находится по формуле:
ε 6 M (pf + pk S).
Здесь
(21)
n
1 X
|∆ki |,
B>
∆
k=1
S — максимум точного решения.
Верхнюю оценку решения можно получить, например, из другого грубого
метода.
Оценку квадратурного метода можно получить и без использования максимума решения [7]. Получение оценки погрешности без использования максимума точного решения является основным направлением дальнейших исследований.
19
Практическая реализация
Пример 1
Рассмотрим уравнение:
ẋ(t) = −12.9x(t) − 12.45x(t − 0.04) + 10x(t − 0.08),
w = 1.
Для этого примера условие сходимости метода последовательных приближений выполняется. Зададим точность ε = 0.001.
20
Методом последовательных приближений было сделано 2 итерации, значение q для данной системы равно 0.2969. Время работы программы 757 секунд
или около 13 минут.
Пример 2
Рассмотрим теперь дифференциально-разностное уравнение с двумя некратными запаздываниями:
ẋ(t) = −x(t) − 2x(t − 1) − 3x(t −
√
2),
w = 1.
Метод последовательных приближений на этом примере не сходится. Точность вычислений ε = 0.005.
21
Стоит отметить, что, так как полуаналитический метод не работает на
√
некратных запаздываниях, то в данном примере вместо 2 было взято число 1.4, а запаздывания 0.1, 0.2,. . . , 1.4. Таким образом, уравнение для полуаналитического метода имело 14 запаздываний. При увеличении точности
вычисления операции извлечения корня на каждый знак после запятой количество запаздываний в полуаналитическом методе будет увеличиваться в 10
√
раз. Так, например, взяв 2 = 1.4142, количество запаздываний будет 14000,
что существенно замедляет работу программы.
Время работы квадратурного метода составляет 47 секунд.
22
Заключение
В данной работе для систем дифференциально-разностных уравнений с
запаздывающим аргументом была построена матрица Ляпунова методом последовательных приближений и методом квадратур. Также в работе были
приведены оценки погрешностей численных методов и рассмотрен численный
пример.
В качестве направлений дальнейших исследований отметим нахождение
более точной оценки погрешности для квадратурного метода, а также реализацию решения интегрального уравнения Фредгольма 2-го рода для систем
дифференциально-разностных уравнений.
23
Список литературы
1. Kharitonov V. L. Time-delay systems: Lyapunov functionals and matrices //
Birkhauser, Basel. 2013.
2. Egorov A.V., Mondie S. Necessary stability conditions for linear delay systems
// Automatica. 2014. Vol. 50(12). 3204-3208.
3. Kharitonov V.L., Zhabko A.P. Lyapunov-Krasovskii approach to the robust
stability analysis of time-delay systems // Automatica. 2003. Vol. 39(1). P.
15-20.
4. Infante, E.F., Castelan, W. V. A Lyapunov functional for a matrix differencedifferential equation // J. Differ. Equat. 29, 439-451. 1978.
5. Huesca E., Mondie S., Santos O. Polynomial approximations of the Lyapunov
matrix of a class of time delay systems // Sinaia, Romania. 2009.
6. Jarlebring E., Vanbiervliet J., Michiels W. Characterizing and computing
H2 norm of time delay systems by solving the delay Lyapunov Equaion //
Proceedings of the 49th IEEE Conference on Decision and Control. 2010.
7. Хатсон В., Пим Дж. Приложения функционального анализа и теории операторов // М., Мир. с. 215-223. 1983.
8. Беллман Р., Кук К. Дифференциально-разностные уравнения // М., Мир.
1967.
9. Чижова О. Н. Методы исследования дифференциальных уравнений с запздывающим аргументом // СПб, Соло. с. 5-11. 2011.
24
10. Колмогоров А. Н., Фомин. С. В. Элементы теории функций и функционального анализа // 4-е изд. М.: Наука. с. 544. 1976.
11. Крылов В. И., Бобков В. В., Монастырный П. И. Вычислительные методы
Том 2 // М., Наука. с. 267-273. 1977.
25
Отзывы:
Авторизуйтесь, чтобы оставить отзыв