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

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

Определение

Интерполяционный многочлен Лагранжа для четырёх точек (-9,5), (-4,2), (-1,-2) и (7,9), а также полиномы $y_i{l}_i(x)$, каждый из которых проходит через одну из выделенных точек, и принимает нулевое значение в остальных.

Пусть задана $n+1$ пара чисел $(x_0,y_0),(x_1,y_1),\dots ,(x_n,y_n),$ где все $x_j$ различны. Требуется построить многочлен $L(x)$ степени не более $n$, для которого $L(x_j)=y_j$.

Общий случай

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

$L(x)=\sum _{i=0}^n{y}_i{l}_i(x),$

где базисные полиномы $l_i$ определяются по формуле

$l_i(x)=\prod _{j=0,j\ne i}^n\frac{x-x_j}{x_i-x_j}=\frac{x-x_0}{x_i-x_0}\cdots \frac{x-x_{i-1}}{x_i-x_{i-1}}\cdot \frac{x-x_{i+1}}{x_i-x_{i+1}}\cdots \frac{x-x_n}{x_i-x_n}$

Для любого $i=0,\dots ,n$ многочлен $l_i$ имеет степень $n$ и

$l_i(x_j)=\{\begin{array}{cc}0,& j\ne i,\\ 1,& j=i.\end{array}$

Отсюда следует, что $L(x)$, являющийся линейной комбинацией многочленов $l_i(x)$, имеет степень не больше $n$ и $L(x_i)=y_i$.

Случай равноотстоящих узлов интерполяции

Пусть узлы интерполяции $x_j$ являются равноотстоящими, то есть выражаются через начальную точку $x_0$ и некоторую фиксированную положительную величину $h$ следующим образом:

$x_j=x_0+jh.$

Отсюда следует, что

$x_j-x_i=(j-i)h.$

Подставляя эти выражения в формулу для базисного полинома и вынося $h$ за знаки произведения в числителе и знаменателе, получим

$l_j(x)=\prod _{i=0,i\ne j}^n\frac{(x-x_i)}{(x_j-x_i)}=\frac{\prod _{i=0,i\ne j}^n(x-x_0-ih)}{h^n\prod _{i=0,i\ne j}^n(j-i)}$

Теперь можно ввести замену переменной

$y=\frac{x-x_0}{h}$

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

$l_j(x)=l_j(x_0+yh)=(-1{)}^{n-j}{C}_n^j\frac{y(y-1)\dots (y-n)}{(y-j)n!}.$

Данные величины называются коэффициентами Лагранжа. Они не зависят ни от $y_0,\dots ,y_n$, ни от $h$ и потому могут быть вычислены заранее и записаны в виде таблиц. Недостатком данного подхода является факториальная сложность числителя и знаменателя, что требует использования длинной арифметики.

Остаточный член

Если считать числа $y_0,\dots ,y_n$ значениями некоторой функции $f$ в узлах $x_0,\dots ,x_n$, то ошибка интерполирования функции $f$ многочленом $L$ равна

$f(x)-L(x)=\frac{f^{(n+1)}(\xi )}{(n+1)!}(x-x_0)\dots (x-x_n),$

где $\xi$ — некоторая средняя точка между наименьшим и наибольшим из чисел $x_0,\dots ,x_n$. Полагая $M_{n+1}=\underset{x\in [x_0,x_n]}{\sup }|f^{(n+1)}(x)|$, можно записать

$|f(x)-L(x)|⩽\frac{M_{n+1}}{(n+1)!}|(x-x_0)\dots (x-x_n)|.$

Единственность

Шаблон:Основная статья Существует единственный многочлен степени не превосходящей $n$, принимающий заданные значения в $n+1$ заданной точке. Шаблон:Доказательство Это утверждение является обобщением того факта, что через любые две точки проходит единственная прямая.

С точки зрения линейной алгебры

На единственность интерполяционного многочлена можно также взглянуть с точки зрения СЛАУ. Рассмотрим систему уравнений $P(x_0)=y_0,P(x_1)=y_1,\dots ,P(x_n)=y_n$. В явном виде она записывается как

$\{\begin{array}{c}{a}_0+a_1{x}_0+a_2{x}_0^2+\dots +a_n{x}_0^n=y_0,\\ a_0+a_1{x}_1+a_2{x}_1^2+\dots +a_n{x}_1^n=y_1,\\ \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots ,\\ a_0+a_1{x}_n+a_2{x}_n^2+\dots +a_n{x}_n^n=y_n\end{array}$

Её можно переписать в виде системы уравнений $Xa=y$ с неизвестным вектором $a=(a_0,a_1,\dots ,a_n)$:

$\left[\begin{array}{ccccc}1& x_0& x_0^2& \dots & x_0^n\\ 1& x_1& x_1^2& \dots & x_1^n\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_n& x_n^2& \dots & x_n^n\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_n\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_n\end{array}\right]$

Матрица $X$ в такой системе является матрицей Вандермонда и её определитель равен $\prod _{i<j}(x_i-x_j)$. Соответственно, если все точки $x_0,x_1,\dots ,x_n$ различны, то матрица невырождена и система обладает единственным решением.

С точки зрения китайской теоремы об остатках

По теореме Безу остаток от деления $P(x)$ на $(x-a)$ равен $P(a)$. Таким образом, всю систему можно воспринимать в виде системы сравнений:

$\{\begin{array}{c}P(x)\equiv y_0(\mathrm{mod}x-x_0),\\ P(x)\equiv y_1(\mathrm{mod}x-x_1),\\ \dots ,\\ P(x)\equiv y_n(\mathrm{mod}x-x_n),\end{array}$

По китайской теореме об остатках у такой системы есть единственное решение по модулю $(x-x_0)(x-x_1)\dots (x-x_n)$, то есть, заданная система однозначно определяет многочлен степени не выше $n$. Такое представление многочлена в виде наборов остатков по модулям мономов аналогично представлению числа в виде остатков от деления на простые модули в системе остаточных классов. При этом явная формула для многочлена Лагранжа также может быть получена в соответствии с формулами китайской теоремы: $P(x)=\sum _{i=0}^n{y}_i{M}_i{M}_i^{-1}$, где $M_i=\prod _{j\ne i}(x-x_j)$ и $M_i^{-1}=\prod _{j\ne i}(x-x_j{)}^{-1}(\mathrm{mod}x-x_i)=\prod _{j\ne i}(x_i-x_j{)}^{-1}$.

Пример

Приближение функции $f(x)=\text{tg}(x)$ (синяя линия) многочленом $L(x)=4,834848x^3-1,477474x$ (зелёная линия).

Найдем формулу интерполяции для $f(x)=\mathrm{t}\mathrm{g}\mathrm{\backslash nolimits}(x)$ имеющей следующие значения:

$\begin{array}{cccccccc}{x}_0& =-1.5& & & & & f(x_0)& =-14,1014\\ x_1& =-0.75& & & & & f(x_1)& =-0,931596\\ x_2& =0& & & & & f(x_2)& =0\\ x_3& =0.75& & & & & f(x_3)& =0,931596\\ x_4& =1.5& & & & & f(x_4)& =14,1014.\end{array}$

$l_0(x)=\frac{x-x_1}{x_0-x_1}\cdot \frac{x-x_2}{x_0-x_2}\cdot \frac{x-x_3}{x_0-x_3}\cdot \frac{x-x_4}{x_0-x_4}=\frac{1}{243}x(2x-3)(4x-3)(4x+3)$

$l_1(x)=\frac{x-x_0}{x_1-x_0}\cdot \frac{x-x_2}{x_1-x_2}\cdot \frac{x-x_3}{x_1-x_3}\cdot \frac{x-x_4}{x_1-x_4}=-\frac{8}{243}x(2x-3)(2x+3)(4x-3)$

$l_2(x)=\frac{x-x_0}{x_2-x_0}\cdot \frac{x-x_1}{x_2-x_1}\cdot \frac{x-x_3}{x_2-x_3}\cdot \frac{x-x_4}{x_2-x_4}=\frac{3}{243}(2x+3)(4x+3)(4x-3)(2x-3)$

$l_3(x)=\frac{x-x_0}{x_3-x_0}\cdot \frac{x-x_1}{x_3-x_1}\cdot \frac{x-x_2}{x_3-x_2}\cdot \frac{x-x_4}{x_3-x_4}=-\frac{8}{243}x(2x-3)(2x+3)(4x+3)$

$l_4(x)=\frac{x-x_0}{x_4-x_0}\cdot \frac{x-x_1}{x_4-x_1}\cdot \frac{x-x_2}{x_4-x_2}\cdot \frac{x-x_3}{x_4-x_3}=\frac{1}{243}x(2x+3)(4x-3)(4x+3).$

Получим

$\begin{array}{cc}L(x)& =\frac{1}{243}(f(x_0)x(2x-3)(4x-3)(4x+3)\\ & -8f(x_1)x(2x-3)(2x+3)(4x-3)\\ & +3f(x_2)(2x+3)(4x+3)(4x-3)(2x-3)\\ & -8f(x_3)x(2x-3)(2x+3)(4x+3)\\ & +f(x_4)x(2x+3)(4x-3)(4x+3))\\ & =4,834848x^3-1,477474x.\end{array}$

Реализация общего случая на языке программирования Python

import numpy as np

# данные для примера
xi = np.array([-1.5, -0.75, 0, 0.75, 1.5])
yi = np.array([-14.1014, -0.931596, 0, 0.931596, 14.1014])


def get_coefficients(_pl: int, _xi: np.ndarray):
    '''
    Определение коэффициентов для множителей базисных полиномов l_i
    :param _pl: индекс базисного полинома
    :param _xi: массив значений x
    :return:
    '''
    n = int(_xi.shape[0])
    coefficients = np.empty((n, 2), dtype=float)
    for i in range(n):
        if i == _pl:
            coefficients[i][0] = float('inf')
            coefficients[i][1] = float('inf')
        else:
            coefficients[i][0] = 1 / (_xi[_pl] - _xi[i])
            coefficients[i][1] = -_xi[i] / (_xi[_pl] - _xi[i])
    filtered_coefficients = np.empty((n - 1, 2), dtype=float)
    j = 0
    for i in range(n):
        if coefficients[i][0] != float('inf'):
            # изменение последовательности, степень увеличивается
            filtered_coefficients[j][0] = coefficients[i][1]
            filtered_coefficients[j][1] = coefficients[i][0]
            j += 1
    return filtered_coefficients


def get_polynomial_l(_xi: np.ndarray):
    '''
    Определение базисных полиномов
    :param _xi: массив значений x
    :return:
    '''
    n = int(_xi.shape[0])
    pli = np.zeros((n, n), dtype=float)
    for pl in range(n):
        coefficients = get_coefficients(pl, _xi)
        for i in range(1, n - 1):  # проходим по массиву coefficients
            if i == 1:  # на второй итерации занимаются 0-2 степени
                pli[pl][0] = coefficients[i - 1][0] * coefficients[i][0]
                pli[pl][1] = coefficients[i - 1][1] * coefficients[i][0] + coefficients[i][1] * coefficients[i - 1][0]
                pli[pl][2] = coefficients[i - 1][1] * coefficients[i][1]
            else:
                clone_pli = np.zeros(n, dtype=float)
                for val in range(n):
                    clone_pli[val] = pli[pl][val]
                zeros_pli = np.zeros(n, dtype=float)
                for j in range(n-1):  # проходим по строке pl массива pli
                    product_1 = clone_pli[j] * coefficients[i][0]
                    product_2 = clone_pli[j] * coefficients[i][1]
                    zeros_pli[j] += product_1
                    zeros_pli[j+1] += product_2
                for val in range(n):
                    pli[pl][val] = zeros_pli[val]
    return pli

def get_polynomial(_xi: np.ndarray, _yi: np.ndarray):
    '''
    Определение интерполяционного многочлена Лагранжа
    :param _xi: массив значений x
    :param _yi: массив значений y
    :return:
    '''
    n = int(_xi.shape[0])
    polynomial_l = get_polynomial_l(_xi)
    for i in range(n):
        for j in range(n):
            polynomial_l[i][j] *= _yi[i]
    L = np.zeros(n, dtype=float)
    for i in range(n):
        for j in range(n):
            L[i] += polynomial_l[j][i]
    return L

# результат в виде массива коэффициентов многочлена при x в порядке увеличения степени
# [ 0.         -1.47747378  0.          4.8348476   0.        ]
# т.е. результирующий многочлен имеет вид: y(x) = -1.47747378*x + 4.8348476*x^3

Применения

Многочлены Лагранжа степеней от нулевой до пятой для функции $\cos (5\pi x)$

Численное интегрирование

Пусть для функции $f(x)$ известны значения $y_i=f(x_i)$ в некоторых точках. Тогда можно интерполировать эту функцию методом Лагранжа:

$f(x)\approx \sum _{i=0}^{n}f(x_i)l_i(x).$

Полученное выражение можно использовать для приближённого вычисления определённого интеграла от функции $f$:

$\underset{a}{\overset{b}{\int }}f(x)dx\approx \sum _{i=0}^{n}f(x_i)\underset{a}{\overset{b}{\int }}{l}_i(x)dx$

Значения интегралов от $l_i$ не зависят от $f(x)$ и их можно вычислить заранее с использованием последовательности $x_j$.

Литература

• Березин, И. С., Жидков Н. П. Методы вычислений. Том I. — 2-е изд., стереотипное — М.: Физматлит. 1962.

Ссылки

• Шаблон:КнигаШаблон:Недоступная ссылка
• А. Г. Хованский. Полиномы Лагранжа и их применения. Видео-лекция. VI Летняя школа «Современная математика», Дубна, 2006.

См. также

• Схема Эйткена
• Интерполяционные формулы Ньютона
• Интерполяция с кратными узлами
• Схема разделения секрета Шамира
• Комбинаторная теорема о нулях
• Интерполяция алгебраическими многочленами