#import "utils.typ": *

#set page("a5", numbering: "1")
#set text(lang: "ru")
#set par(justify: true)
#set heading(numbering: "1.")


`24 Сентября 2024`
= Вводная

<NAME>, D947

Чтобы получить экзамен нужно:
- допуск от преподавателя по практике
- сдача теории (экзамен)  

Можно сдавать частями (три коллоквиума):
- Прямые методы
- Итерационные методы
- Собственные значения

Курсовой проект. Подробно разобрать конкретный метод, который каждому будет дан. Ещё будет две теоретические задачи.

#linebreak()

= Первая лекция

В конечномерном пространстве все нормы эквивалентны. 

#def[
  $norm(x)$ --- число, для которого выполняются 3 аксиомы:\
  1. $norm(x) >= 0, norm(x) = 0 <=> x = 0$
  2. $norm(lambda x) = lambda norm(x)$
  3. $norm(x+y) <= norm(x) + norm(y)$
]

#def[
  Нормы вектора.
  $
    &norm(x)_1 = sum_(i=1)^n abs(x_i)\
    &norm(x)_2 = sqrt((x, x)) - "евклидова норма"\
    &norm(x)_oo = max_i abs(x_i)
  $
]

СЛАУ:
$ A x = f $

#def[
  Нормы матрицы.

  Евклидова норма:

  $ norm(A)_E = sqrt(sum_(i, j) a^2_(i j)) $
]

#def[
  Подчинённая матричная норма:
  $ norm(А) = sup_(x!=0) norm(A x) / norm(x) $
]

Свойства:
1. $norm(I) = 1 $
2. $norm(A B) <= norm(A) norm(B) $

Евклидова норма не является подчинённой, ведь
$norm(I)_E = sqrt(n) != 1$.

Из определения подчинённой матричной нормы вытекает:

$ norm(А) = sup_(x!=0) norm(A x) / norm(x) => norm(A x) <= norm(A) norm(x) $

Обратное 

#def[
  Согласованная матричная норма такая, что:
  $ norm(A x) <= norm(A) norm(x) $
]

$
  &norm(A)_1 = max_j sum_(i=1)^n abs(a_(i j)) \
  &norm(A)_2 = sqrt(max_i lambda_i (A^* A)), " где" A^* - "транспонированная матрица"\
  &norm(A)_oo = max_i sum_(j=1)^n abs(a_(i j))
$

$
  norm(A) = sup_(x!=0) norm(A x) / norm(x)
$

// $ norm(A x)_1 = max_i abs(sum_(j=1)^n a_(i j)x_j)
// <= max_i sum_(j=1)^n abs(a_(i j)) abs(x_j) $

$ norm(A x)_1 =
  sum_(i=1)^n abs(sum_(j=1)^n a_(i j)x_j) <=
  sum_(i=1)^n sum_(j=1)^n abs(a_(i j)) abs(x_j) = \
  = sum_(j=1)^n sum_(i=1)^n abs(a_(i j)) abs(x_j) =
  sum_(j=1)^n (abs(x_j) sum_(i=1)^n abs(a_(i j))) <= \
  <= sum_(j=1)^n ((max_j sum_(i=1)^n abs(a_(i j)))abs(x_j))  =
  max_j sum_i abs(a_(i j)) dot sum_j abs(x_j) =\
  = (max_j sum_i abs(a_(i j))) norm(x)_1
  
$

Доказали, что

$ norm(A x)_1 / norm(x)_1 <= max_j sum_(i=1)^n abs(a_(i j)) $

$k -$  номер столбца, где достигается максимум:

$ sum_(i=1)^n abs(a_(i k)) = max_j sum_(i=1)^n abs(a_(i j)) $

$ e_k = vec(0, dots, underbrace(1, k-"й индекс"), dots, 0, 0) 
  = max_j sum_(i=1)^n abs(a_(i j)) norm(x) ???
$

$ norm(A e_k)_1 = sum_i abs(a_(i k)) $

Для второй матричной нормы  $norm(A)_2 = sqrt(max_i lambda_i (A^* A))$:

$ lambda_i$ должны быть действительными и положительными

$A^* A - $ неотрицательно определена и симметрична

$ (A^* A)^* = A^* (A^*)^* = A^* A $

#def[
  Матрица $C$ называется неотрицательно определённой, если
  $ forall (y, y) >= 0 $
  Похожим образом определяются положительно определённая матрица и другие.
]

Свойство скалярного произведения:

$ (B x, y) = (x, B^* y) $

$ (A^*A x, x) = (A x, A x)) >= 0 $

Из курса линейной алгебры известно, что все собственные значения симметрической матрицы вещественны.

Значит, $norm(A)_2 = sqrt(max_i lambda_i (A^* A))$ имеет смысл для симметрических матриц.

У симметрических матриц существует полная система ортонормированных собственных векторов.

$ (A^* A) u^((i)) = lambda_i u^((i))\
(u^((i)), u^((j))) = cases(
  1 when i = j,
  0 when i != j
)
$

$ x = sum_(i=1)^n alpha_i u^((i)) $

$ norm(x)_2 = sqrt((x, x)) = sqrt((sum_(i=1)^n alpha_i u^((i)), sum_(j=1)^n alpha_j u^((j)))) = \
= sqrt(sum_(i=1)^n alpha_i^2)
$

$
  norm(A x)_2 = sqrt((A x, A x)) =
  sqrt((A^* A x, x)) =\
  = sqrt((A^*A sum_i alpha_i u^((i)), sum_j alpha_j u^((j)) ))
  = sqrt((sum_i alpha_i A^* A u^((i)) comma sum_j alpha_j u^((j)))) =\
  = sqrt((sum_i alpha_i lambda_i u^((i)) comma sum_j alpha_j u^((j))))    
  = sqrt(sum_i alpha_i^2 lambda_i) <=\
  <= sqrt(sum_i alpha_i^2 max_j lambda_i (A^* A))
  = sqrt(max_i lambda_i (A^*A)) norm(x)_2
$

$ lambda_k (A^* A) = max_i lambda_i (A^* A) $

Посчитаем вторую норму. Предположим, что матрица симметричная:
$ A^* = A $

$ norm(A)_2 = sqrt(max_i lambda_i (A^* A)) =
  sqrt(max_i lambda_i (A^2)) =
  sqrt(max_i lambda_i^2 (A)) = 
  max_i abs(lambda_i (A))
$

#linebreak()

`1 Октября 2024`
= Обусловленность матрицы систем

#def[
  Число обусловленности:
$ mu(A) = sup_(x, y != 0)(norm(A x)/norm(x) div norm(A y)/norm(y)) $
]

$ A x = f $

$ A(x+ xi) = f + phi $

$ A xi = phi $

Попробуем оценить следующую величину:

$ norm(xi) / norm(x) dot norm(f) / norm(phi) = 
  norm(xi) / norm(x) dot norm(A x) / norm( A xi) = 
  norm(A x) / norm(x) div norm(A xi) / norm(xi) <=
  mu(A)
$

$ norm(xi) / norm(x) <= mu(a) norm(phi) / norm(f) $

Представим, что $mu(A)$ невелико, меньше единицы оно быть не может, тогда маленькое возмущение правой части гарантирует, что возмущение в решении тоже невелико.

Как влияет погрешность матрицы на результат:

$ (A + Sigma)(x + xi) = f + phi $

$ norm(xi)/norm(x) <= mu(A) / (1 - mu(A) norm(Sigma)/norm(A)) (norm(phi)/norm(f) + norm(Sigma)/norm(A)) $

При условии:

$ norm(A^(-1)) norm(Sigma) < 1 $

Нам нужна формула, по которой мы будем искать число обусловленности, выведем её.

$ mu(A) = sup_(x!=0)(norm(A x)/norm(x)) / (inf_(y!=0)(norm(A y)/norm(y))) 
  = norm(A) / (inf_(z!=0) 1/(norm(A^(-1) z) / norm(z))) =\
  =norm(A) / (1/(sup_(z!=0) norm(A^(-1) z) / norm(z)))
  = norm(A) norm(A^(-1))
$

$ z := A y $
#showybox(breakable: true, title: "Пример")[
  $ A = mat(1, 1.0001; 1.0001, 1) wide f = vec(2.0001, 2.0001) $

  $ x = vec(1, 1) $

  $ norm(A) = max_(i = 1,2) abs(lambda_i (A)) $

  $ matdet(1-lambda, 1.0001; 1.0001, 1-lambda) = 0 $

  $ (1-lambda)^2 - 1.0001^2 = 0 $

  $ lambda_1 = 0.0001 wide lambda_2 = -10^(-4)  $

  $ norm(A) = 2.0001 wide norm(A^(-1)) = 10^4 $

  $ mu(A) = 2001 "большое число" $

  Рассмотрим следующую систему:

  $ (A+Sigma)y = f + phi $
  
  $ Sigma = mat(0, 0; 0, 0.00005) wide
    phi = vec(0.0001, 0)  
  $
  
  Решение системы:
  
  $
    y = vec(0, 2)
  $
]

#linebreak()
= Прямые методы

== Метод Гаусса (метод исключения неизвестных)

Метод в основе которого лежит использование элементарных преобразований матрицы с целью применения её к треугольному или диагональному виду.

$ A x = f $

$ cases(
    a_(11)x_1 + ... + a_(1 n)x_n = f_1,
    dots,
    a_(n 1)x_1 + ... + a_(n n)x_n = f_n,
  )
$

Исключаем элемент $x_1$. Пусть $a_11 != 0$:

$ x_1 + a_(12) / a_11 x_2  + ... + a_(1 n)/a_(11) x_n = f_1 / a_11 | dot a_21 $

$ a_21 x_1 + a_21 a_(12) / a_11 x_2  + ... + a_21 a_(1 n)/a_(11) x_n = a_21 f_1 / a_11 $

Теперь из второго уравнения системы вычитаем вот это, таким образом получим уравнение, в котором нет $x_1$.

$ cases(
    a_(11)x_1 + a_(22) x_2 +... + a_(1 n)x_n = f_1,
    #hide[$a_(11)x_1 + $] a_22^((1))x_2 + ... + a_(2 n)^((1))x_n = f_2^((1)),
    dots,
    a_(n 1)x_1 + a_(n 2)x_2 +... + a_(n n)x_n = f_n,
  )
$

$ cases(
    a_(11)x_1 + a_(22) x_2 +... + a_(1 n)x_n = f_1,
    #hide[$a_(11)x_1 + $] a_22^((1))x_2 + ... + a_(2 n)^((1))x_n = f_2^((1)),
    dots,
    #hide[$a_(11)x_1 + $] a_(n 2)^((1))x_2 +... + a_(n n)^((1))x_n = f_n^((1)),
  )
$

Второй шаг. $a_22^((1)) != 0$ и делаем аналогично.

После $n-1$ шагов получим матрицу, которая имеет верхний треугольный вид.

Теперь запишем формулы всех этих шагов.

1 шаг:

$ d_(i 1) = a_(i 1) / a_11 $

$ a_(i j)^((1)) = a_(i j) - d_(i 1) a_(1 j), quad i, j = 2, ..., n $

$ f_i^((1)) = f_i - d_(i 1) f_1 $

2 шаг:

$ d_(i 2) = a_(i 2)^((1)) / a_22^((1)) $

$ a_(i j)^((2)) = a_(i j)^((1)) - d_(i 2) a_(2 j)^((1)), quad i, j = 3, ..., n $

$ f_i^((2)) = f_i^((1)) - d_(i 2) f_2^((1)) $

$ A_(n-1) x = f^((n -1)) $
$ A_(n-1) =: U "верхняя треугольная матрица" $

$ f^((n-1)) =: g $

$ U x = g $

$ u_(n n) x_n = g_n $

$ x_n = g_n / u_(n n) $

$k-й$ шаг обратного хода:

$ u_(k k) x_k + sum_(j=k+1)^n u_(k j) x_j = g_k $

$ x_k = (g_k - sum_(j=k+1)^n u_(k j)x_j) / u_(k k) $

Условие, при котором ведущий элемент не равен нулю. Во-первых матрица должна быть невырожденной.

#bbox[Критерий работы метода Гаусса][
Отличие от нуля всех главных миноров.]

$ det mat(
  a_11, a_12, ..., a_(1 k);
  a_21, a_22, ..., a_(2 k);
  ..., ..., ..., ... ;
  a_(k 1), a_(k 2), ..., a_(k k)
) != 0 $

$ a_11 != 0 $

$ matdet(a_11, a_12; 0, a_22^((1))) != 0 $

$ a_11 a_22^((1)) != 0 $

$ a_22^((1)) != 0 $

$ matdet(
  a_11, a_12, a_13;
  0, a_22^((1)), a_23^((1));
  0, 0, a_33^((2))
) != 0 $

$ a_11 a_22^((1)) a_33^((2)) != 0 $

$ a_33^((2)) != 0 $

Все главные элементы отличны от нуля.

== LU-разложене

Тот же метод Гаусса, но записанный через разложение матриц.

$ D_1 := mat(
  1, 0, 0;
  -d_21, 1, 0;
  -d_31, 0, 1;
  dots.v;
  -d_(n 1), 0, 0, dots, 1
) $

$ A_1 = D_1 A $

$ D_2 := mat(
  1;
  0, 1, ..., 0;
  0, -d_32, ...;
  0, -d_42, ...;
  dots.v;
  0, -d_(n 2), ..., 1;
  
) $

$ A_2 = D_2 A_1 $

$ A_(n -1) $ 

$ U = A_(n - 1) = D_(n-1) A_(n-2) = D_(n-1) D_(n-2) A_(n-3) = ... $

$ U = D_(n-1) D_(n-2) dots D_1 A $

$ D := D_(n-1) D_(n-2) dots D_1 $

$ U = D A $

$ L := D^((-1)) $

$ A = L U $

Как выглядит эта матрица:

$ L = mat(
  1, ; 
  d_21, 1;
  d_31, d_32, 1;
  dots.v, dots.v, dots.v, dots.down;
  d_(n 1), d_(n 2), d_(n 3), dots, 1
) $

$ L U x = f $

$ U x := y $

$ L y = f $

Докажем необходимость в критерии работы метода Гаусса.

$ A = L U $

$ A^((k)) = L^((k)) U^((k)) $

Подматрица $k$-го порядка.

$ det A^((k)) = det L^((k)) det U^((k)) $

$ det L^((k)) = 1 $

$ det U^((k)) = u_11 u_12 dots u_(k k) $

$ u_(i i) != 0, quad i=overline(1 comma k) $

$ det A^((k)) != 0 $

`08 Октября 2024`
== Метод квадратного корня

Считаем, что все главные миноры не равны нулю.

$ A = L U $
$ U = D V $

$ D = mat(
  d_1;
  ,d_2;
  ,,dots.down;
  ,,,d_n;
) wide  
 V = mat(
  1, v_12, v_13, dots, v_(1 n);
  ,1,v_23,dots, v_(2 n);
  ,,,,dots.v;
  ,,,,1
) $


$ U = mat(
  d_1, d_1 v_12, d_1 v_13, dots, d_1 v_(1 n);
  ,    d_2,      d_2 v_23, dots, d_2 v_(2 n);
  ,            ,         , ,dots.v;
  ,,,,d_n
) $

$ d_i =u_(i i) quad v_(i j) = u_(i j)/(u_(i i)) $

Есть далее понятие LDV-разложение, оно определяется единственным образом.

Пусть $A$ - симметрическая:


$ A = A^* $

$ A = L D V $

$ A^* = (L D V)^* = V^* D^* L^* = (V^* D L^*) = A $

$ => L = V^*, V = L^* $

$ A = V^* D V $

Пусть матрица $A$ положительно определённая. Тогда по критерию Сильвестра все её угловые миноры положительны.

$ (A x, x) > 0 $

$ (D x, x) = (D V y, V y) = (A y, y) > 0 $

$ x = V y $


Все элементы на диагонали в $D$ больше нуля.

$ D^(1/2) = "diag"(sqrt(d_1), ..., sqrt(d_n)) $

$ A = V^* D^(1/2) D^(1/2) V $

$ W = D^(1/2) V $

$ W^* = V^* D^(1/2) $

$ A = W^* W $

$ A x = f $

$ W^* W x = f $

$ y:= W x $

$ W^* y = f $

$ a_(i j) = sum_(k = 1)^n W^*_(i k) W_(k j) = sum_(k = 1)^n W_(k i) W_(k j) = sum_(k = 1)^i W_(k i) W_(k j) $

$ a_11 = W_11 W_11 => W_11 = sqrt(a_11) $

$ i = 1 quad j = 2, 3, ..., n:\ 
  a_(1 j) = W_(1 1)W_(1 j) \
  W_(1 j) = a_(1 j) / W_11
$

$ a_(i i) = sum_(k=1)^i W_(k i) W_(k i) = sum_(k=1)^(i-1) W^2_(k i) + W^2_(i i) $

$ W_(i i) = sqrt(a_(i i) - sum_(k=1)^(i-1) W_(k i)^2) $

$ a_(i j) = sum_(k=1)^(i-1) W_(k i) W_(k j) + W_(i i) W_(i j) $

$ W_(i j) = (a_(i j) - sum_(k=1)^(i-1) W_(k i) W_(k j))/ W_(i i) $

$ j = i+1, ..., n $

Это был метод квадратного корня.

Если матрица симметрическая, но не является положительно определённой, значит решением будут комплексные числа, ничего страшного).

== QR разложение (Пуэр разложение)

#def[
Гиперплоскость с нормалью $p = vec(p_1, dots.v, p_n)$, проходящая через начало координат
$ (p, x) = p_1 x_1 + p_2 x_2 + dots + p_n x_n = 0 $
]

Ортогональное преобразование отражения относительно гиперплоскости, проходящей через начало координат:

$ y = x - 2 ((p, x))/((p, p)) p $


$ y = x - 2 ((p, x))/((p, p)) p = x - 2/((p, p)) (p,x) p = x - 2/((p, p)) p p^* x =\
  = (I - 2/((p,p)) p p^* )x 
$

$ P := I - 2/((p,p)) p p^* $

$ y = P x $

*Свойства:*
+ $P^2 = I$

  $ P^2 x = P(P x) = P( x - 2((p,x))/((p,p)) p) = P x - 2((p,x))/((p,p))P p = \
    = P x - 2 ((p/x))/((p,p)) (p - 2 ((p,p))/((p,p)) p)
    = x - cancel(2 ((p, x))/((p, p))) + cancel(2 ((p,x))/(p, p)) = x  $

+ $P^* =P$
    
  $ P^* = ( I - 2/((p,p)) p p^*)^* = I - 2/((p, p)) p^(**) p^* = p $

+ Матрица ортогональная $P^* P = P^2 = I$

+ Преобразование отражения не изменится, если вместо вектора нормали $p$ взять вектор $beta p$

  $ P = I - 2/((beta p, beta p)) (beta p) (beta p)^* = I - 2/((p, p)) p ^* $

+ Вектор нормали $p$ можно определить как разность между исходным и отражённым вектором: $p = x - y$

  $ x - y = 2 ((p,x))/((p,p)) p $

+ $p_1 = p_2 = dots = p_k = 0$

  $=> y_1 = x_1, y_2 = x_2, dots, y_k = x_k$

  Докажем:
  
  $ y_i = x_i - 2p_i / (sum_(l=1)^n p_l^2) sum_(l=1)^n p_l x_l $

  Из этой формулы очевидно (я хз??? он так сказал на лекции)

+ $p_1 = p_2 = dots = p_k = 0 " и " x_(k+1) = dots = x_n = 0 $

  $=> y_(k+1) = dots = y_n = 0$

$ A = Q R,\ "где" Q "ортогональная матрица, а" R "верхняя треугольная" $

Получим матрицу $A_1$ следующим образом:

$ A_1 = P_1 A $

Как должен выглядеть первый столбец:

$ vec(a_(11)^((1)), 0, 0, dots.v, 0) $

Введём обозначение столбцов матрицы $A$ и $A_1$:
$ a_(j), a^((1))_j$.

$ 
  vec(a_11, a_21, dots.v, a_(n 1)) = a_1
  quad
  vec(a_(11)^((1)), 0, dots.v, 0) = a_1^((1)) 
$

$ P_1 a_1 = vec( a_11^((1)), 0, 0, dots.v, 0) $

Вектор нормали $p^((1)) = a_1 - a_1^((1))$ (разность между исходным и отражённым)

$ &p^((1))_1 = a_11 - a_11^((1)) \
  &p^((1))_2 = a_21 \
  &dots.v \
  &p^((1))_n = a_(n 1)
$

При ортогональном преобразовании длина вектора не меняется:

$ norm(a^((1))_1)_2^2 = norm(a_1)_2^2 $

$ (a_11^((1)))^2  = sum_(l=1)^n a_11^2 $

$ a_11^((1)) = plus.minus sqrt(sum_(l=1)^n a_(l i)^2) $

$ p_1^((1)) = a_11 plus.minus sqrt(sum_(l=1)^n a_(l i)^2) $

$ sigma_1 := cases(
  1 when a_11 >= 0,
  -1 when a_11 < 0
) $

$ p_1^((1)) = a_11 + sigma_1 sqrt(sum_(l=1)^n a_(l i)^2) $

$ a_j^((1)) = a_j - 2 ((p^((1)), a_j)) / ((p_1^((1)), p_1^((1)))) p^((1)) where j = 2, ..., n $

По этой формуле полностью будет определена матрица $A_1$.

`15 Октября 2024`

$ mat(
  a_11^((1)), a_12^((1)), dots, a_(1,k-1)^((1)), a_(1 k)^((1)), a_(1,k+1)^((1)), dots, a_(1 n)^((1));
  ,a_22^((2)), dots, a_(2,k-1)^((2)), a_(2 k)^((1)), a_(2, k+1)^((1)), dots, a_(2 n)^((1));
  ,,dots.down;
  ,,,a_(k-1,k-1)^((k-1)), a_(k-1,k)^((k-1)), a_(k-1, k+1)^((k-1)), dots, a_(k-1, n)^((k-1));
  ,,,,a_(k k)^((k)), a_(k,k+1)^((k)), dots, d_(k,n)^((k));
  ,,,,a_(k+1, k)^((k+1)), a_(k+1,k+1)^((k+1)), dots, d_(k+1,n)^((k+1));
  ,,,,a_(n-1, k)^((n-1)), a_(n-1,k+1)^((n-1)), dots, d_(n-1,n)^((n-1));
) $

$ A_k = P_k A_(k - 1) $

$ a_k^((k)) = P_k a_k^((k-1)) = mat(
  a_(1 k)^((1));
  a_(2 k)^((2));
  dots.v;
  a_(n-1,k)^((k-1));
  a_(k k)^((k));
  0;
  dots.v;
  0
) $

$ p^((k)) = a_k^((k-1)) - a_k^((k)) $

$ p_l^((k)) = 0, quad l = 1, 2, dots, n-1 $

$ p_k^((k)) = a_(k k)^((k-1)) - a_(k k)^((k)) $

$ p_l^((k)) = a_(l k)^((k-1)), quad l = k+1, dots, n $

$ norm(a_k^((k+1)))_2^2 = norm(a_k^((k)))_2^2 $

$ a_(k k)^((k)) = plus.minus sqrt(sum_(l=k)^n (a_(l k)^((k-1)))^2)  $

$ sigma_k = cases(
  1 comma space a_(k k)^((k-1)) >= 0,
  -1 comma space a_(k k)^((k-1)) < 0
) $

$ p_k^((k)) = a_(k k)^((k-1)) + sigma_k sqrt(sum_(l=k)^n (a_(l k)^((k-1)))^2) $

$ a_j^((k)) = a_j^((k-1)) - 2 (p^((k)), a_j^((k-1))) / ((p^((k)), p^((k)))) p^k $

$ a_(i j)^((k)) = a_(i j)^((k-1)) - 2 p_i^((k)) / (sum_(l=n)^n (p_l^((k)))^2) sum_(l=k)^n p_l^((k)) a_(l j)^((k-1))  $

$ j = k+1, dots, n $

$ p_1^((k)) = p_2^((k)) = dots = p_(k-1)^((k)) = 0 $

Нужно сделать $n-1$ шаг. Рассмотрим весь процесс в целом.

$ R = A_(n-1) = P_(n-1) A_(n-2) = P_(n-1) P_(n-2) A_(n-3) = dots =\
  = P_(n-1) P_(n-2) dots P_2 P_1 A
$

$ Q := P_1 P_2 dots P_(n-1) $

Матрица $Q$ ортогональная (как произведение ортогональных матриц).

$ Q^* = P_(n-1)^* dots P_1^* = P_(n-1) dots P_1 $

$ R = Q^* A $

$ A = Q R $

Требования для существования этого разложения:

Рассмотрим случай, когда мы делим на ноль.

$ sum_(l=k)^n (p_l^((k)))^2 = 0  <=> p_l^((k)) = 0, space l = k, dots, n $

Подберём вектор нормали, который будет ненулевым, но оставит вектор на месте. Чтобы формально в алгоритме шаг выполнился. 

Например, $ p^((k)) = vec(0, dots.v, 0, sqrt(2), 0, dots.v, 0) $

Теперь можем говорить, что QR разложение существует для любых матриц.

// $ a_(l k)^((k-1)) = 0, quad l = k, ..., n $

Рассмотрим вопрос решения уравнений с помощью QR разложений.

$ A x = f $
$ Q R x = f quad $
$ R x = Q^* f $

Распишем $Q^*$:

$ Q^* = P_(n-1) dots P_1 $

$ Q^* f := g $

$ f^((1)) = P_1 f $
$ f^((2)) = P_2 f^((1)) $
$ g= f^((n-1)) = P_(n-1) f^((n-2)) $

$ f^((1)) = f - 2 ((p^((1)), f)) / ((p^((1)), p^((1)))) p^((1)) $

$ f^((k)) = f^((k-1)) - 2 p^((k)) / (sum_(l=k)^(n) (p_l^((k)))^2 ) sum_(l = k)^n p_l^((k)) f_l^((k-1)), quad i = k, k+1, dots, n $

== Метод окаймления
Способ отыскания обратной матрицы.

$ A = A_n = mat(augment: #(hline:-1, vline:-1),
  a_11, a_12, dots, a_(1 n-1), a_(1 n);
  a_21, a_22, dots, a_(2 n-1), a_(2 n);
  dots.v;
  a_(n-1, 1), a_(n-1, 2), dots, a_(n-1, n-1), a_(n-1, n);
  a_(n,1), a_(n,2), dots, a_(n, n-1), a_(n, n)
) $

$ A_(n - 1) $

$ v_n = (a_(n 1), a_(n 2), dots, a_(n,n-1)) $

$ u_n = (a_(1 n), a_(2,n), dots, a_(n-1, n))^* $

$ A_n = mat(A_n-1, u_n; v_n, a_(n n)) $

$ A_n^(-1) = mat(P_(n-1), r_n; q_n, 1/alpha_n) = D_n $

Предполагаем, что $A_(n-1)^(-1)$ известна 

$ mat(A_n-1, u_n; v_n, a_(n n)) mat(P_(n-1), r_n; q_n, 1/alpha_n) = mat(E, 0; 0, 1) $

$ cases(
 A_(n-1) P_(n-1) + u_n q_n = E,

 v_n P_(n-1) + a_(n n) q_n = 0,

 A_(n-1) r_n + u_n/alpha_n = 0,

 v_n r_n + a_(n n)/alpha_n = 1 
) $

Возмьём третье и выразим:
#cbox[
$ r_n =  -A_(n-1)^(-1) u_n / alpha_n $]

Подставляем в четвёртое:

$ - (v_n A_(n-1)^(-1) u_n) / alpha_n + a_(n n) / alpha_n  = 1 $

#cbox[
$ alpha_n = a_(n n) - v_n A_(n-1)^(-1) u_n $]

// Перейдём ко второму:

#cbox[
$ P_(n-1) = A_(n-1)^(-1) - A_(n-1)^(-1) u_n q_n $]

$ v_n (A_(n-1)^(-1) - A_(n-1)^(-1) u_n q_n) + a_(n n)q_n = 0 $

$ v_n A_(n-1)^(-1) + (a_(n n)- v_n A_(n-1)^(-1) u_n) q_n = 0 $

$ v_n A_(n-1)^(-1) + alpha_n q_n = 0 $

#cbox[
$ q_n = - (v_n A_(n-1)^(-1))/alpha_n $]

+ $ -A_(n-1)^(-1) u_n wide (beta_(1 n), dots, beta_(n-1, n))^* $
+ $ -v_n A_(n-1)^(-1) wide (gamma_(n 1), dots, gamma_(n, n-1)) $
+ $ alpha_n = a_(n n) + sum_(i=1)^(n-1) a_(n i) beta_(i n) = a_(n n) + sum_(i=1)^(n-1) a_(i n) gamma_(n i) $
+ $ d_(i k) = d_(i k) + (beta_(i k) gamma_(n k))/ alpha_n, quad i, k <= n-1 $
  $ d_(i n) = beta_(i n) / alpha_n quad d_(n k) = gamma_(n k) / alpha_n $

Работа в первом блоке завершена. Можно сдавать первый коллоквиум по прямым методам. Д947.

= Итерационные методы

$ A x = f $

$ { x^((k)) } $

$ lim_(k->oo) x^((k)) = x^* $

$ lim_(k->oo) norm(x^((k)) - x^*) = 0 $

Каноническая форма записи:

$ B (x^((k-1)) - x^((k)))/(tau_k) + A x^((k)) = f, space k = 0, 1, dots $

$x^((0)) "начальное приближение"$

$B - "невырожденная матрица"\ tau - "итерационный параметр (может меняться на разных шагах)" $
Если $tau$ не зависит от $k$, то метод называется стационарным.

Если методы сходятся, то они сходятся к точному решению. Поэтому нужно доказывать только сходимость.

Скорость сходимости. Надо добиться наивысшей скорости сходимости. Зависит от параметра и от матрицы $B$.

Сейчас нас будут интересовать только стационарные методы:

$ B (x^((k-1)) - x^((k)))/(tau) + A x^((k)) = f, space k = 0, 1, dots $

$ B x^((k+1)) = B x^((k)) - tau A x^((k)) + tau f $

$ B x^((k+1)) = (B - tau A)x^((k)) + tau f $

Матрица $B$ должна быть простая, чтобы систему можно было просто решать, а то игра не стоит свеч.

Показатели:

+ Вектор погрешности или вектор ошибки $z^((k)) = x^((k)) - x^*$.
+ Вектор невязки: $r^((k)) = A x^((k)) - f$

$ lim_(k->oo) norm(z^((k))) = 0 - "сходимость к точному решению" $

Добавим и отнимем точное решение:

$ B (x^((k-1)) -x^* - x^((k)) + x^*)/tau  + A x^((k)) - A x^* = 0 $

$ B (z^((k+1)) - z^(k)) / tau + A z^((k)) = 0 $

$ B z^((k+1)) = B z^((k)) + tau A z^((k)) $

$ z^((k+1)) = (I - tau B^(-1)A) z^((k)) $

Матрица перехода $S = (I - tau B^(-1)A)$:

$ z^((k+1)) = S z^((k)) $

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

$ norm(z^((k))) <= 1/e norm(z^((0))) $

Сколько шагов нужно сделать, чтобы начальная ошибка уменьшилась в $e$ раз.

$ z^((k+1)) = S z^((k)) => z^((k)) = S^k z^((0)) $

$ norm(z^((k))) <= norm(S)^k norm(z^((0))) $

Самый грубый подход. Потребуем, чтобы:

$ norm(S)^k <= 1/e $

Вытащим отсюда $k$:

$ k ln norm(S) <= -1 $

Для сходящихся методов норма матрицы всегда меньше единицы. Мы это докажем.

$ k >= 1 / (-ln norm(S)) $

$ R = - ln norm(N) - "ассимптотическая скорость сходимости" $

Теперь рассмотрим основополагающую теорему о сходимости стационарных методов. 

#bbox[Критерий сходимости][
 Для того, чтобы стационарный метод сходился при любом начальном приближении, небходимо и достаточно, чтобы все собственные значения матрицы перехода были по модулю меньше единицы.

 $ abs((lambda(S))) < 1 $
]

Критерий плох тем, что надо постоянно искать собственные значения. Поэтому это теорема базовая.

Докажем необходимость. Пусть метод сходится при любом начальном приближении. Докажем, что все собственные значения по модулю меньше единицы. От противного. Пусть метод сходится, но существует хотя бы одно $abs(lambda(S)) >= 1$. Обозначим за $u$ собственный вектор, который соответствует этому собственному значению. $S u = lambda u$. В качестве начального приближения возьмём $x^((0)) = x^* + u$. Тогда $z^((0)) = u$. 

$ z^((k)) = S^k z^((0)) = S^k u = lambda^k u $

$ norm(z^((k))) = norm(lambda^k u) = abs(lambda)^k norm(u) $

$ u!=0, quad k->oo $

$ norm(z^((k))) arrow.not.r_(k->oo) 0 $

Получили противоречие.

Сформулируем лемму, которая понадобится для доказательства достаточности.

#bbbox[Лемма][
  Пусть все собственные значения матрицы $S$:

  $ abs(lambda_i (s)) < q, space i = 1, dots, n $

  Тогда существует такая невырожденная матрица $T$, что матрица $Lambda = T S T^(-1)$
  удовлетворяет такому условию: $norm(Lambda)_oo <= q$
]

Напоминание:
$ Lambda tilde S => "собственные значения у них совпадают" $

Теперь докажем достаточность.

$ abs(lambda (S)) < 1 $

$ q < 1 $

$ max_i abs(lambda_i (s)) < q $

Надо показать, что $lim_(k->oo) norm(z^((k))) = 0$.

Теперь привлекаем лемму. Существует такая невырожденная матрица $T$, что $norm(Lambda)_oo <= q $.

$ S = T^(-1) Lambda T $

Подставим $S$ в определение $z^((k))$.

$ S^k = T^(-1) Lambda^k T $

$ z^((k)) = S^k z^((0)) = T^(-1) Lambda^k T z^((0)) $

$ norm(z^((k)))_oo = norm(S^k z^((0)) = T^(-1) Lambda^k T z^((0)))_oo <= \
  <= norm(T^(-1))_oo norm(Lambda)_oo^k  norm(T)_oo norm(z^((0)))_oo <= \
  <= mu_oo(T) norm(Lambda^k)_oo norm(z^((0)))_oo <= mu_oo (T) q^k norm(z^((0)))_oo -->_(k->oo) 0
  $


#bbox[Теорема][
Для сходимости двуслойного стационарного метода при любом начальном приближении достаточно, чтобы хотя бы одна из норм матрицы перехода $S$, согласованная с какой-нибудь векторной, была меньше единицы.

$ norm(S) < 1 $
]

Собственное значение и вектор:

$ S u = lambda u $

$ abs(lambda) norm(u) = norm(lambda u) = norm(S u) <= norm(S) norm(u) $

$ abs(lambda) <= norm(S) $

Собственное значение всегда ограничено нормой.

$ norm(S) < 1 => abs(lambda) < 1 $
