Русская Википедия:Интерполяционный многочлен Лагранжа

Материал из Онлайн справочника
Перейти к навигацииПерейти к поиску

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

Определение

Файл:Lagrange polynomial.svg
Интерполяционный многочлен Лагранжа для четырёх точек (-9,5), (-4,2), (-1,-2) и (7,9), а также полиномы <math>y_i l_i(x)</math>, каждый из которых проходит через одну из выделенных точек, и принимает нулевое значение в остальных.

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

Общий случай

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

<math>L(x) = \sum_{i=0}^n y_i l_i(x),</math>

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

<math>l_i(x)=\prod_{j=0, j\neq 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}}</math>

Для любого <math>i = 0,\ldots,n</math> многочлен <math>l_i</math> имеет степень <math>n</math> и

<math> l_i(x_j) = \left\{ \begin{array}{rl} 0, &j \neq i, \\ 1, &j = i. \end{array} \right. </math>

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

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

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

<math> x_j = x_0 + jh.</math>

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

<math> x_j - x_i = (j - i)h.</math>

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

<math>l_j(x) = { \prod_{i=0,\,i \ne j}^n {(x - x_i) \over (x_j - x_i)}} = {\prod\limits_{i=0,\,i \ne j}^n (x - x_0 - ih) \over h^n\prod\limits_{i=0,\,i \ne j}^n (j - i)}</math>

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

<math>y = {{x - x_0} \over h}</math>

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

<math>l_j(x) = l_j(x_0 + yh) = (-1)^{n-j} C_n^j \dfrac{y (y-1) \ldots (y-n)}{(y-i) n!}. </math>

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

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

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

<math> f(x) - L(x) = \dfrac{f^{(n+1)}(\xi)}{(n+1)!} (x-x_0) \ldots (x-x_n),</math>

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

<math> |f(x) - L(x)| \leqslant \dfrac{M_{n+1}}{(n+1)!} |(x-x_0) \ldots (x-x_n)|.</math>

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

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

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

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

<math>

\begin{cases} 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{cases}</math>

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

<math>

\begin{bmatrix} 1 & x_0 & x_0^2 & \dots & x_0^n \\ 1 & x_1 & x_1^2 & \dots & x_1^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \dots & x_n^n \\ \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_n \end{bmatrix} </math>

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

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

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

<math>

\begin{cases} P(x) \equiv y_0 \pmod{x-x_0},\\ P(x) \equiv y_1 \pmod{x-x_1},\\ \dots,\\ P(x) \equiv y_n \pmod{x-x_n},\\ \end{cases} </math>

По китайской теореме об остатках у такой системы есть единственное решение по модулю <math>(x-x_0)(x-x_1)\dots(x-x_n)</math>, то есть, заданная система однозначно определяет многочлен степени не выше <math>n</math>. Такое представление многочлена в виде наборов остатков по модулям мономов аналогично представлению числа в виде остатков от деления на простые модули в системе остаточных классов. При этом явная формула для многочлена Лагранжа также может быть получена в соответствии с формулами китайской теоремы: <math display=inline> P(x) = \sum\limits_{i=0}^n y_i M_i M_i^{-1} </math>, где <math>M_i = \prod\limits_{j \neq i} (x-x_j)</math> и <math>M_i^{-1} = \prod\limits_{j \neq i}(x-x_j)^{-1}\pmod{x-x_i}= \prod\limits_{j\neq i}(x_i-x_j)^{-1}</math>.

Пример

Файл:LagrangeInterp.png
Приближение функции <math>f(x) = \mbox{tg}(x)</math> (синяя линия) многочленом <math>L(x) = 4,834848x^3 - 1,477474x</math> (зелёная линия).

Найдем формулу интерполяции для <math>f(x) = \mathop{\mathrm{tg}}\nolimits(x)</math> имеющей следующие значения:

<math>

\begin{align} 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{align} </math>

<math>l_0(x)={x - x_1 \over x_0 - x_1}\cdot{x - x_2 \over x_0 - x_2}\cdot{x - x_3 \over x_0 - x_3}\cdot{x - x_4 \over x_0 - x_4}
            ={1\over 243} x (2x-3)(4x-3)(4x+3)</math>
<math>l_1(x) = {x - x_0 \over x_1 - x_0}\cdot{x - x_2 \over x_1 - x_2}\cdot{x - x_3 \over x_1 - x_3}\cdot{x - x_4 \over x_1 - x_4}
            = {} -{8\over 243} x (2x-3)(2x+3)(4x-3)</math>
<math>l_2(x)={x - x_0 \over x_2 - x_0}\cdot{x - x_1 \over x_2 - x_1}\cdot{x - x_3 \over x_2 - x_3}\cdot{x - x_4 \over x_2 - x_4}
            ={3\over 243} (2x+3)(4x+3)(4x-3)(2x-3) </math>
<math>l_3(x)={x - x_0 \over x_3 - x_0}\cdot{x - x_1 \over x_3 - x_1}\cdot{x - x_2 \over x_3 - x_2}\cdot{x - x_4 \over x_3 - x_4}
            =-{8\over 243} x (2x-3)(2x+3)(4x+3)</math>
<math>l_4(x)={x - x_0 \over x_4 - x_0}\cdot{x - x_1 \over x_4 - x_1}\cdot{x - x_2 \over x_4 - x_2}\cdot{x - x_3 \over x_4 - x_3}
            ={1\over 243} x (2x+3)(4x-3)(4x+3).</math>

Получим

<math> \begin{align}L(x) &= {1\over 243}\Big(f(x_0)x (2x-3)(4x-3)(4x+3) \\

& {} \qquad {} - 8f(x_1)x (2x-3)(2x+3)(4x-3) \\ & {} \qquad {} + 3f(x_2)(2x+3)(4x+3)(4x-3)(2x-3) \\ & {} \qquad {} - 8f(x_3)x (2x-3)(2x+3)(4x+3) \\ & {} \qquad {} + f(x_4)x (2x+3)(4x-3)(4x+3)\Big)\\ & = 4,834848x^3 - 1,477474x. \end{align} </math>

Реализация общего случая на языке программирования 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

Применения

Файл:Lagrange polynomials of increasing degrees.gif
Многочлены Лагранжа степеней от нулевой до пятой для функции <math>\cos(5\pi x)</math>

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

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

<math>f(x) \approx \sum_{i=0}^n f(x_i) l_i(x).</math>

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

<math>\int\limits_a^b f(x)dx \approx \sum_{i=0}^n f(x_i) \int\limits_a^b l_i(x) dx</math>

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

Литература

Ссылки

См. также

Шаблон:Викиучебник

Партнерские ресурсы
Криптовалюты
Магазины
Хостинг
Разное

Источник — http://wikihandbk.com/ruwiki/index.php?title=Русская_Википедия:Интерполяционный_многочлен_Лагранжа&oldid=10004786