Алгоритм обчислення власних значень

Однією з найважливіших задач обчислювальної математики є створення ефективних і стійких алгоритмів знаходження власних значень матриці. Ці алгоритми обчислення власних значень можуть також знаходити власні вектори.

Власні значення і власні векториРедагувати

Якщо задано n × n квадратну матрицю A над дійсними або комплексними числами, власне значення λ і відповідний йому кореневий вектор v - це пара, що задовольняє рівності[1]

 

де v ненульовий n × 1 вектор-стовпець, E — n × n одинична матриця, k — додатне ціле, а λ і v можуть бути комплексними, навіть якщо A дійсне. Якщо k = 1, вектор просто називається власним вектором. У цьому випадку Av = λv. Будь-яке власне значення λ матриці A має простий[прим. 1] власний вектор, відповідний йому так, що якщо k — найменше ціле, за якого (A - λE)k v = 0 для кореневого вектора v, то (A - λE)k-1 v буде простим власним вектором. Значення k завжди можна взяти меншим або рівним n. Зокрема, (A - λE)n v = 0 для всіх кореневих векторів v відповідних λ.

Для будь-якого власного значення λ матриці A ядро ker(A - λE) складається з усіх власних векторів, відповідних λ, (разом з 0) і називається власним підпростором числа λ, а векторний підпростір ker((A - λE)n) складається з усіх кореневих векторів (доповнений нульовим вектором) і називається кореневим підпростором. Геометрична кратність значення λ є розмірністю його власного підпростору. Алгебрична кратність значення λ є розмірністю його кореневого підпростору. Подальші терміни пов'язані з рівністю

 

Тут det — визначник, λi — всі різні власні значення матриці A, а αi — відповідні алгебричні кратності. Функція pA(z) — це характеристичний многочлен матриці A. Отже, алгебрична кратність є кратністю власних значень як коренів характеристичного многочлена. Оскільки будь-який власний вектор є кореневим вектором, геометрична кратність менша або дорівнює алгебричній кратності. Сума алгебричних кратностей дорівнює n степеню характеристичного многочлена. Рівняння pA(z) = 0 називають характеристичним рівнянням, оскільки його корені є точно власними значеннями матриці A. В теоремі Гамільтона — Келі сама матриця A задовольняє тому самому рівнянню: pA(A) = 0[прим. 2]. Як наслідок, стовпці матриці   мають бути або 0, або кореневими векторами, відповідними власному значенню λj, оскільки вони знищуються матрицею  

Будь-який набір кореневих векторів різних власних значень лінійно незалежний, так що базис для всього C n можна вибрати з набору кореневих векторів. Точніше цей базис {vi}ni=1 можна вибрати і вибудувати так, що

  • якщо vi і vj мають одне і те саме власне значення, то це виконуватиметься для будь-якого vk при k між i і j;
  • якщо vi не є простим власним вектором і якщо λi — його власне значення, то (A - λiE )vi = vi-1 (зокрема v1 має бути простим власним вектором).

Якщо ці базисні вектори розташувати як стовпці матриці V = [ v1 v2vn ], то V можна використовувати для перетворення матриці A в її нормальну жорданову форму:

 

де λi — власне значення, βi = 1 якщо (A - λi+1)vi+1 = vi і βi = 0 в інших випадках.

Якщо W є оборотною матрицею і λ — власне значення матриці A з відповідним кореневим вектором v, то (W -1AW - λE )k W -kv = 0. Отже, λ є власним значенням матриці W -1AW з відповідним кореневим вектором W -kv. Таким чином, подібні матриці мають однакові власні значення.

Нормальні, ермітові і дійсні симетричні матриціРедагувати

Ермітово-спряжена матриця M* до комплексної матриці M — це траспонована матриця із заміною всіх елементів спряженими значеннями: M * = M T. Квадратну матрицю A називають нормальною, якщо вона комутує з ермітово-спряженою: A*A = AA*. Матрицю називають ермітовою, якщо вона дорівнює своїй спряженій: A* = A. Всі ермітові матриці нормальні. Якщо A має тільки дійсні елементи, то спряжена до неї — це просто транспонована матриця і вона буде ермітовою в тому і тільки в тому випадку, коли вона симетрична. Якщо застосувати це до стовпців, спряженість можна використовувати для визначення канонічного скалярного добутку в C n: w • v = w* v[прим. 3]. Нормальні, ермітові і дійсні симетричні матриці мають низку корисних властивостей:

  • Кожен кореневий власний вектор нормальної матриці є простим власним вектором.
  • Будь-яка нормальна матриця подібна до діагональної, оскільки її нормальна жорданова форма є діагональною матрицею.
  • Власні вектори, відповідні різним власним значенням нормальної матриці, ортогональні.
  • Для будь-якої нормальної матриці A C n має ортонормальний базис, що складається з власних векторів матриці A. Відповідна матриця власних векторів є унітарною.
  • Власні значення ермітової матриці є дійвсними числами, оскільки (λ — λ)v = (A* — A)v = (A — A)v = 0 для ненульового власного вектора v.
  • Якщо матриця A дійсна, існує ортонормальний базис для R n, що складається з власних векторів матриці A, в тому і тільки в тому випадку, коли A симетрична.

Як дійсні, так і комплексні матриці, які не є при цьому ермітовими матрицями, можуть мати всі власні значення дійсними. Наприклад, дійсна трикутна матриця має всі свої власні значення на діагоналі, але, в загальному випадку, не симетрична.

Число обумовленостіРедагувати

Будь-яку задачу обчислювальної математики можна розглядати як обчислення деякої функції ƒ від деякого аргументу x. Число обумовленості κ(ƒ, x) задачі — це відношення відносної похибки результату обчислення до відносної похибки параметра функції і залежить як від функції, так і від параметра. Число обумовленості описує, наскільки зростає похибка під час обчислень. Десятковий логарифм цього числа говорить про кількість знаків, які ми втрачаємо відносно вихідних даних. Число обумовленості стосується найкращого сценарію, відбиваючи нестабільність самої задачі, незалежно від способу розв'язування. Ніякий алгоритм не може дати результату кращого, ніж визначений числом обумовленості, хіба що випадково. Однак погано розроблений алгоритм може дати істотно гірші результати. Наприклад, як буде згадано нижче, задача знаходження власних значень нормальної матриці завжди добре обумовлена, проте задача знаходження коренів многочлена може бути дуже погано обумовленою[en]. Такі алгоритми обчислення власних значень, які працюють шляхом знаходження коренів характеристичного многочлена, можуть виявитися погано обумовленими, навіть якщо сама задача добре обумовлена.

Для задачі розв'язування системи лінійних рівнянь Av = b, де A є оборотною, число обумовленості κ(A-1, b) визначається виразом ||A||op||A−1||op, де || ||op — операторна норма, підпорядкована звичайній евклідовій нормі на C n. Оскільки це число не залежить від b і є тим самим як для A, так і для A-1, його зазвичай називають числом обумовленості κ(A) матриці A. Це значення κ(A) є також абсолютним значенням відношення найбільшого власного значення матриці A до її найменшого власного значення. Якщо A є унітарною, то ||A||op = ||A−1||op = 1, так що κ(A) = 1. У загальному випадку для матриць часто складно обчислити операторну норму. З цієї причини для оцінення числа обумовленості зазвичай використовують інші норми матриці.

Для задачі обчислення власних значень Бауер і Файк довели[en], що якщо λ є власним значенням діагоналізовної n × n матриці A з матрицею власних векторів V, то абсолютна похибка обчислення λ обмежена добутом κ(V) і абсолютною похибкою в A:  [2]. Як наслідок, число обумовленості для обчислення λ дорівнює κ(λ, A) = κ(V) = ||V ||op ||V −1||op. Якщо матриця A нормальна, то V є унітарною і κ(λ, A) = 1. Таким чином, задача обчислення власних значень нормальних матриць добре обумовлена.

Показано, що число обумовленості задачі обчислення власного підпростору нормальної матриці A, відповідного власного значення λ, обернено пропорційне мінімальній відстані між λ та іншими, відмінними від λ, власними значеннями матриці A[3]. Зокрема, задача визначення власного підпростору для нормальних матриць добре обумовлена для ізольованих власних значень. Якщо власні значення не ізольовані, найкраще, на що можна розраховувати, це визначення підпростору всіх власних векторів довколишніх власних значень.

АлгоритмиРедагувати

Будь-який нормований многочлен є характеристичним многочленом супутньої матриці, тому алгоритм для обчислення власних значень можна використовувати для знаходження коренів многочленів. Теорема Абеля — Руффіні показує, що будь-який такий алгоритм для розмірності більшої від 4 має або бути нескінченним, або залучати функції складніші, ніж елементарні арифметичні операції або дробові степені. З цієї причини алгоритми, які обчислюють точно власні значення за скінченне число кроків, існують тільки для спеціальних класів матриць. У загальному випадку алгоритми є ітеративними і дають на кожній ітерації чергове наближення до розв'язку.

Деякі алгоритми дають усі власні значення, інші дають кілька значень або навіть всього одне, однак і ці алгоритми можна використовувати для обчислення всіх власних значень. Як тільки власне значення λ матриці A визначено, його можна використовувати або для зведення алгоритму до отримання іншого власного значення, або для зведення задачі до такої, для якої λ не є розв'язком.

Зведення зазвичай виконують зсувом: A замінюється на A - μE для деякої константи μ. Щоб отримати власне значення матриці A, потрібно власне значення, знайдене для A - μE, додати до μ. Наприклад, у степеневому методі μ = λ. Ітерація степеневого методу знаходить найбільше за абсолютною величиною значення, так що навіть якщо λ є наближенням до власного значення, ітерація степеневого методу навряд чи знайде його вдруге. І навпаки, методи, засновані на зворотних ітераціях знаходять найменше власне значення, так що μ вибирається подалі від λ в надії опинитися ближче до якогось іншого власного значення.

Зведення можна виконати звуженням матриці A до простору стовпців матриці A - λE. Оскільки A - λE вироджена, простір стовпців має меншу розмірність. Алгоритм обчислення власних значень можна тоді застосувати до цієї звуженої матриці. Процес можна продовжувати, поки не буде знайдено всі власні значення.

Якщо алгоритм знаходження власних значень не дає власного вектора, поширеним є застосування алгоритму, заснованого на зворотній ітерації, з прирівнюванням μ до найближчої апроксимації власного значення. Це швидко приводить до власного вектора найближчого до μ власного значення. Для невеликих матриць альтернативою є використання стовпцевого підпростору добутку A - λ́E для кожного з решти власних значень λ́.

Матриці Гессенберга і тридіагональні матриціРедагувати

Оскільки власними значеннями трикутної матриці є діагональні елементи, в загальному випадку не існує скінченного методу, подібного до виключень Гауса, для зведення матриці до трикутної форми, зберігаючи при цьому власні значення, однак можна отримати щось близьке до трикутної матриці. Верхня матриця Гессенберга — це квадратна матриця, в якій усі елементи нижче від першої піддіагоналі дорівнюють нулю. Нижня матриця Гессенберга — це квадратна матриця, в якій усі члени вище від першої наддіагоналі дорівнюють нулю. Матриці, які є як нижніми, так і верхніми матрицями Гессенберга — це тридіагональні матриці. Матриці Гессенберга і тридіагональні матриці є початковими для багатьох алгоритмів обчислення власних значень, оскільки нульові значення зменшують складність задачі. Існує кілька методів зведення матриць до матриць Гессенберга з тими ж власними значеннями. Якщо початкова матриця симетрична або ермітова, то кінцева матриця буде тридіагональною.

Якщо потрібні тільки власні значення, не потрібно обчислювати матрицю подібності, оскільки перетворена матриця має ті самі власні значення. Якщо також потрібні і власні вектори, матриця подібності необхідна для перетворення власних векторів матриці Гессенберга на власні вектори початкової матриці.

Метод Застосовний до матриць Результат Ціна без матриці подібності Ціна з матрицею подібності Опис
Перетворення Хаусхолдера загального вигляду матриця Гессенберга 2n33 + O(n2)[4] 4n33 + O(n2)[4] Відбиття кожного стовпця відносно підпростору для обнулення нижніх елементів стовпця
Повороти Ґівенса загального вигляду матриця Гессенберга 4n33 + O(n2)[4] Здійснюється плоске обертання для обнулення окремих елементів. Обертання впорядковані так, що наступні обертання не зачіпають нульових елементів
Ітерації Арнольді загального вигляду матриця Гессенберга Здійснюється ортогоналізація Грама — Шмідта на підпросторах Крилова
Алгоритм Ланцоша[en][5] ермітова тридіагональна матриця Ітерації Арнольді для ермітових матриць

Ітеративні алгоритмиРедагувати

Ітеративні алгоритми розв'язують задачу обчислення власних значень побудовою послідовностей, що збігаються до власних значень. Деякі алгоритми дають також послідовності векторів, що збігаються до власних векторів. Найчастіше послідовності власних значень виражаються через послідовності подібних матриць, які збігаються до трикутної або діагональної форми, що дозволяє потім просто отримати власні значення. Послідовності власних векторів виражаються через відповідні матриці подібності.

Метод Застосовний до матриць Результат Ціна за один крок Збіжність Опис
Степеневий метод загального вигляду найбільше власне значення і відповідний вектор O(n2) лінійна Багаторазове множення матриці на довільно вибраний початковий вектор з подальшою нормалізацією
Зворотний степеневий метод загального вигляду найближче до μ власне значення і відповідний вектор лінійна Степенева ітерація з матрицею (A - μE )-1
Метод ітерацій Релея загального вигляду найближче до μ власне значення і відповідний вектор кубічна Степенева ітерація з матрицею (A - μiE )-1, де μi — відношення Релея від попередньої ітерації
Передобумовлена зворотна ітерація[6] або LOBPCG[en] додатноозначена дійсна симетрична найближче до μ власне значення і відповідний вектор Зворотна ітерація з передобумовленням (наближене звернення матриці A)
Метод поділу навпіл[7] дійсна симетрична тридіагональна будь-яке власне значення лінійна Використовує метод бісекції для пошуку коренів характеристичного многочлена і властивості послідовності Штурма
Ітерації Лагерра дійсна симетрична тридіагональна будь-яке власне значення кубічна[8] Використовує метод Лагерра[en] обчислення коренів характеристичного многочлена і властивості послідовності Штурма
QR-алгоритм[9] Гессенберга всі власні значення O(n2) кубічна Розкладання A = QR, де Q ортогональна, R — трикутна, потім використовується ітерація до RQ
всі власні значення 6n3 + O(n2)
Метод Якобі дійсна симетрична всі власні значення O(n3) квадратична Використовує поворот Ґівенса у спробі позбутися від недіагональних елементів. Спроба не вдається, але підсилює діагональ
Розділяй і владарюй[en] ермітова тридіагональна всі власні значення O(n2) Матриця розбивається на підматриці, які діагоналізуються, потім возз'єднуються
всі власні значення (43)n3 + O(n2)
Метод гомотопії дійсна симетрична тридіагональна всі власні значення O(n2)[10] Будується обчислювана гомотопія
Метод спектральної згортки[en] дійсна симетрична найближче до μ власне значення і відповідний власний вектор Передобумовлена зворотна ітерація, застосована до (A - μE )2
Алгоритм MRRR[11] дійсна симетрична тридіагональна деякі або всі власні значення і відповідні власні вектори O(n2) «Multiple Relatively Robust Representations» — здійснюється зворотна ітерація з розкладанням LDLT зміщеної матриці.

Пряме обчисленняРедагувати

Не існує простих алгоритмів прямого обчислення власних значень для матриць у загальному випадку, але для багатьох спеціальних класів матриць власні значення можна обчислити прямо. Це:

Трикутні матриціРедагувати

Оскільки визначник трикутної матриці є добутком її діагональних елементів, то для трикутної матриці T  . Отже, власні значення T — це її діагональні елементи.

Розкладані поліноміальні рівнянняРедагувати

Якщо p — будь-який многочлен і p(A) = 0, то власні значення матриці A задовольняють тому ж рівнянню. Якщо вдається розкласти многочлен p на множники, то власні значення A містяться серед його коренів.

Наприклад, проєктор — це квадратна матриця P, що задовольняє рівнянню P2 = P. Коренями відповідного скалярного поліноміального рівняння λ2 = λ будуть 0 і 1. Таким чином, проєктор має 0 і 1 власними значеннями. Кратність власного значення 0 — це дефект P, тоді як кратність 1 — це ранг P.

Інший приклад — матриця A, що задовольняє рівнянню A2 = α2E для деякого скаляра α. Власні значення мають бути рівними ±α. Оператори проєктування

 
 

задовольняють рівностям

 

і

 

Простори стовпців матриць P+ і P- є підпросторами власних векторів матриці A, які відповідають і , відповідно.

Матриці 2×2Редагувати

Для розмірностей від 2 до 4 існують формули з використанням радикалів, які можна використовувати для обчислення власних значень. Для матриць 2х2 і 3х3 це звичайна практика, але для матриць 4х4 через зростання складності формул коренів[en] цей підхід менш привабливий.

Для матриць 2×2

 

характеристичний многочлен дорівнює

 

Власні значення можна знайти як корені квадратного рівняння:

 

Якщо визначити   як відстань між двома власними значеннями, легко обчислити

 

з подібними формулами для c і d. Звідси випливає, що обчислення добре обумовлене, якщо власні значення ізольовані.

Власні вектори можна отримати, використовуючи теорему Гамільтона — Келі. Якщо λ1, λ2 — власні значення, то (A - λ1E )(A - λ2E ) = (A - λ2E )(A - λ1E ) = 0, так що стовпці (A - λ2E ) обнуляються матрицею (A - λ1E ) і навпаки. Припускаючи, що жодна з матриць не дорівнює нулю, стовпці кожної матриці мають містити власні вектори для іншого власного значення (якщо ж матриця нульова, то A є добутком одиничної матриці на константу і будь-який ненульовий вектор є власним).

Наприклад, нехай

 

тоді tr(A) = 4 - 3 = 1 і det(A) = 4(-3) - 3(-2) = -6, так що характеристичне рівняння

 

а власні значення рівні 3 і -2. Тепер

 ,  

В обох матрицях стовпці відрізняються скалярними коефіцієнтами, так що можна вибирати будь-який стовпець. Так, (1, -2) можна використовувати як власний вектор, відповідного власному значенню -2, а (3, -1) — як власний вектор для власного значення 3, що легко можна перевірити множенням на матрицю A.

Матриці 3 х 3Редагувати

Якщо A — матриця 3 х 3, то характеристичним рівнянням буде:

 

Це рівняння можна розв'язати за допомогою методів Кардано або Лагранжа, але афінне перетворення матриці A істотно спрощує вираз і веде прямо до тригонометричного розв'язку. Якщо A = pB + qE, то A і B мають однакові власні вектори і β є власним значенням матриці B тоді й лише тоді, коли α = + q є власним значенням матриці A. Якщо покласти   і  , одержимо

 

Підстановка β = 2cos θ і деяке спрощення з використанням тотожності cos 3θ = 4cos3 θ - 3cos θ зводить рівняння до cos 3θ = det(B) / 2. Отже,

 

Якщо det(B) є комплексним числом або більше 2 за абсолютною величиною, арккосинус для всіх трьох величин k слід брати з однієї гілки. Проблема не виникає, якщо A дійсна і симетрична, приводячи до простого алгоритму:[12]

% Given a real symmetric 3x3 matrix A, compute the eigenvalues

p1 = A(1,2)^2 + A(1,3)^2 + A(2,3)^2
if (p1 == 0) 
  % A is diagonal.
  eig1 = A(1,1)
  eig2 = A(2,2)
  eig3 = A(3,3)
else
  q = trace(A)/3
  p2 = (A(1,1) - q)^2 + (A(2,2) - q)^2 + (A(3,3) - q)^2 + 2 * p1
  p = sqrt(p2 / 6)
  B = (1 / p) * (A - q * E)    % E is the identity matrix
  r = det(B) / 2

  % In exact arithmetic for a symmetric matrix -1 <= r <= 1
  % but computation error can leave it slightly outside this range.
  if (r <= -1) 
   phi = pi / 3
  elseif (r >= 1)
   phi = 0
  else
   phi = acos(r) / 3
  end

  % the eigenvalues satisfy eig3 <= eig2 <= eig1
  eig1 = q + 2 * p * cos(phi)
  eig3 = q + 2 * p * cos(phi + (2*pi/3))
  eig2 = 3 * q - eig1 - eig3   % since trace(A) = eig1 + eig2 + eig3
end

Власні вектори також A можна отримати, скориставшись теоремою Гамільтона — Келі. Якщо α1, α2, α3 — різні власні значення матриці A, то (A - α1E)(A - α2E)(A - α3E) = 0. Тоді стовпці добутку будь-яких двох з цих матриць містять власні вектори третього власного значення. Однак якщо a3 = a1, то (A - α1E)2(A - α2E) = 0 і (A - α2E)(A - α1E)2 = 0. Таким чином, кореневий власний підпростір α1 натягнуий на стовпці A - α2E, тоді як звичайний власний підпростір натягнутий на стовпці (A - α1E)(A - α2E). Звичайний власний підпростір α2 натягнутий на стовпці (A - α1E)2.

Наприклад, нехай

 

Характеристичне рівняння буде

 

з власними значеннями 1 (кратності 2) і −1. Обчислюємо

 ,

а потім

 .

Тоді (-4, -4, 4) є власним вектором для −1, а (4, 2, -2) є власним вектором для 1. Вектори (2, 3, -1) і (6, 5, -3) є кореневими векторами, відповідними значенню 1, будь-який з яких можна скомбінувати з (-4, -4, 4) і (4, 2, -2), утворюючи базис кореневих векторів матриці A.

Див. такожРедагувати

КоментаріРедагувати

  1. Термін «простий» тут лиш підкреслює відмінність між «власним вектором» і «кореневим вектором».
  2. де сталий член множиться на одиничну матрицю E.
  3. Такому порядку в скалярному добутку (зі спряженим елементом ліворуч) надають перевагу фізики. Алгебристи часто застосовують запис w • v = v* w.

ПриміткиРедагувати

  1. Sheldon Axler. Down with Determinants! // American Mathematical Monthly. — 1995. — Вип. 102 (19 вересня). — С. 139—154. Архівовано з джерела 13 вересня 2012. Процитовано 11 вересня 2021.
  2. F. L. Bauer, C. T. Fike. Norms and exclusion theorems // Numer. Math. — 1960. — Вип. 2 (19 вересня). — С. 137—141.
  3. S. C. Eisenstat, I. C. F. Ipsen. Relative Perturbation Results for Eigenvalues and Eigenvectors of Diagonalisable Matrices // BIT. — 1998. — Т. 38, вип. 3 (19 вересня). — С. 502—9. — DOI:10.1007/bf02510256.
  4. а б в William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery. Numerical Recipes in C. — 2nd. — Cambridge University Press, 1992. — ISBN 0-521-43108-5.
  5. Х. Д. Икрамов. Разреженные матрицы. — 1982. — Т. 20. — (Итоги науки и техники. Сер. Мат. анал)
  6. K. Neymeyr. A geometric theory for preconditioned inverse iteration IV: On the fastest convergence cases. // Linear Algebra Appl. — 2006. — Т. 415, вип. 1 (19 вересня). — С. 114—139. — DOI:10.1016/j.laa.2005.06.022.
  7. (Уилкинсон, 1970), стр. 274, Метод деления пополам
  8. T. Y Li, Zhonggang Zeng. Laguerre's Iteration In Solving The Symmetric Tridiagonal Eigenproblem - Revisited // SIAM Journal on Scientific Computing. — 1992. — 19 вересня.
  9. (Парлетт, 1983), стр. 156, глава 8, Алгоритмы QR и QL
  10. Moody T. Chu. A Note on the Homotopy Method for Linear Algebraic Eigenvalue Problems // Linear Algebra Appl. — 1988. — Т. 105. — С. 225—236. — DOI:10.1016/0024-3795(88)90015-8.
  11. Inderjit S. Dhillon, Beresford N. Parlett, Christof Vömel. The Design and Implementation of the MRRR Algorithm // ACM Transactions on Mathematical Software. — 2006. — Т. 32, вип. 4 (19 вересня). — С. 533—560. — DOI:10.1145/1186785.1186788.
  12. Oliver K. Smith. Eigenvalues of a symmetric 3 × 3 matrix // Communications of the ACM. — Т. 4, вип. 4. — С. 168. — DOI:10.1145/355578.366316.

ЛітератураРедагувати

  • Дж. Голуб, Ч. Ван Лоун. Матричные вычисления. — М. : «Мир», 1999. — ISBN 5-03-002406-9.
  • Б. Парлетт. Симметричная проблема собственных значений. — М. : «Мир», 1983.
  • Дж. Х. Уилкинсон. Алгебраическая проблема собственных значений. — М. : «Наука» Главная редакция физико-математической литературы, 1970.

Додаткова літератураРедагувати