Материал из eSyr's wiki.

Предыдущая лекция | Следующая лекция

Глава 1. Численные методы линейной алгебры

Параграф 12. Предварительное преобразование к верхней почти треугольной форме

A_k = Q_kR_k, Q_k⁻¹ = Q_k^T, R — ВТФ, A_{k + 1} = R_kQ_k, k = 0 ,1, … (1)

Оценка количества действий:

• A — произвольная: O(m³)
• A — ВПТФ: O(m²)
• A — трёхдиагональная: O(m)

Лемма 1. Пусть C = B × A, B — ВТФ, A — ВПТФ. Тогда C — ВПТФ.
Доказательство. Распишем c_ij = ∑_{α = 1}^mb_iαa_αj = (*). Так как B — ВТФ, то b_iα = 0 при i > α. Тогда (*) = ∑_{α = i}^mb_iαa_αj = (**). Так как A — ВПТФ, то a_αj = 0 при j + 1 < α. Тогда (**) = ∑_{α = i}^{j + 1}b_iαa_αj. Отсюда следует, что c_ij = 0, i > j + 1 ⇒ C — ВПТФ.

Задача. Пусть C = B × A, B — ВПТФ, A — ВТФ. Доказать, что C — ВПТФ.

Теперь мы можем сказать, что QR-алгоритм инвариантен относительно ВПТФ, он её не портит. Это не значит, что он её сохраняет, так как в результате мы должны получить ВТФ. Докажем это:

Пусть A = A₀ — ВПТФ. Тогда мы можем выразить Q_k = R_k⁻¹A_k, Q₀ = R₀⁻¹A₀. Но так как A₀ — ВПТФ, то Q₀ тоже ВПТФ. Рассмотрим A_{k + 1} = R_kQ_k. Но R_k — ВТФ, Q_k — ВПТФ. Тогда A_{k + 1} тоже ВПТФ. В результате A_k, k = 0, 1, … — ВПТФ.

На этом завершим первую главу.

Глава 2. Интерполирование и приближение функций

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

Сейчас мы вспомним и полином Лагранжа, затем рассмотрим интерполяционный полином Ньютона. Далее начнём строить полином Эрмита и покажем, для чего нужны они. Ну и наилучшее квадратичное приближение функции.

Это наиболее древняя область, но она актуальна и поныне, поскольку без интерполирования математик обойтись не может. Это актуально для специальных функций, которые не выражаются через элементарные, но используются для решения дифуров. Обычно используются таблицы, но возникает вопрос, что делать, когда нужно получить значение для аргумента, которого нет в таблице. Первое решение — сгущать таблицы. Ещё одно — пользоваться интерполированием.

Где ещё требуется интерполирование — есть домна, за ней надо следить. Надо наставить датчиков, но много их поставить нельзя, иначе она разрушится, а информацию надо получать по всему профилю. Поэтому надо использовать интерполяцию.

Ешё одно применение — разностные схемы.

Ещё одно применение — экстраполирование.

Обычно будем использовать приближение полиномом.

Параграф 1. Введение

Суть задачи: есть отрезок x ∈ [a, b]. Мы разбиваем его несовпадающими узлами a ≤ x₀ < x₁ < x₂ < … < x_n ≤ b — {x_i}₀ⁿ — узлы интерполирование (n + 1 штук). Есть f(x_i) = f_i, i = 0…n. Если полином P_n(x) = a₀ + a₁x + … a_nxⁿ такой, что P_n(x_i) = f(x_i), i = 0…n (1), то данный полином называется интерполяционным для f (интерполирует f).

Докажем, что интерполирующий полином существует и единственный. Рассмотрим систему:

• a₀ + a₁x₀ + … a_nx₀ⁿ = f₀
• a₀ + a₁x₁ + … a_nx₁ⁿ = f₁
• ...
• a₀ + a₁x_n + … a_nx_nⁿ = f_n

Определитель этой системы D_{n + 1} = ∏_{n ≥ i > j ≥ 1}(x_i − x_j) ≠ 0 (x_i ≠ x_j) — определитель Вандермонда.

Параграф 2. Интерполяционная формула Лагранжа. Погрешность формулы Лагранжа

Есть {x_i}₀ⁿ — узлы интерполирования. Заданы f(x_i) = f_i, i = 0…n. Задача построить полином L_n(x): L_n(x_i) = f(x_i), i = 0…n (1). Будем искать полином в виде L_n(x) = ∑_{k = 0}ⁿc_k(x)f(x_k) (2), c_k(x) — полином n-й степени. Введём многочлен ω(x) = (x − x₀)(x − x₁)…(x − x_n) = ∏_{i = 0}ⁿ(x_i − x_j). Найдём производную:

• ω(x) = [ ](x − x_k)
• ω'(x) = [ ] + [ ]'(x − x_k)
• ω'(x_k) = ∏_{k ≠ j}ⁿ(x_k − x_j)

Возьмём c_k(x) = (ω(x))/((x − x_k)ω'(x_k))

Тогда L_n(x) = ∑_{k = 0}ⁿ(ω(x))/((x − x_k)ω'(x_k)) × f(x_k) (3)

Пример

Для трёх точек x₀, x₁, x₂ L₂(x) = ^{(x − x₀)(x − x₁)(x − x₂)}/_{(x − x0)(x0 − x1)(x0 − x2)} + ^{(x − x₀)(x − x₁)(x − x₂)}/_{(x − x1)(x1 − x0)(x1 − x2)} + ^{(x − x₀)(x − x₁)(x − x₂)}/_{(x − x2)(x2 − x0)(x2 − x1)} = ^{(x − x₁)(x − x₂)}/_{(x0 − x1)(x0 − x2)} + ^{(x − x₀)(x − x₂)}/_{(x1 − x0)(x1 − x2)} + ^{(x − x₀)(x − x₁)}/_{(x2 − x0)(x2 − x1)}.

Погрешность полинома

Погрешность интерполяционной формулы Лагранжа ψ_Ln(x) = f(x) − L_n(x). Во всех узлах она будет совпадать: ψ_Ln(x_i) = 0.

Функция настолько гладкая, насколько требуется. В данном случае функция должна иметь n + 1-ю производную.

В данном случае рассматриваем g(s) = f(s) − L_n(s) − Kω(s), ω(s) = ∏_{i = 0}ⁿ(s − x_i)

f(s) − L_n(s) = Kω(s)

В данном случае получали оценку f(s) − L_n(s) = (f^{(n + 1)}(ξ)/(n + 1)!) × ω(x)

|f(x) − L_n(x)| ≤ (M_{n + 1}/(n + 1)!) × |ω(x)|

Замечание 1. Если f(x) — полином степени ≤ n, то полином Лагранжа даёт точное приближение, то есть погрешность тождественно равна 0.

Замечание 2. По большому числу узлов интерполяцию проводят крайне редко. Можно подобрать узлы так, что при увеличении количества узлов сходимости не будет.

Параграф 3. Разделённые разности

Разделённые разности — подводка к интерполяционным полиномам Ньютона.

Есть отрезок x ∈ [a, b]. Мы разбиваем его несовпадающими узлами a ≤ x₀ < x₁ < x₂ < … < x_N ≤ b — {x_i}₀^N — узлы интерполирования (n + 1 штук). Есть интерполируемая функция f(x), f(x_i) = f_i, i = 0…N.

Разделённой разностью первого порядка, построенной по узлам x_i, x_j называется отношение f(x_i, x_j) = (f(x_j) − f(x_i)/(x_j − x_i))

Рассмотрим разделённые разности между соседними узлами:

• f(x₀, x₁) = (f(x₁) − f(x₀)/(x₁ − x₀))
• f(x₁, x₂) = (f(x₂) − f(x₁)/(x₂ − x₁))

Разделённая разность второго порядка:

• f(x_i, x_j, x_k) = (f(x_j, x_k) − f(x_i, x_j)/(x_k − x_i))
• f(x₀, x₁, x₂) = (f(x₁, x₂) − f(x₀, x₁)/(x₂ − x₀))
• f(x₁, x₂, x₃) = (f(x₂, x₃) − f(x₁, x₂)/(x₃ − x₁))

Разделённая разность k+1-го порядка: пусть есть f(x_j, x_{j + 1}, …, x_{j + k}), f(x_{j + 1}, x_{j + 2}, …, x_{j + k + 1}) — разности k-го порядка. Тогда f(x_j, x_{j + 1}, …, x_{j + k + 1}) = (f(x_{j + 1}, x_{j + 2}, …, x_{j + k + 1}) − f(x_j, x_{j + 1}, …, x_{j + k}))/(x_{j + k + 1} − x_j) — разделённая разность k + 1-го порядка.

Докажем, что разделённые разности выражаются через значения функций в точках. Рассмотрим:

• ω(x) = (x − x₀)…(x − x_N) = ω_{0, N}(x)
• ω_{α, β}(x) = (x − x_α)…(x − x_β)
• ω_{α, β}'(x) = ∏_{σ, (α ≠ k)}^β(x_k− x_σ), k = 1

Утверждение. Разделённая разность k-го порядка, построенная по узлам x₀, … x_k может быть записана как f(x₀, x₁, …, x_k) = ∑_{i = 0}^k(f(x_i))/(ω_{0, k}'(x_i)) (1)

Доказательство (методом полной математической индукции).

1. f(x₀, x₁) = (f(x₁) − f(x₀))/(x₁ − x₀) = f(x₀)/(x₀ − x₁) + f(x₁)/(x₁ − x₀)
2. По предпололжению индукции:
   • f(x₀, x₁, …, x_l) = ∑_{i = 0}^l(f(x_i))/(ω_{0, l}'(x_i))
   • f(x₁, x₂, …, x_{l + 1}) = ∑_{i = 1}^{l + 1}(f(x_i))/(ω_{1, l + 1}'(x_i))
3. Рассмотрим f(x₀, x₁, …, x_{l + 1})
   • f(x₀, x₁, …, x_{l + 1}) = 1/(x_{l + 1} − x₀) × (∑_{i = 1}^{l + 1}(f(x_i))/(ω_{1, l + 1}'(x_i)) − ∑_{i = 0}^l(f(x_i))/(ω_{0, l}'(x_i))) = 1/(x_{l + 1} − x₀) × ((f(x_{l + 1}))/(ω_{1, l + 1}'(x_{l + 1})) − (f(x₀))/(ω_{0, l}'(x₀)) + ∑_{i = 0}^l((f(x_i))/(ω_{1, l + 1}'(x_i)) − (f(x_i))/(ω_{0, l}'(x_i)))) = (***)
     • (f(x_{l + 1}))/((x_{l + 1} − x₀)ω_{1, l + 1}'(x_{l + 1})) = (f(x_{l + 1}))/(ω_{0, l + 1}'(x_{l + 1}))
     • (f(x₀))/((x_{l + 1} − x₀)ω_{0, l}'(x₀)) = (f(x₀))/(ω_{0, l + 1}'(x₀))
     • (f(x_i)(x_i − x₀))/(ω_{1, l + 1}'(x_i)) − (f(x_i)(x_i − x₀))/(ω_{0, l}'(x_i)) = (f(x_i))/(ω_{0, l + 1}'(x_i))
   • (***) = ∑_{l = 0}^{l + 1}f(x_i/(ω_{0, l + 1}'(x_i)))