Чем отличается полином лагранжа от полинома ньютона

3. Интерполяция

Пусть функция задана набором точек на интервале :

, , (3.1)

Задача интерполяции – найти функцию , принимающую в точках те же значения . Тогда, условие интерполяции:

(3.2)

При этом предполагается, что среди значений нет одинаковых. Точки называют узлами интерполяции.

Если ищется только на отрезке – то это задача интерполяции, а если за пределами первоначального отрезка, то это задача экстраполяции.

Задача нахождения интерполяционной функции имеет много решений, так как через заданные точки можно провести бесконечно много кривых, каждая из которых будет графиком функции, для которой выполнены все условия интерполяции. Для практики важен случай интерполяции функции многочленами:

, (3.3)

При этом искомый полином называется интерполяционным полиномом.

При построении одного многочлена для всего рассматриваемого интервала для нахождения коэффициентов многочлена необходимо решить систему уравнений, построенную на основе полинома (3.3). Данная система содержит уравнение, следовательно, с ее помощью можно определить коэффициент. Поэтому максимальная степень интерполяционного многочлена , и многочлен принимает вид

, (3.4)

3.2. Локальная и глобальная интерполяция

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

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

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

3.3. Кусочно-линейная интерполяция

Простейшим и часто используемым видом локальной интерполяции является линейная (или кусочно-линейная) интерполяция. Она заключается в том, что узловые точки соединяются отрезками прямых (рис.3.1), то есть через каждые две точки и проводится прямая, то есть составляется полином первой степени:

, при (3.5)

Рис. 3.1. Кусочно-линейная интерполяция

Коэффициенты и разные на каждом интервале , и находятся из выполнения условий интерполяции на концах отрезка:

(3.6)

Из системы уравнений (3.6) можно найти коэффициенты:

, (3.7)

При использовании кусочно-линейной интерполяции сначала нужно определить интервал, в который попадает значение x, а затем подставить его в выражение (3.5), используя коэффициенты для данного интервала.

3.4. Кусочно-квадратичная интерполяция

В случае квадратичной интерполяции, для каждых трех узловых точек , , , строится уравнение параболы:

, при (3.8)

Рис.3.2. Кусочно-квадратичная интерполяция

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

(3.9)

Из системы уравнений (3.9) можно найти коэффициенты:

(3.10)

3.5. Многочлен Лагранжа

При глобальной интерполяции на всем интервале строится единый многочлен. Одной из форм записи интерполяционного многочлена для глобальной интерполяции является многочлен Лагранжа:

(3.11)

где – базисные многочлены степени n:

(3.12)

То есть многочлен Лагранжа можно записать в виде:

(3.13)

Многочлен удовлетворяет условию . Это условие означает, что многочлен равен нулю при каждом кроме , то есть – корни этого многочлена. Таким образом, степень многочлена равна n и при обращаются в ноль все слагаемые суммы, кроме слагаемого с номером , равного .

Выражение (3.11) применимо как для равноотстоящих, так и для не равноотстоящих узлов. Погрешность интерполяции методом Лагранжа зависит от свойств функции , от расположения узлов интерполяции и точки x. Полином Лагранжа имеет малую погрешность при небольших значениях n (n<20). При больших n погрешность начинает расти, что свидетельствует о том, что метод Лагранжа не сходится (то есть его погрешность не убывает с ростом n).

Многочлен Лагранжа в явном виде содержит значения функций в узлах интерполяции, поэтому он удобен, когда значения функций меняются, а узлы интерполяции неизменны. Число арифметических операции, необходимых для построения многочлена Лагранжа, пропорционально и является наименьшим для всех форм записи. К недостаткам этой формы записи можно отнести то, что с изменением числа узлов приходится все вычисления проводить заново.

Кусочно-линейная и кусочно-квадратичная локальные интерполяции являются частными случаями интерполяции многочленом Лагранжа.

3.6. Многочлен Ньютона

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

Разделенные разности нулевого порядка совпадают со значениями функции в узлах. Разделенные разности первого порядка определяются через разделенные разности нулевого порядка:

(3.14)

Разделенные разности второго порядка определяются через разделенные разности первого порядка:

(3.15)

Разделенные разности k-го порядка определяются через разделенные разности порядка :

(3.16)

Используя понятие разделенной разности интерполяционный многочлен Ньютона можно записать в следующем виде:

(3.17)

За точностью расчета можно следить по убыванию членов суммы (3.17). Если функция достаточно гладкая, то справедливо приближенное равенство . Это приближенное равенство можно использовать для практической оценки погрешности интерполяции: .

Интерполяционные полиномы Лагранжа и Ньютона

где значения х_п и у_п заданы. Такой полином определяется единственным образом и может быть представлен в следующем виде:

где сомножители 1_п(х) задаются как

Полином (8.2) называется интерполяционным полиномом Лагранжа степени N. Для вычисления значения интерполяционного полинома в некорой точке х, не прибегая к явному вычислению коэффициентов полинома, применяется схема Невилла (Neville). На первом этапе зададим значения

Далее значения вспомогательных полиномов в точке х определяются согласно рекуррентному соотношению

и окончательно мы получаем

Пример 8.1 (интерполяционный полином Лагранжа)

В табл. 8.1 представлены значения у_п в соответствующих точках х_п и требуется провести гладкую кривую через эти точки.

Определим сомножители 1_п(х), которые для данного набора данных записываются как

Тогда интерполяционный полином Лагранжа имеет вид

Предположим, что нам нужно вычислить значение Р$(х) в промежуточной точке х = 3, тогда схема Невилла (8.3) представляется в виде следующей таблицы:

Pf >(.v), я = 0. 3 — к

Рис. 8.1. Гладкая кривая, описывающая данные из табл. 8.1:

О — данные;—интерполяционный полином Лагранжа

Представление интерполяционного полинома в форме (8.2) имеет простой вид, и поэтому оно часто используется для теоретического анализа. Однако с точки зрения практических вычислений такое представление применимо только для небольших значений N. С увеличением N сомножители 1_п(х) возрастают и начинают сильно осциллировать. Существует другое представление интерполяционного полинома, которое лучше подходит для вычислений, и оно основано на разделенных разностях. Пусть заданы различные точки х₀. x_N и значения в этих точках у₀, . у_jV. Разделенные разности нулевого порядка есть просто значения у_п:

Разделенные разности порядка k в точке х_п вычисляются но следующим рекурентным формулам:

Тогда, используя введенные обозначения, интерполяционный полином, удовлетворяющий условию (8.1), записывается в виде

Полином (8.4) называется интерполяционным полиномом Ньютона. Его можно представить и в другой форме:

Используя это представление, значение интерполяционного полинома Ньютона в некоторой точке х может быть вычислено простым и эффективным способом. Этот алгоритм называется схемой Горнера (Horner) и задается следующими рекурентными формулами:

Тогда требуемое значение полинома равно последнему элементу этой последовательности, т.е. Р^(х) = я₀.

Пример 8.2 (интерполяционный полином Ньютона) Рассмотрим данные, представленные в табл. 8.1. Разделенные разности, вычисленные на основе этих данных, показаны в табл. 8.2.

D^Xx), п = 0. 3 — к

VMath

Для понимания материалов настоящего раздела крайне желательно ознакомиться с разделом ПОЛИНОМ ОДНОЙ ПЕРЕМЕННОЙ.

Интерполяция

Интерполяция или интерполирование — приближенное или точное нахождение какой-либо величины по известным отдельным значениям этой же величины, или других величин, с ней связанных.

Различают интерполяцию и экстраполяцию: образно говоря, интерполяция занимается восстановлением неизвестной функции в промежутках между точками, где известны ее значения (например, при известных значениях $ f(1) $ и $ f(2) $ требуется оценить $ f(1.5) $), в то время как экстраполяция предполагает, что мы пытаемся «растянуть» множество задания функции (например, при известных $ f(1) $ и $ f(2) $ оцениваем $ f(2.5) $ или $ f(0.5) $).

Пример. Вставить пропущенные буквы: ИНТ. РПО. ЯЦИЯ ; . КСТРАПОЛЯЦ. .

Происхождение слова «интерполяция» ☞ ЗДЕСЬ.

Задача интерполяции решается в разных классах функций — полиномов алгебраических или тригонометрических, комбинаций экспонент; в классе рациональных функций. Начинаем изложение материала с самого простого случая —

Полиномиальная интерполяция

Задача. Построить полином $ y_<>=f(x) $, принимающий значения согласно следующей таблице: $$ \begin x & x_1 & x_2 & \dots & x_n \\ \hline y & y_1 & y_2 &\dots & y_n \end $$ Здесь числа $ \< x_, y_ \>_^n $ — рациональные (т.е. из $ \mathbb Q_<> $), вещественные (из $ \mathbb R_<> $) или комплексные (из $ \mathbb C_<> $); числа $ \ < x_\>_^n $ называются узлами интерполяции; искомый полином называется интерполяционным. Для однообразия любое из упомянутых множеств будем обозначать через $ \mathbb A_<> $.

Геометрическая интерпретация для случая $ \mathbb A_<> = \mathbb R_<> $: построить алгебраическую кривую $ y_<> = f(x) $, проходящую через заданные точки $ (x_,y_) $ плоскости $ (x_<>,y_<>) $.

Теорема. При $ x_,\dots,x_ $ различных существует единственный полином $ f(x) \in \mathbb A_<>[x] $ степени $ \le n_<>-1 $ такой, что

Доказательство. Коэффициенты полинома $ f(x)=A_+A_1x+\dots+A_x^ $ можно определить из системы уравнений $$ \left\ A_0+A_1x_1+\dots+A_x_1^&=&y_1 \\ A_0+A_1x_2+\dots+A_x_2^&=&y_2 \\ \dots & & \dots \\ A_0+A_1x_n+\dots+A_x_n^&=&y_n, \end \right. $$ которая имеет единственное решение поскольку определитель системы (см. ☞ определитель Вандермонда, формулы Крамера ) отличен от нуля. ♦

Иногда интерполяционным полиномом называют полином из теоремы — т.е. полином минимально возможной степени. Но я в дальнейшем не ставлю такого ограничения.

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

$$ f(x)+(x-x_1)\times \dots \times(x-x_n)q(x) \ , $$ где $ f_<>(x) $ – полином из предыдущей теоремы, а $ q_<>(x) \in \mathbb A[x] $ – произвольный полином.

Придумать упрощение схемы вычисления интерполяционного полинома для таблицы с набором узлов симметричным относительно $ 0_<> $:

$$ \begin x & -x_n & -x_ & \dots & -x_1 & x_1 & \dots & x_ & x_n \\ \hline y & y_ & y_ &\dots & y_ & y_1 & \dots & y_ & y_n \end $$

Решение ☞ ЗДЕСЬ.

Пусть имеется интерполяционная таблица

$$ \begin x & x_1 & x_2 & \dots & x_n \\ \hline y & F(x_1) & F(x_2) &\dots & F(x_n) \end $$ построенная для полинома $ F_<>(x) $ степени $ \ge n $. Какое отношение имеет интерполяционный полином, построенный по этой таблице, к полиному $ F_<>(x) $?

Вычисление интерполяционного полинома посредством решения системы линейных уравнений из теоремы в общем случае затрудняется тем обстоятельством, что матрица Вандермонда является плохо обусловленной, т.е. небольшая погрешность в представлении ее элементов или значений $ y_1,\dots,y_n $ может приводить к существенному изменению решений системы уравнений.

Пример. Интерполяционный полином для таблицы

$$ \begin x & 1 & 2 & 3 & 4 & 5 \\ \hline y & 1 & -2 & 33 & 166 & 481 \end $$ имеет вид $$ x^4-6\,x^2 + 6 .$$ Если немного изменить значения $ y $: $$ \begin x & 1 & 2 & 3 & 4 & 5 \\ \hline y & 1.05 & -2.05 & 32.95 & 166.05 & 480.95 \end $$ то полином изменится не очень существенно: $$ 0.9875x^4+0.125x^3-6.3875x^2+0.375x+5.95 \, . $$ Но вот если осуществить следующие изменения того же порядка малости $$ \begin x & 1 & 2 & 3 & 4 & 5 \\ \hline y & 1.02 & -2.05 & 33.05 & 165.95 & 481 \end $$ то полином меняется гораздо существенней: $$ 1.03x^4-0.361(6)x^3-4.495x^2-2.50(3)x+7.35 \, . $$

Следовательно, применение вычислительных методов решения систем линейных уравнений — типа метода Гаусса — к системе с матрицей Вандермонда столкнется с необходимостью строгого контроля округлений. ♦

Практическое построение интерполяционного полинома производится альтернативными алгоритмами — посредством вспомогательных промежуточных представлений полинома в специальных, сравнительно просто вычисляемых, видах. Самыми распространенными являются формы Лагранжа и Ньютона.

Интерполяционый полином в форме Лагранжа

Уравнение для определения интерполяционного полинома можно записать в детерминантной форме: $$ 0\equiv \left| \begin 1 & x_1 & x_1^2 & \dots & x_1^ & -y_1 \\ 1 & x_2 & x_2^2 & \dots & x_2^ & -y_2 \\ \vdots & & & & & \vdots \\ 1 & x_n & x_n^2 & \dots & x_n^ & -y_n \\ 1 & x & x^2 & \dots & x^ & — f(x) \end \right|_ \ . $$ Представив последний столбец в виде суммы, получим отсюда: $$ f(x) \equiv \left| \begin 1 & x_1 & x_1^2 & \dots & x_1^ & -y_1 \\ 1 & x_2 & x_2^2 & \dots & x_2^ & -y_2 \\ \dots & & & & & \dots \\ 1 & x_n & x_n^2 & \dots & x_n^ & -y_n \\ 1 & x & x^2 & \dots & x^ & 0 \end \right| \Big/ \left| \begin 1 & x_1 & x_1^2 & \dots & x_1^ \\ 1 & x_2 & x_2^2 & \dots & x_2^ \\ \dots & & & & \dots \\ 1 & x_n & x_n^2 & \dots & x_n^ \end \right| \ . $$ Теперь разложим определитель из числителя по последнему столбцу и вспомним выражение определителя Вандермонда. Обозначим $$ W(x) = (x-x_1)\times \dots \times (x-x_n) \ , $$ $$ W_j(x) = \frac =(x-x_1)\times \dots \times (x-x_) (x-x_)\times \dots \times (x-x_n) . $$ Тогда интерполяционный полином в форме Лагранжа записывается в виде: $$ f(x) \equiv \sum_^n \frac= $$ $$ = y_1\frac+ $$ $$ +y_2\frac+ \dots+ $$ $$ +y_n\frac<(x-x_1)\times \dots \times (x-x_)><(x_n-x_1)\times \dots \times (x_n-x_)> . $$

Образно говоря, каждое слагаемое в форме Лагранжа «отвечает» исключительно только за свой узел интерполяции (т.е. обеспечивает в нем нужное значение полинома) — и при этом «не портит» остальные узлы (обращается в них в нуль). Полиномы $$ \_j(x):= W_j(x)/W_j(x_j) \>_^n $$ называются базисными интерполяционными полиномами. Эта система полиномов действительно образует базис линейного пространства полиномов степеней, не превосходящих $ n-1 $.

а) $ W_j(x_j)=W^<\prime>(x_j) $;

б) $ \displaystyle \sum_^n _j(x) \equiv 1 $ .

Пример. Построить интерполяционный полином по таблице

$$ \begin x & -2 & -1 & 1 & 2 \\ \hline y & -29 & -8 & -2 & 7 \end $$ и с его помощью интерполировать значение неизвестной функции при $ x_<>=0 $.

Решение. Имеем: $ W_<>(x)=(x+2)(x+1)(x-1)(x-2) $, и полином в форме Лагранжа: $$ f(x)=\frac(x+1)(x-1)(x-2)+ \frac(x+2)(x-1)(x-2)+ \frac(x+2)(x+1)(x-2) + $$ $$ + \frac (x+2)(x+1)(x-1) \ . $$ Подставляем сюда $ x_<>=0 $:

Ответ. $ f_<>(0)=-3 $.

Рекурсивное вычисление коэффициентов

Заметим, что представление интерполяционного полинома в форме Лагранжа еще не дает его канонической формы, т.е. явного выражения для коэффициентов $ A_0,A_1,\dots, A_ $. Для получения последних требуется разложить по степеням переменной базисные интерполяционные полиномы. $$ _j(x)\equiv w_+w_x+\dots+ w_x^ \, . $$ Кстати, осуществив такие разложения, мы получим и возможность обращения матрицы Вандермонда: $$ \left( \begin 1& x_1 &x_1^2 & \ldots& x_1^\\ 1& x_2& x_2^2 & \ldots& x_2^\\ \vdots& \vdots & \vdots & & \vdots \\ 1 & x_& x_^2& \ldots& x_n^ \end \right)^= \left( \begin w_& w_& \ldots& w_\\ w_& w_& \ldots& w_\\ \dots & \dots & & \dots \\ w_& w_& \ldots& w_\\ \end \right)^ \, . $$ Подробности ☞ ЗДЕСЬ; решение этой задачи востребовано не только для целей интерполяции.

В настоящем пункте мы произведем «доводку» метода Лагранжа до коэффициентов интерполяционного полинома

Теорема. Пусть числа $ x_,\dots,x_n $ все различны. Для полинома

$$ W(x)=(x-x_)\times \dots \times (x-x_n) $$ справедливы следующие равенства Эйлера-Лагранжа : $$ \sum_^n \frac=\left\< \begin 0 & npu \ k < n-1; \\ 1 & npu \ k= n-1. \end\right. $$

Доказательство. Построим интерполяционный полином по следующей таблице: $$ \begin x & x_1 & \dots & x_n \\ \hline y & x_1^k & \dots & x_n^k \end $$ С одной стороны, ответ известен заранее: $ f(x)\equiv x^k $. С другой стороны, формула Лагранжа дает его же в виде суммы: $$ x^k \equiv \sum_^n x_j^k \frac \equiv x^\sum_^n \frac+ \dots \ , $$ где в правой части стоит разложение по убывающим степеням $ x_<> $. В этом тождестве степени полиномов слева и справа должны быть одинаковыми. Если $ k< n-1 $, то старший коэффициент правого полинома должен обратиться в нуль. Если же $ k= n-1 $, то должны совпасть старшие коэффициенты обоих полиномов. ♦

Для любого полинома $ g(x)\in \mathbb C[x], \deg g < n-1 $ выполняется равенство

В следующем результате изменен порядок нумерации коэффициентов интерполяционного полинома!

Теорема. Обозначим

$$ \sigma_k = \sum_^ \frac>(x_j)> \ , \ \tau_k = \sum_^ y_j \frac>(x_j)> \ . $$ Имеют место равенства, связывающие коэффициенты интерполяционного полинома $$ f(x)=A_x^+\dots+A_ $$ с величинами $ \sigma_<> $ и $ \tau_<> $: $$ \tau_0=A_0,\ \tau_k=A_0\sigma_k+A_1\sigma_+\dots+A_\sigma_1 + A_k \ npu \ k\in \ \ . $$

Формулы позволяют рекурсивно, начиная со старшего, вычислить коэффициенты интерполяционного полинома по величинам $ \sigma_<> $ и $ \tau_<> $.

Пример. Найти корни интерполяционного полинома, заданного таблицей

$$ \begin x & 1 & 2 & 3 & 4 & 5 \\ \hline y & 1 & -2 & 33 & 166 & 481 \end $$

Ответ. $ \sqrt>,\, — \sqrt \pm \sqrt> $.

Величины $ \sigma_k $,рассматриваемые как функции $ x_1,\dots,x_n $, известны в классической алгебре XIX века под названием алеффункции. Несмотря на их рациональный вид, они являются полиномами (симметрическими) по $ x_1,\dots,x_n $. Подробности ☞ ЗДЕСЬ.

Можно организовать и альтернативную схему вычисления коэффициентов интерполяционного полинома, начав с его свободного члена. Так, при условии $ \_^n $, имеем: $$ A_=(-1)^ \tau_ \prod_^n x_ \quad npu \quad \tau_= \sum_^ y_j \frac(x_j)> \, . $$

Интерполяционный полином в форме Ньютона

Основной недостаток построения интерполяционного полинома по методу (в форме) Лагранжа заключается в том, что при добавлении в таблицу нового узла (новых результатов измерений), в формуле приходится пересчитывать все слагаемые. От этого недостатка свободен метод Ньютона, в котором добавление нового узла ведет к добавлению лишь одного слагаемого к построенному ранее полиному.

Вывод представления этого полинома с помощью аппарата определителей ☞ ЗДЕСЬ. А само представление имеет вид $$ f(x)=B_1+(x-x_1)B_2+(x-x_1)(x-x_2)B_3 + \dots + (x-x_1)\times \dots \times (x-x_)B_n \ ; $$ он и называется интерполяционным полиномом в форме Ньютона.

Коэффициенты $ B_,B_2,\dots,B_n $ в этой форме определяются с помощью интерполяционной таблицы: $$ \begin y_1 & = & B_1, \\ y_2 & = & B_1 +(x_2-x_1)B_2, \\ y_3 & = & B_1 +(x_3-x_1)B_2+ (x_3-x_1)(x_3-x_2)B_3, \\ \dots & & \dots \\ y_n&= & B_1+(x_n-x_1)B_2+(x_n-x_1)(x_n-x_2)B_3 + \dots + (x_n-x_1)\times \dots \times (x_n-x_)B_n \end $$ Получаем последовательно: $$ B_1 = y_1,\ B_2= \frac,\ B_3=\frac-\frac>,\dots $$ Схему вычисления коэффициентов можно оформить в виде таблицы, если ввести следующее обозначение. Выражение $$ [x_j,x_k]=\frac \quad npu \quad j,k \in \, j\ne k $$ называется разделенной разностью первого порядка. Выражение $$[x_,x_,\dots,x_,x_j]=\frac<[x_,x_,\dots,x_]-[x_,\dots,x_,x_j]> $$ называется разделенной разностью порядка $ \mathbf p $.

Теорема. Интерполяционный полином в форме Ньютона записывается в виде:

$$ f(x)=y_1+[x_2,x_1](x-x_1)+[x_3,x_2,x_1](x-x_1)(x-x_2)+\dots+ $$ $$ +[x_n,x_,\dots,x_1] (x-x_1)\times \dots \times (x-x_) $$

Поясним схему вычисления разделенных разностей для случая $ n_<>=5 $. В первых двух столбцах стоят данные интерполяционной таблицы (со вставленными пустыми строками между соседними): $$ \begin x & y & \\ \hline x_1 & y_1 & \\ & & [x_2,x_1] \\ x_2 & y_2 & & [x_3,x_2,x_1] & \\ & & [x_3,x_2] & & [x_4,x_3,x_2,x_1] \\ x_3 & y_3 & & [x_4,x_3,x_2] & & [x_5,x_4,x_3,x_2,x_1] \\ & & [x_4,x_3] & & [x_5,x_4,x_3,x_2] \\ x_4 & y_4 & & [x_5,x_4,x_3] \\ & & [x_5,x_4] \\ x_5 & y_5 \end $$ Для вычисления третьего столбца мы вычитаем соседние значения $ y_<> $ и делим на соответствующую разность $ x_<> $: $$ [x_2,x_1] = \frac,\ [x_3,x_2]= \frac,\dots $$ и ставим получившиеся числа между соответствующими значениями $ x_<> $. Четвертый столбец заполняется аналогично: вычитаются соседние числа третьего столбца и делятся на разность крайних задействованных значений $ x_<> $: $$ [x_3,x_2,x_1] = \frac,\ [x_4,x_3,x_2] = \frac, \dots $$ И так далее. После получения треугольной таблицы выбираются числа из верхней стороны треугольника: именно они и являются искомыми коэффициентами интерполяционного полинома в форме Ньютона.

Пример. Построить интерполяционный полином по таблице

$$ \begin x & 1 & -2 & 3 & 0 & -1 \\ \hline y & 2 & 17 & 82 & 1 & 2 \end $$

Решение. $$ \begin x & y & \\ \hline 1 & \underline & \\ & & \underline \\ -2 & 17 & & \underline & \\ & & 13 & & \underline \\ 3 & 82 & & 7 & & \underline \\ & & 27 & & 0 \\ 0 & 1 & & 7 \\ & & -1 \\ -1 & 2 \end $$

Ответ. $$ f(x)=2 -5(x-1)+9(x-1)(x+2)+2(x-1)(x+2)(x-3)+(x-1)(x+2)(x-3)x \equiv x^4+1 \ . $$

Вычисления особенно упрощаются в случае равноотстоящих узлов интерполяции. Если $$ x_2-x_1=x_3-x_2=\dots=x_n-x_=h , $$ то $$ f(x)=y_1+\Delta y_1 \frac+\Delta^2 y_1 \frac+\dots+ \Delta^ y_1 \frac<(x-x_1)(x-x_2)\times \dots \times (x-x_)><(n-1)!h^> $$ где $$ \Delta^k y_1=y_-C_^1y_+ C_^2y_+\dots+(-1)^y_1 \ , $$ и $ C_^j $ означает биномиальный коэффициент.

Предельным переходом при $ h_<> \to 0 $ из этой формулы может быть получена формула Тейлора.

Применение полиномиальной интерполяции в задаче о разделении секрета

Обратная интерполяция

В одном из предшествующих ☝ пунктов решалcя следующий

Пример. Найти корни интерполяционного полинома, заданного таблицей

$$ \begin x & 1 & 2 & 3 & 4 & 5 \\ \hline y & 1 & -2 & 33 & 166 & 481 \end $$

Способ решения этой задачи предлагался стандартный: сначала строился интерполяционный полином $ f_<>(x) $, потом искались его корни. Можно, однако, попробовать идти следующим обходным путем. Именно, исходя из заданной таблицы попытаться восстановить обратную функцию для неизвестной, т.е. функцию $ g_<>(y) $ такую, что $ g(f(x)) \equiv x $ или же, эквивалентно, $ f(g(y)) \equiv y $. Разумеется, точно эту функцию — для полинома $ f_<>(x) $ — построить крайне сложно 2) , но ведь и сам интерполяционный полином является лишь некоторым приближением реальности! Имеет смысл решить поставленную задачу непосредственным использованием результатов таблицы, которая задает соответствие $ x\leftrightarrow y_<> $ для некоторых пар значений; это соответствие можно использовать не только в направлении $ x\rightarrow y_<> $, но и в направлении $ x \leftarrow y_<> $, т.е. строить интерполяционный полином для $ x_<> $ как функции от $ y_<> $ (фактически, «перевернуть» интерполяционную таблицу).

Для рассматриваемого примера форма Лагранжа дает представление полинома в виде: $$ g(y) =1 \times \frac+ 2\times \frac+ $$ $$ +3\times \frac+4\times \frac+ $$ $$ +5\times \frac \ . $$ Чтобы найти какому значению $ x_<> $ соответствует значение $ y_<>=0 $ мы просто подставляем это последнее в формулу и получаем $$ g(0) =\frac \approx 1.304553 \ . $$ Сравнивая этот ответ с полученным ☝ ВЫШЕ, мы видим некоторое соответствие с корнем $ \sqrt> \approx 1.126032 $. ♦

Понятно, что исключительным случаем для обратной интерполяции будет случай «коллизии» 3) : двум разным узлам интерполяции соответствует одно и то же значение $ y_<> $.

Ряд проведенных численных экспериментов показал, что бывают очень сильные расхождения в результатах «прямого» и обратного интерполирования — например, при сравнении корней $ f(x)=c $ со значением $ g(c_<>) $. Кажется, что результаты более адекватны реальности если величины аргументов и значений интерполируемой функции не слишком различаются по величине (хотя это еще надо проверять). Еще один пример см. в пункте Замена переменных и обратное отображение.

Интерполяционный полином Эрмита

Задача [Эрмит]. Построить полином $ y_<>=f(x) $, имеющий заданные значения своих производных в узлах интерполяции $ \_^s $: $$ \begin P(x_1), & \dots &, P^(x_1), \\ P(x_2), & \dots &, P^(x_2), \\ \dots ,& & \dots, \\ P(x_s), & \dots & , \ P^(x_s) . \\ \end $$ При $ m_j=1 $ узел $ x_ $ называется простым узлом интерполяции, при $ m_j>1 $ узел $ x_ $ называется кратным узлом.

Для случая вещественной интерполяционной таблицы ( $ \mathbb A = \mathbb R_<> $ ) задаче можно придать следующую геометрическую интерпретацию: требуется провести алгебраическую кривую $ y_<>=f_<>(x) $ через заданные точки $ (x_, y_j) $ так, чтобы в каждой точке обеспечить заданные наклоны касательных (а также, возможно, кривизны и т.п.).

Интерполяционный полином Эрмита используется в задаче разложения дробно-рациональной функции на простейшие дроби. Сам Эрмит применял его для оценки величины определенного интеграла по значениям функции и ее производных на концах интервала. Еще одно приложение — в задаче вычисления функции от матрицы.

Интерполяционная таблица дает $ m_+\dots + m_ $ условий на коэффициенты неизвестного полинома. По аналогии со стандартной задачей полиномиальной интерполяции, можно ожидать, что искомый полином будет существовать среди полиномов степени $ \le m_+\dots + m_s -1 $. Будем искать этот полином методом неопределенных коэффициентов. Обозначим $$ \widehat(x) = (x-x_1)^\times \dots \times (x-x_s)^ \ , $$ $$ \widehat(x) = \frac<\widehat(x)>> =(x-x_1)^\times \dots \times (x-x_)^ (x-x_)^ \times \dots \times (x-x_s)^ . $$ Пусть $ f_<>(x) $ — произвольный полином степени $ \le m_1+\dots + m_ -1 $. Разложим дробь $ f_<>(x)\big/\widehat(x) $ на сумму простейших над множеством $ \mathbb A_<> $: $$ \begin \displaystyle \frac<\widehat(x)> &\equiv & \displaystyle \frac>+\frac>+ \dots+ \frac><(x-x_1)^>+ & \\ \\ &+& \displaystyle \frac>+\frac>+ \dots+ \frac>>+ & \\ \\ &+&\dots + & \\ \\ &+& \displaystyle \frac>+\frac>+ \dots+ \frac><(x-x_s)^> & \\ && & \displaystyle = \sum_^ \sum_^<<\mathfrak m>_j>\frac><(x-\lambda_j)^<\ell>> \end $$ Определим числители дробей $ A_^<> $ с помощью интерполяционных данных. Домножим обе части тождества на $ (x-x_)^ $, получим: $$ \frac(x)> \equiv A_(x-x_1)^+A_(x-x_1)^+ \dots+ A_ + (x-x_1)^ G(x) \ ; $$ здесь через $ G_<>(x) $ обозначена дробно-рациональная функция по $ x_<> $, знаменатель которой не обращается в нуль при $ x=x_ $. Подставим это значение $ x_<> $ в обе части последнего равенства: $$ A_ = \frac(x_1)> \ . $$ Теперь продифференцируем последнее тождество по $ x_<> $, подставим $ x=x_ $ и воспользуемся данными интерполяционной таблицы: $$ A_ = \frac \left[ \frac(x)> \right] \Bigg|_> = \frac(x_1)\widehat(x_1)-P(x_1)\widehat^<\prime>(x_1)><\left[\widehat(x_1)\right]^2> \ . $$ Снова продифференцируем по $ x_<> $ и т.д. В результате получаем: $$ A_=\frac\frac \left[ \frac(x)> \right] \Bigg|_> \ . $$ Аналогично поступаем и с другими узлами интерполяции. В результате, получаем решение задачи в виде интерполяционного полинома Эрмита: $$ f(x)=\sum_^s \left[ A_(x-x_j)^+A_(x-x_j)^+ \dots+ A_ \right] \widehat(x) \ . $$

В литературе имеется неоднозначность терминологии — этот же полином называется и интерполяционным полиномом Лагранжа-Сильвестра, и интерполяционным полиномом Серре.

Теорема. Подмножество всевозможных полиномов из $ \mathbb A[x] $, принимающих значения по таблице, можно представить в виде

$$ \< f(x)+\widehat(x)q(x) \ \big| \ q_<>(x) \in \mathbb A[x] \> \ ; $$ здесь $ f_<>(x) $ — интерполяционный полином Эрмита.

Интерполяционный полином Эрмита является естественным обобщением обычного интерполяционного полинома в форме Лагранжа ( $ s=n, m_=1,\dots,m_s=1 $) и формулы Тейлора ( $ s=1, m_=n $).

Можно указать и явное представление для этого полинома — с использованием формализма определителей. На основании правила дифференцирования дробно-рациональной функции, получаем: $$ f(x) = $$ $$ =- \sum_^s \frac(x)> <\left[ \widehat(x_j)\right]^<\mathfrak m_j>> \left|\begin \widehat (x_j) & 0 & 0 & \dots & 0 & P(x_j) \\ \widehat^ <\prime>(x_j) & \widehat (x_j) & 0 & \dots & 0 & P^<\prime>(x_j) \\ \widehat^<\prime \prime>(x_j) & 2 \widehat^<\prime>(x_j) & \widehat(x_j) & \dots & 0 & P^<\prime \prime>(x_j) \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \widehat^<[\mathfrak m_j-1]>(x_j) & C_<\mathfrak m_j-1>^1 \widehat^<[\mathfrak m_j-2]>(x_j) & C_<\mathfrak m_j-1>^2 \widehat^<[\mathfrak m_j-3]>(x_j) & \dots & \widehat (x_j) & f^<[\mathfrak m_j-1]>(x_j) \\ 1 & \displaystyle> & \displaystyle> & \dots & \displaystyle <\frac<(x-x_j)^<\mathfrak m_j - 1>><(\mathfrak m_j-1)!>>& 0 \end \right|_ <(\mathfrak m_j+1)\times (\mathfrak m_j+1)>\ . $$ Здесь $ C_N^M $ — биномиальный коэффициент. Еще один подход к построению полинома см. ☞ [8].

Пример. Построить интерполяционный полином по таблице

$$ \begin x & -1 & 0 & 1 & 2 \\ \hline f(x) & 16 & 7 & 8 & 217 \\ \hline f^<\prime>(x) & & -1 & -4 & 1375 \\ \hline f^<\prime \prime>(x) & & 6 & -44 & \\ \hline f^<\prime \prime \prime>(x) & & & -126 & \end $$

Решение. Здесь $$ \widehat(x) = (x+1)x^3(x-1)^4(x-2)^2 \ . $$ Для $ x_<>=-1 $ имеем $ m_=1 $ и в формуле Эрмита этому узлу соответствует одно слагаемое: $$ 16 \frac \ . $$ Для $ x_<>=0 $ имеем $ m_=3 $ и этому узлу соответствует полином $$ \left[ A_x^+A_x+ A_ \right] (x+1)(x-1)^4(x-2)^2 \ , $$ значения которого — вместе с первой и второй производной — в точке $ x_<>=0 $ должны совпадать с табличными: $$ A_ = \frac=\frac \ ; $$ $$ 4\,A_-16 A_ = -1 \quad \Rightarrow \quad A_= 27/4 \ ; $$ $$ 8\,A_-32\,A_+42\,A_= 6 \quad \Rightarrow \quad A_ = 297/16 \ . $$ Для $ x_<>=1 $ имеем $ m_=4 $: $$ \left[ A_(x-1)^+A_(x-1)^2+ A_(x-1) +A_ \right] \underbrace_(x)> \ , $$ и для определения четырех коэффициентов этого полинома мы имеем четыре условия из таблицы: $$ 2\, A_=8 \quad \Rightarrow \quad A_=4 \ , $$ $$ 2\, A_+3\, A_=-4 \quad \Rightarrow \quad A_=-8 \ , $$ $$ 4\,A_+6\,A_-6\,A_= -44 \quad \Rightarrow \quad A_=7 \ , $$ $$ 12\,A_+18\,A_-18\,A_-36\,A_ =-126 \quad \Rightarrow \quad A_=-21 . $$ Этот же результат можно получить и в виде альтернативного — детерминантного — представления: $$ \widehat(1)=2,\ \widehat^<\prime>(1)=3,\ \widehat^<\prime \prime>(1)=-6,\ \widehat^<\prime \prime \prime>(1)=-36, $$ и $$ -\frac<\widehat(x)><[\widehat(1)]^4> \left| \begin \widehat(1) & 0 & 0 & 0 & f(1) \\ \widehat^<\prime>(1) & \widehat(1) & 0 & 0 & f^<\prime>(1) \\ \widehat^<\prime \prime>(1) & 2 \, \widehat^<\prime>(1) & \widehat(1) & 0 & f^<\prime \prime>(1) \\ \widehat^<\prime \prime \prime>(1) & 3\, \widehat^<\prime \prime>(1) & 3\, \widehat^<\prime>(1) & \widehat(1) & f^<\prime \prime \prime>(1) \\ 1 & x-1 & \frac(x-1)^2 & \frac(x-1)^3 & 0 \end \right|= $$ $$ =- \frac<\widehat(x)> \left| \begin 2 & 0 & 0 & 0 & 8 \\ 3 & 2 & 0 & 0 & -4 \\ -6 & 6 & 2 & 0 & -44 \\ -36 & -18 & 9 & 2 & -126 \\ 1 & x-1 & \frac(x-1)^2 & \frac(x-1)^3 & 0 \end \right| \ . $$

Наконец, для $ x_<>=2 $ имеем $ m_=2 $: $$ \left[\frac(x-2)+ \frac \right](x+1)x^3(x-1)^4 \ . $$

Ответ. $ 2\,x^-3\,x^8-4\,x^5+5\,x^4-x^3+3\,x^2-x+7 $.

Построить уравнение «горки»: найти полином из условий

Рациональная интерполяция

Задача. Для таблицы значений $$ \begin x & x_1 & x_2 & \dots & x_N \\ \hline y & y_1 & y_2 &\dots & y_N \end $$ построить рациональную функцию в виде $ f(x)=p(x)/q(x) $ при $ p_<>(x) $ — полиноме степени $ \le n $, $ q_<>(x) $ — полиноме степени $ \le m $, так, чтобы $ \< f(x_j)=y_j \>_^N $. При этом $ N,n,m $ связаны соотношением: $$ N=m+n+1 \, . $$

Общее число неопределенных коэффициентов полиномов $ p_<> $ и $ q_<> $ равно $ n+m+2 $; но мы не будем различать решения задачи отличающиеся на общий числовой множитель числителя и знаменателя.

Первое решение задачи было предложено Коши в 1821 г. [9].

Теорема [Коши]. Обозначим:

При $ m=0 $ получаем представление интерполяционного полинома в форме Лагранжа.

Пример. $ n =1, m=2, N=4 $

Биографические заметки о Коши ☞ ЗДЕСЬ.

Другие примеры на применение теоремы Коши ☞ ЗДЕСЬ

Теорема Коши дает решение задачи в смысле «как правило». Дело в том, что задача рациональной интерполяции (в указанной постановке) не всегда разрешима.

Пример. Для таблицы

$$ \begin x & -1 & 0 & 1 & 2 & 3 \\ \hline y & 1 & 1 & 1/3 & 3 & 1/13 \end $$ построить удовлетворяющую ей рациональную функцию $ f(x)=p(x)/q(x) $ при $ \deg p(x)=1,\deg q(x) =3 $.

Решение. Теорема Коши дает решение в виде $$ p(x)\equiv x-2, \quad q(x)\equiv x^3-x^2-x-2 \, . $$ Однако $ p(2)=0 $ и $ q(2)=0 $, и условие $ f(2)=3 $ не выполняется. Оно не выполнится даже если мы разделим полученные полиномы на общий линейный делитель.

Вместе с тем, той же таблице удовлетворяют (безо всяких проблем) рациональные дроби $$ \frac \ , \quad \frac, \dots $$ для которых нарушено поставленное ограничение $ N=m+n+1 $, связывающее степени числителя и знаменателя с количеством узлов интерполяции. ♦

Альтернативный подход к решению задачи основывается на следующей теореме, развивающей результат К.Якоби [10,11]; он основан на идее из пункта ☝ о рекурсивном вычислении коэффициентов интерполяционного полинома.

Теорема. Пусть $ \left\< y_j \ne 0 \right\>_^N $. Вычислим последовательности

Подробнее о методе Якоби (в том числе и об эффективном способе вычисления ганкелевых полиномов) ☞ ЗДЕСЬ.

Тригонометрическая интерполяция

Потребность в подобной интерполяции возникает в случае, когда приближаемая функция по своей природе предполагается периодической с известным периодом, например $ 2 \pi_<> $.

Задача. Построить периодическую функцию периода $ 2 \pi_<> $, принимающую значения согласно следующей таблице: $$ \begin x & x_1 & x_2 & \dots & x_n \\ \hline y & y_1 & y_2 &\dots & y_n \end $$ Здесь узлы интерполяции предполагаются различными вещественными $ \,\dots,x_n\> \subset [0,2 \pi] $ (хотя все приведенные в этом пункте результаты будут справедливы и для комплексных узлов при различных $ \ < \mathfrak(x_k)\>_^ \subset [0,2\pi[ $).

Теорема [Эрмит ?]. Функция

Доказательство тривиально, если обратить внимание на аналогию с интерполяционным полиномом в форме Лагранжа. ♦

Проанализируем полученный ответ, сначала преобразовав его. Каждое произведение в числителе можно представить в виде комбинации синусов и косинусов от аргументов кратных $ x_<> $. Так, например, $$ \sin(x-x_1) \sin (x-x_2) \equiv \frac \left( \cos (x_2-x_1) — \cos(2\,x-x_1-x_2) \right)\equiv $$ $$ \equiv \frac \cos (x_2-x_1) — \frac \cos(x_1+x_2) \cos\, 2\,x — \frac \sin(x_1+x_2) \sin \,2\,x \ ; $$ $$ \sin(x-x_1) \sin (x-x_2) \sin (x-x_3) \equiv \frac \left( \cos (x_2-x_1) \sin (x-x_3) — \cos(2\,x-x_1-x_2) \sin (x-x_3) \right) \equiv $$ $$ \equiv \frac\left[\cos(x_1-x_2)\cos(x_3)+\cos(x_1-x_3)\cos \, x_2+\cos(x_2-x_3)\cos(x_1)\right]\sin\,x- $$ $$ — \frac\left[\cos(x_1-x_2)\sin\, x_3+\cos(x_1-x_3)\sin\, x_2+\cos(x_2-x_3)\sin\,x_1\right]\cos\,x — $$ $$ -\frac\cos (x_1+x_2+x_3) \sin 3\,x +\frac \sin(x_1+x_2+x_3)\cos\, 3\, x \ , $$ и т.д. Таким образом, функцию $ F_<>(x) $ можно представить в виде линейной комбинации тригонометрических функций $ \sin \,x, \cos \, x $, $ \sin 2\,x, \cos 2\,x , \dots, \sin \, (n-1) x,\ \cos \, (n-1) x $.

Функция от $ x_<> $ вида $$ \begin f_n(x)=a_0 & + & a_1 \cos x + a_2 \cos 2x+\dots + a_n \cos nx + \\ \ &+&b_1 \sin x + b_2 \sin 2x+\dots + b_n \sin nx . \end $$ при постоянных $ \ $ называется тригонометрическим полиномом по переменной $ x_<> $. Если хотя бы одно из чисел $ a_ $ или $ b_ $ отлично от нуля, то говорят, что тригонометрический полином имеет порядок $ \mathbf n $.

Тригонометрический полином порядка $ n_<> $ определяется заданием $ (2\, n+1) $-го своего коэффициента $ \ $. Возникает гипотеза, что, по аналогии со случаем полиномиальной интерполяции, задание $ (2\,n+1) $-го значения интерполяционной таблицы определит такой полином однозначно. Полином, построенный по теореме Эрмита, не удовлетворяет этому условию: указанные выше преобразования, проведенные для случая $ (2\, n+1) $-го узла интерполяции приведут к тригонометрическому полиному порядка $ 2\, n $. Тем не менее, идея, лежащая в основе той теоремы, допускает естественное обобщение.

Теорема. Функция

$$ \begin F(x)&=&y_1\frac\times\dots \times\sin \frac>>\times \dots \times \sin \frac>> + \\ &+&y_2\frac\sin\frac\times \dots \times \sin\frac>>\sin\frac\times \dots \times \sin\frac>> + \\ & + & \dots + \\ &+& y_\frac \times \dots \times \sin \frac>>< \displaystyle \sin \frac \times \dots \times \sin \frac> \end $$ удовлетворяет условиям $ \< F(x_j)=y_j \>_^ $. Эта функция является тригонометрическим полиномом порядка не выше $ n_<> $. При таком ограничении на порядок полином, удовлетворяющий заданным интерполяционным условиям, определяется единственным образом.

Задача. Найти явные выражения для коэффициентов тригонометрического полинома из последней теоремы.

«Лобовое» решение аналогично решению задачи полиномиальной интерполяции — сведением ее к подходящей системе линейных уравнений.

Пример. Построить интерполяционный полином второго порядка по следующей таблице

$$ \begin x & 0 & 2\pi/5 & 4\pi/5 & 6\pi/5 & 8\pi/5 \\ \hline y & 1 & -1 & 3 & -2 & 2 \end $$

Решение. Будем искать коэффициенты полинома традиционным методом неопределенных коэффициентов: полином $$ a_0 + a_1 \cos x + b_1 \sin x + a_2 \cos 2\, x + b_2 \sin 2\, x $$ должен удовлетворять таблице. Получаем систему линейных уравнений, которую перепишем в матричном виде: $$ \left( \begin 1 & 1 & 0 & 1 & 0 \\ 1 & \cos 2\pi/5 & \sin 2\pi/5 & \cos 4\pi/5 & \sin 4\pi/5 \\ 1 & \cos 4\pi/5 & \sin 4\pi/5 & \cos 8\pi/5 & \sin 8\pi/5 \\ 1 & \cos 6\pi/5 & \sin 6\pi/5 & \cos 12\pi/5 & \sin 12\pi/5 \\ 1 & \cos 8\pi/5 & \sin 8\pi/5 & \cos 16\pi/5 & \sin 16\pi/5 \end \right) \left( \begin a_0 \\ a_1 \\ b_1 \\ a_2 \\ b_2 \end \right) = \left( \begin 1 \\ -1 \\ 3 \\ -2 \\ 2 \end \right) \ . $$ Согласно теореме эта система должна иметь решение относительно $ a_0,a_1,b_1,a_2,b_2 $, причем единственное. Так оно и оказывается: на акцентируя внимания на методе решения этой системы (см. ☞ ЗДЕСЬ ) , а также на том, как можно выразить в виде квадратичных иррациональностей косинусы и синусы углов кратных $ \pi/5 $ (см. ☞ ЗДЕСЬ ), приведу лишь конечный результат — $$ a_0 =3/5,\ a_1 =1/5,\ a_2 =1/5,\ $$ $$ b_1=-\frac\sqrt<2(5+\sqrt)>(11-5\sqrt)\approx 0.034302,\ b_2= -\frac\sqrt<2(5+\sqrt)>(7+3\sqrt)\approx -2.607455 \ . $$ ♦

$$ \begin x =\frac\pi \times & \scriptstyle 1 & \scriptstyle 2 & \scriptstyle 3 & \scriptstyle 4 & \scriptstyle 5 & \scriptstyle 6 & \scriptstyle 7 & \scriptstyle 8 & \scriptstyle 9 & \scriptstyle 10 & \scriptstyle 11 & \scriptstyle 12 & \scriptstyle 13 & \scriptstyle 14 & \scriptstyle 15 & \scriptstyle 16 & \scriptstyle 17 & \scriptstyle 18 & \scriptstyle 19 & \scriptstyle 20 & \scriptstyle 21 \\ \hline y & \scriptstyle 175 & \scriptstyle 196& \scriptstyle 230& \scriptstyle 253& \scriptstyle 245& \scriptstyle 205& \scriptstyle 135& \scriptstyle 100& \scriptstyle 82& \scriptstyle 85& \scriptstyle 82& \scriptstyle 64& \scriptstyle 29& \scriptstyle 10& \scriptstyle 15& \scriptstyle 50& \scriptstyle 110& \scriptstyle 158& \scriptstyle 174& \scriptstyle 173& \scriptstyle 175 \end $$ построить (алгебраический) полином $ 20 $-й степени и тригонометрический полином $ 10_<> $-го порядка и сравнить их значения при $ x=\pi/7 $ и $ x=\pi/5 $.

Решение ☞ ЗДЕСЬ.

Для случая системы равноотстоящих узлов решение задачи значительно упрощается.

Задача. Построить полином $ f_ (x) $ такой, что $$ f_n\left( \frac<2\pi j> \right) = y_j \quad npu \quad j \in\ \ . $$

Теорема. Коэффициенты полинома $ f_ (x) $ определяются формулами

$$ a_0= \frac \sum_^n y_j, $$ $$ a_m= \frac \sum_^n y_j \cos \left( \frac<2\pi jm> \right),\quad b_m= \frac \sum_^n y_j \sin \left(\frac<2\pi jm> \right) $$ при $ m\in \ $.

При бесконечном увеличении числа узлов интерполяции $ n_<> $ в пределе получаем коэффициенты ряда Фурье для функции $ y_<>=f(x) $: $$ a_0= \frac<2> <\pi>\int_<-\pi>^ <\pi>f(x) d\, x ,\ $$ $$ a_m = \frac <\pi>\int_<-\pi>^ <\pi>f(x) \cos mx d\, x ,\ b_m = \frac <\pi>\int_<-\pi>^ <\pi>f(x) \sin mx d\, x \ \mbox < при >\ m\in \mathbb N . $$

Подробное изложение теории тригонометрической интерполяции (и дискретного преобразования Фурье) ☞ ЗДЕСЬ

Интерполяция суммами экспонент

Если Вы решили упражнение из предыдущего пункта, то могли обнаружить что интерполяционный алгебраический полином $ A_0+A_1x+A_2x^2+\dots $ не всегда является подходящим средством приближения неизвестной функции, особенно если та, по своей «скрытой природе», принципиально неполиномиальна. В последней фразе слово «приближение» пока не формализовано — оно связано с понятием аппроксимации, которое дается ☟ НИЖЕ. Забегая вперед, заметим, что алгебраический полином, построенный в задаче интерполяции, подходит не для всех задач,связанных с приближением неизвестной функции. Так, интуитивно понятно, почему для приближения периодической функции в предыдущем пункте использовались периодические функции с соизмеримыми периодами.

Другой класс таких «специфических» функций составляют функции, имеющие определенную асимптотику при неограниченном возрастании их аргументов (по абсолютной величине или в определенном направлении). В этом случае ставят задачу приближения функции $ F_<>(x) $ суммой экспонент: $$ f(x) = c_1e^+c_2e^+\dots+c_ne^ \ ; $$ как видим, эта функция задается $ 2\,n $ параметрами — коэффициентами $ \\>_^n $ и показателями $ \_^n $.

Пример. Пусть имеется вещество, представляющее из себя смесь двух радиоактивных материалов с различными периодами полураспада. По закону радиоактивного распада изменение количества $ N(t) $ радиоактивного материала (измеряемого в граммах или в количестве атомов) во времени вычисляется по формуле

$$ N(t)=N(0) e^ <-\lambda t>\ , $$ где $ N(0) $ — исходное количество материала, а $ \lambda_<> > 0 $ — постоянная распада, различная для различных материалов. Тогда для смеси двух материалов получаем зависимость $$ N_(t)=N_1(0) e^<-\lambda_1 t>+ N_2(0) e^ \ . $$ Требуется определить параметры процесса $ N_1(0), N_2(0), \lambda_1, \lambda_2 $ по показаниям счетчика Гейгера в различные моменты времени. ♦

См. также моделирование отталкивающего потенциала при описании межатомного взаимодействия ☞ Ziegler, Biersack, Littmark (ZBL) repulsive potential.

Задача. Построить функцию $ \displaystyle f(x) = \sum_^n c_k e^ $, удовлетворяющую условиям $$ \_^ ; $$ узлы интерполяции $ \_^ $ предполагаются различными.

В отличие от рассмотренных выше задач алгебраической или тригонометрической интерполяции, поставленная задача является, во-первых, принципиально нелинейной относительно параметров, и, во-вторых, не всегда разрешимой.

Пример. Условиям $ f(x_0)=0, f(x_1)=1 $ функция $ f(x)=ce^ $ не удовлетворяет ни при каких значениях параметров. ♦

Следующий метод решения задачи для случая вещественных равноотстоящих узлов $$ x_j= j h \quad \mbox \quad h>0,\ j\in\ \ . $$ был предложен Прони 4) в статье [1].

Обозначим $$ u_1=e^,\dots, u_n=e^ ; $$ в этих обозначениях условия интерполяции можно переписать в виде: $$\+c_2u_2^+\dots+ c_nu_n^=y_j \>_^ \ . $$ Задача свелась к конечной проблеме моментов, из решения которой и вытаскивается алгоритм

Итак, получившаяся система из $ 2\, n $ уравнений включает в себя $ 2\, n $ неопределенных параметров $ \\>_^n $ и $ \_^n $. Вхождение первых — хорошее, поскольку линейное. Со вторыми приходится обращаться с помощью «обходного манёвра».

1. Составляется система линейных уравнений, ее матричный вид: $$ \left(\begin y_0 & y_1 & y_2 & \dots & y_ \\ y_1 & y_2 & y_3 & \dots & y_n \\ y_2 & y_3 & y_4 & \dots & y_ \\ \vdots & & & & \vdots \\ y_ & y_n & y_ & \dots & y_ \end\right) \left( \begin a_n \\ a_ \\ a_ \\ \vdots \\ a_1 \end \right)= \left( \begin y_n \\ y_ \\ y_ \\ \vdots \\ y_ \end \right) $$ здесь неизвестными являются 5) $ a_1,a_2,\dots,a_ $. На числа $ \_^ $ накладываются два ограничения: $$ \left|\begin y_0 & y_1 & y_2 & \dots & y_ \\ y_1 & y_2 & y_3 & \dots & y_n \\ y_2 & y_3 & y_4 & \dots & y_ \\ \vdots & & & & \vdots \\ y_ & y_n & y_ & \dots & y_ \end\right| \ne 0,\ \left|\begin y_1 & y_2 & y_3 & \dots & y_ \\ y_2 & y_3 & y_4 &\dots & y_ \\ y_3 & y_4 & y_5 &\dots & y_ \\ \vdots & & & & \vdots \\ y_n & y_ & y_ & \dots & y_ \end\right| \ne 0 \ . $$ эти два определителя относятся к типу ганкелевых.

2. При выполнении этих ограничений система линейных уравнений имеет единственное решение $$ a_,a_,\dots,a_1 \ , $$ у которого $ a_ \ne 0 $ (см. ☞ Формулы Крамера ). Находим это решение.

3. Составляем алгебраическое уравнение относительно новой переменной $ \lambda $: $$ g(\lambda)= \lambda^n — a_1 \lambda^-a_2 \lambda^-\dots-a_n = 0 \ , $$ коэффициенты которого найдены в предыдущем пункте. Проверяем это уравнение на отсутствие кратных корней (см. методы ☞ ЗДЕСЬ или ☞ ЗДЕСЬ ). Находим все корни (в общем случае, они — комплексные и, как правило, могут быть найдены разве лишь приближенно). Эти корни и будут искомыми числами $ u_1,u_2,\dots,u_n $.

4. Оставшиеся параметры $ c_1,c_2,\dots,c_n $ определяем из системы линейных уравнений $$\+c_2u_2^+\dots+ c_nu_n^=y_j \>_^ \ . $$ В матричном виде: $$ \left(\begin 1 & 1 & \dots & 1 \\ u_1 & u_2 & \dots & u_n \\ u_1^2 & u_2^2 & \dots & u_n^2 \\ \vdots & & & \vdots \\ u_^ & u_^ & \dots & u_^ \end\right) \left( \begin c_1 \\ c_2 \\ c_3 \\ \vdots \\ c_n \end \right)= \left( \begin y_0 \\ y_ \\ y_ \\ \vdots \\ y_ \end \right) \ , $$ в которую подставлены величины $ u_1,u_2,\dots,u_n $, найденные в предыдущем пункте. Ее определитель — определитель Вандермонда, который, по предположению пункта 3 , отличен от нуля. Следовательно система имеет единственное решение.

5. Показатели $ \_^n $ вычисляются по формулам $$ m_j=(\operatorname u_j)/h \, . $$ Я не объясняю как вычисляется натуральный логарифм мнимого числа, но вот правило вычисления $ e^ $ можно посмотреть ☞ ЗДЕСЬ.

Пример. Построить экспоненциальную функцию вида

$$ f(x)=c_1e^+c_2e^+c_3e^+c_4e^ \, , $$ принимающую значения согласно следующей таблице:

$$ \begin x & 0 & 0.1 & 0.2 & 0.3 & 0.4 & 0.5 & 0.6 & 0.7 \\ \hline y & 0.5 & 0.378907 & 0.232694 & 0.088027 & -0.016793 & -0.028242 & 0.127344 & 0.550507 \end $$

Решение. Для данного примера $ h=0.1, n=4 $. Система из пункта 1 приведенного алгоритма имеет решение $$ a_1 \approx 4.462091,\ a_2 \approx — 7.447287, \ a_3\approx 5.521117, a_4 \approx -1.537257 \ . $$ Полином $$ g(\lambda)=\lambda^4-4.462091\, \lambda^3+7.447287\, \lambda^2-5.521117\, \lambda+1.537257 $$ из пункта 3 алгоритма имеет корни $$ u_1 \approx 1.127497,\ u_2 \approx 1.363425, u_ \approx 0.985585\pm 0.169182\, \mathbf i \ . $$ Составляем на основании этих значений систему линейных уравнений из пункта 4 алгоритма и решаем ее: $$ c_1 \approx -2.399999, c_2 \approx 0.600001, c_ \approx 1.149999\mp 0.0000004 \mathbf i \ . $$ Теперь осталось только перевести найденные величины $ \_^4 $. в искомые показатели $ \_^4 $ по формулам $ \ u_j)/h\>_^4 $: $$ m_1 \approx 1.200001, m_2 \approx 3.0999996, \ m_ \approx -0.0000003\pm 1.700000 \mathbf i \ . $$ Окончательно: $$ f(x)\approx -2.399999\, e^+0.600001\, e^+ 2.299999\, \cos 1.700000 + 0.0000008\, \mathbf i \, \sin 1.700000 \ , $$ при истинном ответе: $$ f(x) = -2.4\, e^+0.6\,e^+2.3\,\cos \, 1.7\,x \ . $$ ♦

Последний пример дает интересный «поворот мысли». Представленный алгоритм отработал правильно и обнаружил периодическую функцию, наличие которой в гипотезе постановки задачи изначально не предполагалось! Алгоритм оказался «умнее» постановщика задачи… Этот эффект является проявлением замечательного свойства функций комплексного аргумента: в терминах комплексных чисел свойства функций $ e^ $ практически совпадают со свойствами $ \cos z $ и $ \sin z $, поскольку две последние линейно выражаются через первую: $$ \cos z= \frace^ <\mathbf i z>+ \frac e^<-\mathbf i z>, \sin z = — \frac<\mathbf i>e^ <\mathbf i z>+ \frac<\mathbf i>e^ <-\mathbf i z>\ .$$ Возникает ужасное подозрение, что весь материал пункта о тригонометрической интерполяции может быть получен средствами настоящего пункта . — Так оно и есть. Однако универсальный подход оказывается более громоздким, и в частном случае когда приближаемая функция предполагается периодической он менее эффективен, чем частный алгоритм построения тригонометрического полинома. Другое дело, если интерполируемая функция имеет вид: $$ \sin x + 3 \cos \sqrt x $$ — т.е. представляет комбинацию двух периодических функции, но с несоизмеримыми периодами! В этом случае для ее идентификации придется применять «тяжелую артиллерию» в виде комбинации экспонент…

Аппроксимация

Задача интерполяции является частным случаем более общей задачи аппроксимации функций, т.е. замены одной неизвестной или сложной для вычисления функции другой, более простой. Здесь существенно понятие «близости» функций, которое может быть различным в конкретных задачах.

Можно определить «расстояние» между функциями $ f_<>(x) $ и $ g_<>(x) $, заданными на отрезке $ [a_<>,b] $, как $$ \max_ |f(x) -g(x)| \ , $$ т.е. как величину максимального «разноса» ординат графиков этих функций. С точки зрения этого определения, интерполяционный полином не всегда дает подходящий способ аппроксимации. Это показывает и упражненение из пункта о тригонометрической интерполяции и следующий пример.

Пример [Рунге]. Пусть на интервале $ [-1,1] $ задана функция

$$F(x)=\frac \ . $$ Построим для нее интерполяционный полином $ f_(x) $ степени $ n_<> $ по равноотстоящим узлам $$ x_=-1+\frac\ npu \ k \in \ \ , $$ и будем оценивать $$ \max_ |f_n(x) -F(x)| \ . $$

Почему это происходит? Вопрос оказывается связанным со сходимостью ряда. Интерполяционный ряд $$ f_1+(f_2-f_1)+(f_3-f_2) + \dots $$ будет сходиться в некоторой окрестности нуля и расходиться при значениях $ x_<> $ близких к краю рассмотренного интервала. Для приведенного примера я, правда, оценок не нашел, но вот для другого примера Рунге $$ F(x)=\frac $$ в [3] указывается 6) , что интерполяционный ряд построенный для равноотстоящих на интервале $ [-5,5_<>] $ узлов сходится при $ x_<> \in ] -3.63,\ 3.63 [ $ и расходится вне этого интервала. ♦

Еще одна иллюстрация проблемы — аппроксимация ступенчатой функции

Можно определить «расстояние» между функциями $ f_<>(x) $ и $ g(x_<>) $, непрерывными на отрезке $ [a_<>,b] $, как $$ \int_a^b |f(x) -g(x)| d\,x \ , $$ или же $$ \int_a^b [f(x) -g(x)]^2 d\,x \ . $$ Для выяснения геометрического смысла этих определений, вспомним, что геометрический смысл интеграла $ \int_a^ f(x) d\, x $ от неотрицательной на $ [a_<>,b] $ функции $ f_<>(x) $ — площадь криволинейной трапеции на плоскости $ (x_<>,y) $, ограниченной прямыми $ x=a_<>,\,x=b,\,y=0 $ и графиком $ y=f(x_<>) $. Следовательно, интеграл $$ \int_a^b \left| f(x) -g(x) \right| d\, x $$ определяет площадь фигуры, ограниченной прямыми $ x=a_<>,\,x=b $ и графиками $ y=f(x_<>), y=g(x) $ (закрашена голубым на рисунке). Чем меньше эта площадь, тем «теснее» друг к другу расположены графики $ y=f_<>(x) $ и $ y=g_<>(x) $ на отрезке $ [a_<>,b] $. Величина $$ \int_a^b [f(x) -g(x)]^2 d\,x , $$ вообще говоря, не совпадает с предыдущей, но смысл ее тот же: она позволяет оценивать близость графиков на всем отрезке $ [a_<>,b] $. С вычислительной точки зрения, проще вычислять интеграл от квадрата функции, чем от ее модуля.

Задача. Для функции $ f(x_<>) $, непрерывной на $ [0_<>,1] $, построить полином $ f_(x_<>)=a_0+a_1x+\dots+a_x^ $ такой, чтобы величина $$ \int_0^1 \left[f(x)-f_(x) \right]^2 d\, x $$ была наименьшей.

Теорема. Коэффициенты полинома $ f_(x_<>) $ находятся из системы уравнений

$$ \left( \begin 1 & \frac & \frac & \dots & \frac \\ &&&& \\ \frac & \frac & \frac & \dots & \frac \\ &&&& \\ \frac & \frac & \frac & \dots & \frac \\ &&&& \\ \dots & & & & \dots \\ &&&& \\ \frac & \frac & \frac & \dots & \frac \end \right) \left( \begin a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_ \end \right)= \left( \begin b_0 \\ b_1 \\ b_2 \\ \vdots \\ b_ \end \right) \ . $$ Здесь $$ b_=\int_0^1 f(x) x^ d\, x \quad npu \quad j\in \ . $$

Матрица $ <\mathfrak H>_ $ системы называется матрицей Гильберта. Ее определитель равен $$ \frac <[1!\,2!\, 3! \dots (n-1)!]^3> \ , $$ т.е. отличен от нуля, и, следовательно, система имеет единственное решение.

Хотя определитель Гильберта и отличен от нуля, но очень близок к нему: $$ \det <\mathfrak H>_3=\frac,\ \det <\mathfrak H>_4=\frac,\ \det <\mathfrak H>_6=\frac $$ поэтому при решении линейной системы надо внимательно следить за ошибками округлений [4],[5].

Задача. Для функции $ f_<>(x) $, непрерывной на $ [-\pi,\pi_<>] $, построить тригонометрический полином $$ \begin f_n(x)=a_0 & + & a_1 \cos x + a_2 \cos 2x+\dots + a_n \cos nx + \\ \ &+&b_1 \sin x + b_2 \sin 2x+\dots + b_n \sin nx . \end $$ такой, чтобы величина $$ \int_<-\pi>^ <\pi>\left[f(x)-f_(x) \right]^2 d\, x $$ была наименьшей.

Теорема. Коэффициенты тригонометрического полинома $ f_(x_<>) $ являются коэффициентами Фурье функции $ f(x_<>) $:

$$ a_0=\frac <2\pi>\int_<-\pi>^ <\pi>f(x) d\, x \ , $$ $$ a_m = \frac <\pi>\int_<-\pi>^ <\pi>f(x) \cos mx d\, x ,\quad b_m = \frac <\pi>\int_<-\pi>^ <\pi>f(x) \sin mx d\, x \quad $$ при $ m\in \ $ .

Аппроксимация в случае недостоверности данных

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

Все числа в настоящем и следующем пунктах предполагаются вещественными.

Задача. Найти координаты точки $ (x, y) $, сумма квадратов расстояний до которой от точек $ P_1=(x_,y_1),\dots, P_n=(x_n,y_n) $ будет минимальной.

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

Теорема 1. Координаты искомой точки определяются средними значениями координат точек $ \< P_j\>_^n $:

Обозначение среднего значения числовой последовательности через $ \overline> $ является традиционным; не следует его путать с комплексным сопряжением. Точку $ (\overline,\overline) $ иногда называют центроидом для системы точек $ \< P_j\>_^n $.

Координаты точки, для которой величина

Иными словами: искомая точка совпадает с центром масс (барицентром) системы масс $ \_^n $, расположенных в точках $ \_^n $.

Задача. Найти уравнение прямой в виде $$ ax+by+1 =0 \ , $$ сумма квадратов расстояний до которой от точек $ \_^n $ будет минимальной.

Если обозначить $ \delta_j $ расстояние от $ P_j $ до некоторой прямой, то речь идет о поиске такой прямой, которая бы обеспечивала $$ \min (\delta_1^2+\delta_2^2+\dots + \delta_n^2) \, . $$

Теорема 2 [3],[4]. Обозначим

Геометрическая интерпретация. Линейное уравнение из теоремы означает, что искомая прямая обязательно проходит через центр тяжести системы точек. Уравнение для определение тангенса угла наклона этой прямой к оси $ Ox_<> $ — квадратное, оно имеет два корня $ \mu_, \mu_2 $, т.е. имеется два кандидата на решение задачи. По формулам Виета $ \mu_=-1/\mu_1 $, что означает: две прямые взаимно перпендикулярны. На рисунке эти прямые показаны для случая $$ \begin x & 0.5 & 1 & 1.5 & 2 & 2.5 & 3 \\ \hline y & 0.35 & 0.80 & 1.70 & 1.85 & 3.51 & 1.02 \end $$ Их уравнения: $$-1.59\, x +1.16\, y +1=0, \ -0.26\, x — 0.35\, y +1 =0 \ ; $$ они пересекаются в центроиде системы точек $ \ $: $$ \overline =1.75, \overline \approx 1.54 \ . $$ Сумма квадратов расстояний от данных точек до этих прямых: $$ 2.25 \ \mbox < и >\ 8.36 $$ соответственно. Следовательно, решением задачи является первая прямая (зеленая). ♦

Теорема 3 [Пирсон]. Пусть матрицы $ X $ и $ Y $ определены как в теореме $ 2 $. Тогда минимум суммы квадратов расстояний точек $ \<(x_,y_j)\>_^n $ до прямой

$$ ax+by+c=0 $$ равен минимальному положительному корню $ \lambda_ $ уранения $$ \left| \begin \sum_^n x_j^2 — \lambda & \sum_^n x_jy_j & \sum_^n x_j \\ \sum_^n x_jy_j & \sum_^n y_j^2 — \lambda & \sum_^n y_j \\ \sum_^n x_j & \sum_^n y_j & n \end \right|=0 \, . $$ Коэффициенты $ a, b $ и $ c $ уравнения искомой прямой определяются из системы уравнений $$ \left( \begin \sum_^n x_j^2 — \lambda_ <\min>& \sum_^n x_jy_j & \sum_^n x_j \\ \sum_^n x_jy_j & \sum_^n y_j^2 — \lambda_ <\min>& \sum_^n y_j \end \right) \left( \begin a \\ b \\ c \end \right)= \left( \begin 0 \\ 0 \\ 0 \end \right) \, . $$

Пример [Пирсон]. Построить прямую из теоремы $ 3 $ для точек

$$ \begin x & 0.0 & 0.9 & 1.8 & 2.6 & 3.3 & 4.4 & 5.2 & 6.1 & 6.5 & 7.4 \\ \hline y & 5.9 & 5.4 & 4.4 & 4.6 & 3.5 & 3.7 & 2.8 & 2.8 & 2.4 & 1.5 \end $$

Решение. Здесь $ n=10 $ и $$ \sum_^ x_j= 38.2,\ \sum_^ y_j= 37.0 \quad \Rightarrow \quad \overline x =3.82, \overline y = 3.70 \, . $$ Далее, $$ \sum_^ x_j^2= 202.32,\ \sum_^ y_j^2= 154.12, \ \sum_^ x_jy_j =110.91 \, . $$ Вычисляем определитель $$ \left| \begin 202.32 -\lambda & 110.91 & 202.32 \\ 110.91 & 154.12 — \lambda & 37.0 \\ 202.32 & 37.0 & 10 \end \right|\equiv 10\lambda^2-736.16\lambda+451.5422 \, . $$ Корни полинома: $$ \lambda_1 \approx 0.618573,\ \lambda_2 \approx 72.997427 \, . $$ Следовательно, минимальное значение суммы квадратов расстояний от заданных точек до наилучшей прямой равна $ \lambda_1 $ Подставляем это значение в систему линейных уравнений относительно $ a, b, c $. Она имеет бесконечное множество решений, из которого мы выберем, например, $$ a\approx -0.094322,\ b \approx -0.172889,\ c=1 \, . $$ Уравнение прямой совпадает с решением Пирсона, но вот значение для минимума у него приведено другое: $ 0.2484 $. Прямая проверка подтверждает значение $ \lambda_1 $. ♦

Пример [Гальтон]. В каждой семье, имеющей взрослых детей, замерим рост родителей и рост взрослых детей. Именно это проделал Ф.Гальтон в последней четверти XIX века. Результаты его измерений по $ 205 $ семьям и $ n= 898 $ детям можно найти ☞ ЗДЕСЬ. Целью исследования было установление зависимости между усредненным ростом обоих родителей и ростом их детей. Если изобразить полученные данные в виде точек $ \_^n $ плоскости, то получим следующую картину:

Здесь размеры указаны в дюймах ($ 1 $ дюйм $ \approx 2.54 $ см), а усредненный рост родителей вычисляется по формуле $$ \frac(\mbox + 1.08 \times \mbox) \, . $$ Координаты центроида: $ \overline\approx 69.22, \overline\approx 66.76 $. Уравнение $$ 898\lambda^2-13002110\lambda+27405100000=0 $$ имеет два положительных корня $$ \lambda_\approx 2560.575, \quad \lambda_\approx 11918.389 \, . $$ Уравнение прямой, обеспечивающей минимум сумме квадратов расстояний от экспериментальных точек: $$ y \approx 4.71\, x-259.40 \, , $$

В случае, когда координаты точек $ \<(x_,y_j)\>_^n $ центрированы к началу координат, т.е.

$$ \sum_^n x_j = \sum_^n y_j =0 \, ,$$ то обе прямые из теоремы $ 3 $ проходят через начало координат ($ c=0 $), числа $ \lambda_ $ и $ \lambda_ $ являются собственными числами (т.е. корнями характеристического полинома) матрицы $$ G=\left( \begin X^ X & X^ Y \\ X^ Y & Y^ Y \end \right) \quad npu \quad X^=[x_1,x_2,\dots,x_n], \ Y^=[y_1,y_2,\dots,y_n] \, . $$ Собственный вектор этой матрицы, соответствующий числу $ \lambda_ $ будет направляющим вектором прямой, обеспечивающей минимум сумме квадратов.

По аналогии хочется утверждать, что собственный вектор, соответствующий числу $ \lambda_ $ будет направляющим вектором прямой, обеспечивающей максимум сумме квадратов. Однако, это утверждение неверно. Оно будет верно, если мы добавим условие «…среди всех прямых, проходящих через начало координат».

Матрица $ G $, составленная для произвольных векторов $ \ \subset \mathbb R^n $, известна как матрица Грама этих векторов. Но более интересна другая ее интерпретация — статистическая. Матрица $$ \frac \left( \begin X^ X & X^ Y \\ X^ Y & Y^ Y \end \right) \, , $$ составленная для векторов $ \ \subset \mathbb R^n $, каждый из которых имеет нулевое среднее значение, называется ковариационной матрицей этих векторов.

Метод наименьших квадратов

Предположим снова, что данные исходной таблицы не являются достоверными: значения обеих переменных подвержены воздействию случайных погрешностей, но — в отличие от предыдущего пункта — эти погрешности имеют разные порядки. Например, будем считать, что в заданной таблице значений $$ \begin x & x_1 & x_2 & \dots & x_m \\ \hline y & y_1 & y_2 &\dots & y_m \end $$ величины $ x_<> $ нам известны практически без искажений (т.е. на входе процесса мы имеем абсолютно достоверные данные), а вот измерения величины $ y_<> $ подвержены случайным погрешностям.

Задача. Построить полином $ f_<>(x) $ такой, чтобы величина $$ \sum_^m (f(x_j)-y_j)^2 $$ стала минимальной.

В случае $ \deg f_<> =m-1 $ мы возвращаемся к задаче интерполяции в ее классической постановке. Практический интерес, однако, представляет случай $ \deg f_<> < m-1 $: экспериментальных данных больше (обычно — много больше) чем значений параметров (коэффициентов полинома $ f_<>$), которые требуется определить.

Так, в случае $ \deg f_<>=1 $ речь идет о построении прямой $ y=ax+b $ на плоскости $ (x,y) $, обеспечивающей $$ \min (\varepsilon_1^2+\varepsilon_2^2+\dots + \varepsilon_m^2) \, . $$ Здесь $ \varepsilon_j = |a\,x_j+b-y_j| $.

Обратите внимание на отличие этого рисунка от его аналога из предыдущего ПУНКТА.

Теорема. Если $ m\ge n_<> $ и $ x_,\dots,x_m $ все различны, то существует единственный набор коэффициентов $ a_,\dots,a_ $, обеспечивающий минимальное значение для

$$ \sum_^m (a_0+a_1x_j+\dots+a_x_j^ -y_j)^2 \ . $$ Этот набор определяется как решение системы нормальных уравнений $$ \underbrace< \left(\begin s_0 &s_1&s_2&\ldots&s_& s_\\ s_1 &s_2&s_3&\ldots&s_& s_\\ s_2 &s_3&s_4&\ldots&s_& s_\\ \vdots& & & && \vdots\\ s_ &s_n&s_&\ldots &s_&s_ \end\right)>_> \left(\begin a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_ \end \right)= \left(\begin t_0 \\ t_1 \\ t_2 \\ \vdots \\ t_ \end \right) $$ при $ s_ = x_1^k+\dots+x_m^k, \ t_ = x_1^ky_1+\dots+x_m^k y_m $.

Если одно из значений $ x_ $ равно $ 0_<> $ , то полагаем $ 0^ = 1 $, так что $ s_=m $ при любых $ \,\dots, x_m\> \subset \mathbb R $.

Доказательство ☞ ЗДЕСЬ.

Показать, что линейный полином $ y=a_+a_1x $, построенный по методу наименьших квадратов, определяет на плоскости $ (x_<>,y) $ прямую, проходящую через центроид

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

$$ \begin x & 0.5 & 1 & 1.5 & 2 & 2.5 & 3 \\ \hline y & 0.35 & 0.80 & 1.70 & 1.85 & 3.51 & 1.02 \end $$

Решение. Имеем $$ s_0=6, s_1=0.5 + 1 + 1.5 + 2 + 2.5 + 3=10.5, $$ $$ s_2=0.5^2 + 1^2 + 1.5^2 + 2^2 + 2.5^2 + 3^2=22.75, $$ $$t_0=0.35 + 0.80 + 1.70 + 1.85 + 3.51 + 1.02=9.23, $$ $$ t_1 =0.5\cdot 0.35 + 1 \cdot 0.80 + 1.5 \cdot 1.70 + 2 \cdot 1.85 + $$ $$ +2.5 \cdot 3.51 + 3 \cdot 1.02=19.06 . $$ Решаем систему нормальных уравнений $$ \left( \begin 6 & 10.5 \\ 10.5 & 22.75 \end \right) \left( \begin a_0 \\ a_1 \end \right)= \left( \begin 9.23 \\ 19.06 \end \right), $$ получаем уравнение прямой в виде $$ y= 0.375 + 0.664\, x\ .$$ На рисунке она изображена в цвете охры (для сравнения оставлена также зеленая прямая, полученная по методу предыдущего пункта). ♦

Дальнейшее развитие идеологии МНК ☞ ЗДЕСЬ.

Многомерная интерполяция

По аналогии с одномерным случаем можно поставить задачу интерполяции в многомерном пространстве. Пусть заданы точки (узлы интерполяции) $$ (x_,\dots,x_) \in \mathbb C^ \ npu \ j\in \ $$ и заданы значения $ \_^N $ функции в этих точках. Требуется построить полином от $ \ell $ переменных, принимающий заданные значения в узлах интерполяции.

Сложности: парадокс Крамера

Задача. Построить полином $ f_<>(x,y) $ с комплексными коэффициентами по его значениям на конечном наборе точек: $$ f(x_1,y_1)=z_1,\dots, f(x_N,y_N)=z_N \ .$$

Задача аналогична одномерной, однако ее решение в двумерном случае имеет одну особенность. Коэффициенты полинома $ f_<>(x) \in \mathbb C[x] $ степени не выше, чем $ n_<>-1 $, однозначно восстанавливаются по значениям этого полинома в произвольном наборе из $ n_<> $ различных точек $ x_,\dots,x_ $. Полином $ f_<>(x,y) \in \mathbb C[x,y] $ третьей степени определяется своими 10-ю коэффициентами: $$ f(x,y)\equiv a_x^3+a_x^2y+a_xy^2+a_y^3+a_x^2+a_xy+a_y^2+a_x+a_y+a_ \ . $$ По аналогии с одномерным случаем, можно было бы ожидать, что задание этого полинома его значениями в 9-ти произвольных узлах всегда позволит установить его коэффициенты, причем бесконечным числом способов. Эти ожидания, однако, опровергаются примером.

Пример [7]. Построить полином $ f_<>(x,y) $ третьей степени такой, что

Решение. В каноническом представлении полинома $ f_<>(x,y) $ имеется 10 коэффициентов, которые мы считаем неопределенными и ищем из заданных условий. Разрешая получающуюся систему из 9-ти линейных уравнений относительно этих коэффициентов, мы приходим к ответу: матрица системы имеет ранг 7, т.е., согласно теореме Кронекера-Капелли, сама система будет совместной только при дополнительном условии «связи» величин $ z_,\dots,z_9 $. Именно, последние должны удовлетворять определенному линейному соотношению $$ p_2z_2+\dots+p_9z_9=0 $$ при целых $ p_,\dots,p_9 $ (мы не указываем эти числа ввиду их громоздкости). Таким образом, в общем случае поставленная задача неразрешима (например, она не имеет решения при $ z_=\dots=z_9=1 $).

При $ z_,\dots,z_9 $, удовлетворяющих упомянутому уравнению «связи», поставленная задача имеет решение. Однако, оно неединственно в силу все той же теоремы Кронекера-Капелли, поскольку число совместных линейных уравнений меньше числа определяемых ими коэффициентов.

Как удалось подобрать узлы настолько «неудачные» для задачи интерполяции? Крамер предложил выбирать узлы как точки пересечения двух кривых третьего порядка (кубик). Согласно теореме Безу, две такие кривые пересекаются как раз в 9 точках. Например, если взять $$ \,x^ \,x^y \,xy^\,y^+ \,x^ \,xy \,y^ \,x \,y=0 $$ и $$ \,x^ \,x^y+ \,xy^\,y^ \,x^\,xy \,y^\,x \,y=0\ . $$ то все точки пересечения этих кривых будут иметь рациональные координаты — как раз те, что приведены выше. Система из 9-ти линейных уравнений для определения 10-ти коэффициентов интерполяционного полинома при $ z_=\dots=z_9=0 $ становится однородной. Поскольку она имеет 2 линейно независимых набора решений (соответствующие коэффициентам двух рассмотренных кубик), то ранг ее матрицы должен быть не выше 8. Если мы «сдвинем» хоть одно из значений $ z_ $ из нуля, то почти наверняка соответствующая неоднородная система станет несовместной, а задача интерполяции — неразрешимой среди полиномов третьей степени. ♦

Прямоугольная сетка

В достаточно важном — с практической точки зрения — частном случае задачу интерполяции можно решить. Пусть узлы интерполяции расположены в точках $$ (x_j,y_k) \ npu \ j \in \, k\in \ \ , $$ и в них известны значения функции $ z_ $. Требуется построить полином $ f_<>(x,y) $ по этим значениям. Такая задача возникает, например, при восстановлении рельефа морского дна с помощью промеров глубин эхолотом.

Теорема. Если все числа $ \ < x_\>_^n $ различны и все числа $ \\>_^m $ различны, то существует единственный полином $ f_<>(x,y) $, степень которого по переменной $ x_<> $ не превосходит $ n_<>-1 $, а по переменной $ y_<> $ — не превосходит $ m_<>-1 $, принимающий значения по заданной таблице.

Доказательство. Построить выражение для полинома можно обобщением формы Лагранжа. Обозначим $$ W(x) = (x-x_1)\times \dots \times (x-x_n) \ ,\quad \mathfrak W(y) = (y-y_1)\times \dots \times (y-y_m) $$ и $$ W_j(x) = \frac =(x-x_1)\times \dots \times (x-x_) (x-x_)\times \dots \times (x-x_n) \quad npu \quad j \in \, $$ и $$ \mathfrak W_k(y) = \frac<\mathfrak W(y)> =(y-y_1)\times \dots \times (y-y_) (y-y_)\times \dots \times (y-y_m) \quad npu \quad k \in \ \, . $$ Тогда $$ f(x,y)=\sum_^n \sum_^m z_ \frac \frac<\mathfrak W_k(y)> <\mathfrak W_k(y_k)>\, . $$ Легко проверить, что этот полином удовлетворяет условиям теоремы: $$ f(x_j,y_k)=z_, \ \deg f_x \le n-1, \ \deg f_y \le m-1 \ . $$

Задачи

Источники

[1]. Baron de Prony Gaspard-Clair-François-Marie Riche, Essai expérimental et analytique sur les lois de la Dilatabilité des fluides élastiques et sur celles de la Force expansive de la vapeur de l’eau et de la vapeur de l’alkool, à différentes températures, J. de l’École Polytechnique, 1, pp. 24–76, (1795). 109

[2]. Henrici P. Applied and Computational Complex Analysis. V. 1. 1974. NY. Wiley

[4]. Hilbert D. Ein Beitrag zur Theorie des Legendreschen Polynoms. Acta Math. Bd.18, 1894, S.155-160

[5]. Форсайт Дж., Молер К. Численное решение систем линейных алгебраических уравнений. М. Мир. 1969

[6]. Линник Ю.В. Метод наименьших квадратов и основы математико-статистической теории обработки наблюдений. М.ГИФМЛ. 1958

[7]. Калинина Е.А., Утешев А.Ю. Теория исключения. Учеб. пособие. СПб.: НИИ Химии СПбГУ, 2002. 72 с.

[8]. Утешев А.Ю., Тамасян Г.Ш. К задаче полиномиального интерполирования с кратными узлами. Вестник СПбГУ. Серия 10. 2010. Вып. 3, С. 76-85. Текст (pdf) ☞ ЗДЕСЬ

[9]. Cauchy A.-L. Cours d’Analyse de l’École Royale Polytechnique: Part I: Analyse Algébrique. Paris, France: L’Imprimerie Royale, 1821, pt. 1.

[10]. Jacobi C.G.J. Űber die Darstellung einer Reihe gegebner Werthe durch eine gebrochne rationale Function. J.reine angew. Math. 1846. Bd. 30, S. 127-156

[11]. Утешев А.Ю., Боровой И.И. Решение задачи рациональной интерполяции с использованием ганкелевых полиномов. Вестник СПбГУ. Серия 10. 2016. Вып. 4, С. 31-43. Текст ☞ ЗДЕСЬ (pdf).

[12]. Pearson K. On lines and planes of closest fit to systems of points in space. Phil. Mag. 1901. V.2, pp. 559-572

Интерполяционный многочлен в форме Ньютона

Выборка экспериментальных данных представляет собой массив данных, который характеризует процесс изменения измеряемого сигнала в течение заданного времени (либо относительно другой переменной). Для выполнения теоретического анализа измеряемого сигнала необходимо найти аппроксимирующую функцию, которая свяжет дискретный набор экспериментальных данных с непрерывной функцией — интерполяционным полиномом n-степени. Данный интерполяционный полином n-степени может быть записан, например, в форме Ньютона (один из способов представления).

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

1. Интерполяционная формула Ньютона для неравноотстоящих значений аргумента

В общем виде интерполяционный многочлен в форме Ньютона записывается в следующем виде:

где n – вещественное число, которое указывает степень полинома;

– переменная, которая представляет собой разделенную разность k-го порядка, которая вычисляется по следующей формуле:

Разделённая разность является симметричной функцией своих аргументов, то есть при любой их перестановке её значение не меняется. Следует отметить, что для разделённой разности k-го порядка справедлива следующая формула:

В качестве примера, рассмотрим построение полинома в форме Ньютона по представленной выборке данных, которая состоит из трех заданных точек . Интерполяционный многочлен в форме Ньютона, который проходит через три заданных точки, будет записываться в следующем виде:

• Разделенная разность 1-го порядка определяется следующим выражением

Следует отметить, что данное выражение может быть переписано в другом виде:

• Разделенная разность 2-го порядка определяется следующим выражением

Следует отметить, что данное выражение может быть переписано в другом виде:

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

Следует отметить, что полином в форме Ньютона может быть представлен в более компактном виде (по схеме Горнера), которая получается путем последовательного вынесения за скобки множителей

2. Интерполяционная формула Ньютона для равноотстоящих значений аргумента

В случае если значения функции заданы для равноотстоящих значений аргумента, которые имеют постоянный шаг измерений , то используют другую форму записи интерполяционного многочлена по формуле Ньютона.

• Для интерполирования функции в конце рассматриваемого интервала (интерполирование назад и экстраполирование вперед) используют интерполяционный полином в форме Ньютона в следующей записи:

где конечные разности k-порядка определяются по следующему выражению

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

• Для интерполирования функции в начале рассматриваемого интервала (интерполирование вперед и экстраполирование назад) используют интерполяционный полином в форме Ньютона в следующей записи:

где конечные разности k-порядка определяются по следующему выражению

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

3. Погрешность интерполяционного полинома в форме Ньютона

Рассмотрим функцию f(x), которая непрерывна и дифференцируема на рассматриваемом отрезке [a, b]. Интерполяционный полином P(x) в форме Ньютона принимает в точках заданные значения функции . В остальных точках интерполяционный полином P(x) отличается от значения функции f(x) на величину остаточного члена, который определяет абсолютную погрешность интерполяционной формулы Ньютона:

Абсолютную погрешность интерполяционной формулы Ньютона определяют следующим образом:

Переменная представляет собой верхнюю границу значения модуля (n +1)-й производной функции f(x) на заданном интервале [a, b]

В случае равноотстоящих узлов абсолютная погрешность интерполяционной формулы Ньютона определяют следующим образом:

Выражение записано с учетом следующей формулы:

Выбор узлов интерполяции

С помощью корректного выбора узлов можно минимизировать значение в оценке погрешности, тем самым повысить точность интерполяции. Данная задача может быть решена с помощью многочлена Чебышева:

В качестве узлов следует взять корни этого многочлена, то есть точки:

4. Методика вычисления полинома в форме Ньютона (прямой способ)

Алгоритм вычисления полинома в форме Ньютона позволяет разделить задачи определения коэффициентов и вычисления значений полинома при различных значениях аргумента:

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

2. Выполняется вычисление разделенных разностей n-порядка, которые будет использоваться для построения полинома в форме Ньютона.

3. Выполняется вычисление полинома n-степени в форме Ньютона по следующей формуле:

Алгоритм вычисления полинома в форме Ньютона представлен на рисунке 1.

Рис.1. Методика вычисления полинома в форме Ньютона

Следует отметить, что разделённые разности k-го порядка в соответствии с представленной методикой перезаписывается в вектор столбец функции , а результирующая разделенная разность всегда находится в первой ячейке функций . Рассмотрим, каким образом будет изменяться вектор столбец функции при выполнении расчета по представленной методике.

В качестве примера рассмотрим следующую практическую задачу. В рамках задачи известен набор шести значений, которые получены методом случайной выборки для различных моментов времени. Следует отметить, что данная выборка значений описывает функция на интервале [0, 10]. Необходимо построить многочлен в форме Ньютона для представленного набора значений. С помощью интерполяционной формулы вычислить приближенное значение функции в точке , а также определить оценку погрешности результата вычислений.

Многочлен в форме Ньютона, который строится на основании шести значений, представляет собой полином 5 степени. Результат построения полинома в форме Ньютона показан в графическом виде.

Рис.2. Исходная функция и полином в форме Ньютона, построенный по шести заданным точкам

С помощью найденного полинома можно определить значение функции в любой точке заданного интервала. Определение промежуточных значений величины по имеющемуся дискретному набору известных значений называется «интерполяцией». В соответствии с условиями задачи полином в форме Ньютона в точке x=9,5 принимает следующее значение: L(9,5)= – 4,121. Из графика видно, что полученное значение не совпадает cо значением функции f(x) на величину абсолютной погрешности интерполяционной формулы Ньютона.

Интерполяционный полином в форме Ньютона часто оказывается удобным для проведения различных теоретических исследований в области вычислительной математики. Так, например, полином в форме Ньютона используются для интерполяции, а также для численного интегрирования таблично-заданной функцией.