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

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

Определение

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

Пусть задана $n+1$ $n+1$ пара чисел $(x_{0},y_{0}),(x_{1},y_{1}),\ldots ,(x_{n},y_{n}),$ $(x_{0},y_{0}),(x_{1},y_{1}),\ldots ,(x_{n},y_{n}),$ где все $x_{j}$ $x_{j}$ различны. Требуется построить многочлен $L(x)$ $L(x)$ степени не более $n$ $n$ , для которого $L(x_{j})=y_{j}$ $L(x_{j})=y_{j}$ .

Общий случай

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

L(x)=\sum _{i=0}^{n}y_{i}l_{i}(x),

L(x)=\sum _{i=0}^{n}y_{i}l_{i}(x),

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

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}}}

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}}}

Для любого $i=0,\ldots ,n$ $i=0,\ldots ,n$ многочлен $l_{i}$ $l_{i}$ имеет степень $n$ $n$ и

l_{i}(x_{j})=\left\{{\begin{array}{rl}0,&j\neq i,\\1,&j=i.\end{array}}\right.

l_{i}(x_{j})=\left\{{\begin{array}{rl}0,&j\neq i,\\1,&j=i.\end{array}}\right.

Отсюда следует, что $L(x)$ $L(x)$ , являющийся линейной комбинацией многочленов $l_{i}(x)$ $l_{i}(x)$ , имеет степень не больше $n$ $n$ и $L(x_{i})=y_{i}$ $L(x_{i})=y_{i}$ .

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

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

x_{j}=x_{0}+jh.

x_{j}=x_{0}+jh.

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

x_{j}-x_{i}=(j-i)h.

x_{j}-x_{i}=(j-i)h.

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

l_{j}(x)={\prod _{i=0,\,i\neq j}^{n}{(x-x_{i}) \over (x_{j}-x_{i})}}={\prod \limits _{i=0,\,i\neq j}^{n}(x-x_{0}-ih) \over h^{n}\prod \limits _{i=0,\,i\neq j}^{n}(j-i)}

l_{j}(x)={\prod _{i=0,\,i\neq j}^{n}{(x-x_{i}) \over (x_{j}-x_{i})}}={\prod \limits _{i=0,\,i\neq j}^{n}(x-x_{0}-ih) \over h^{n}\prod \limits _{i=0,\,i\neq j}^{n}(j-i)}

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

y={{x-x_{0}} \over h}

y={{x-x_{0}} \over h}

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

l_{j}(x)=l_{j}(x_{0}+yh)=(-1)^{n-j}C_{n}^{j}{\dfrac {y(y-1)\ldots (y-n)}{(y-j)n!}}.

l_{j}(x)=l_{j}(x_{0}+yh)=(-1)^{n-j}C_{n}^{j}{\dfrac {y(y-1)\ldots (y-n)}{(y-j)n!}}.

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

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

Если считать числа $y_{0},\ldots ,y_{n}$ $y_{0},\ldots ,y_{n}$ значениями некоторой функции $f$ $f$ в узлах $x_{0},\ldots ,x_{n}$ $x_{0},\ldots ,x_{n}$ , то ошибка интерполирования функции $f$ $f$ многочленом $L$ $L$ равна

f(x)-L(x)={\dfrac {f^{(n+1)}(\xi)}{(n+1)!}}(x-x_{0})\ldots (x-x_{n}),

f(x)-L(x)={\dfrac {f^{(n+1)}(\xi)}{(n+1)!}}(x-x_{0})\ldots (x-x_{n}),

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

|f(x)-L(x)|\leqslant {\dfrac {M_{n+1}}{(n+1)!}}|(x-x_{0})\ldots (x-x_{n})|.

|f(x)-L(x)|\leqslant {\dfrac {M_{n+1}}{(n+1)!}}|(x-x_{0})\ldots (x-x_{n})|.

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

Существует единственный многочлен степени не превосходящей $n$ $n$ , принимающий заданные значения в $n+1$ $n+1$ заданной точке.

Доказательство

Предположим, что существуют два различных многочлена $P(x)$ $P(x)$ и $Q(x)$ $Q(x)$ степени не более $n$ $n$ , для которых верно, что для $n+1$ $n+1$ пар чисел $(x_{0},y_{0}),(x_{1},y_{1}),\ldots ,(x_{n},y_{n}),$ $(x_{0},y_{0}),(x_{1},y_{1}),\ldots ,(x_{n},y_{n}),$ где все $x_{j}$ $x_{j}$ различны, $P(x_{j})=Q(x_{j})=y_{j}.$ $P(x_{j})=Q(x_{j})=y_{j}.$ Рассмотрим многочлен $L(x)=P(x)-Q(x)$ $L(x)=P(x)-Q(x)$ . Подставляя в него $x_{j}$ $x_{j}$ ( $0\leq j\leq n$ $0\leq j\leq n$ ), получаем, что $L(x_{j})=P(x_{j})-Q(x_{j})=y_{j}-y_{j}=0$ $L(x_{j})=P(x_{j})-Q(x_{j})=y_{j}-y_{j}=0$ . Таким образом, многочлен $L(x)$ $L(x)$ имеет $n+1$ $n+1$ корней и все они различны. Следовательно $L(x)\equiv 0$ $L(x)\equiv 0$ , так как ненулевой многочлен степени не превосходящей $n$ $n$ имеет не более $n$ $n$ корней. Следовательно, $P(x)=Q(x)$ $P(x)=Q(x)$ . ■

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

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

На единственность интерполяционного многочлена можно также взглянуть с точки зрения СЛАУ . Рассмотрим систему уравнений $P(x_{0})=y_{0},P(x_{1})=y_{1},\dots ,P(x_{n})=y_{n}$ $P(x_{0})=y_{0},P(x_{1})=y_{1},\dots ,P(x_{n})=y_{n}$ . В явном виде она записывается как

{\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}}

{\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}}

Её можно переписать в виде системы уравнений $Xa=y$ $Xa=y$ с неизвестным вектором $a=(a_{0},a_{1},\dots ,a_{n})$ $a=(a_{0},a_{1},\dots ,a_{n})$ :

{\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}}

{\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}}

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

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

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

{\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}}

{\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}}

По китайской теореме об остатках у такой системы есть единственное решение по модулю $(x-x_{0})(x-x_{1})\dots (x-x_{n})$ $(x-x_{0})(x-x_{1})\dots (x-x_{n})$ , то есть, заданная система однозначно определяет многочлен степени не выше $n$ $n$ . Такое представление многочлена в виде наборов остатков по модулям мономов аналогично представлению числа в виде остатков от деления на простые модули в системе остаточных классов . При этом явная формула для многочлена Лагранжа также может быть получена в соответствии с формулами китайской теоремы : ${\textstyle P(x)=\sum \limits _{i=0}^{n}y_{i}M_{i}M_{i}^{-1}}$ ${\textstyle P(x)=\sum \limits _{i=0}^{n}y_{i}M_{i}M_{i}^{-1}}$ , где $M_{i}=\prod \limits _{j\neq i}(x-x_{j})$ $M_{i}=\prod \limits _{j\neq i}(x-x_{j})$ и $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}$ $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}$ .

Пример

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

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

{\begin{aligned}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{aligned}}

{\begin{aligned}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{aligned}}

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)

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)

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)

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)

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)

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)

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)

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)

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).

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).

Получим

{\begin{aligned}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{aligned}}

{\begin{aligned}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{aligned}}

Реализация общего случая на языке программирования 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)$ $\cos(5\pi x)$

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

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

f(x)\approx \sum _{i=0}^{n}f(x_{i})l_{i}(x).

f(x)\approx \sum _{i=0}^{n}f(x_{i})l_{i}(x).

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

\int \limits _{a}^{b}f(x)dx\approx \sum _{i=0}^{n}f(x_{i})\int \limits _{a}^{b}l_{i}(x)dx

\int \limits _{a}^{b}f(x)dx\approx \sum _{{i=0}}^{n}f(x_{i})\int \limits _{a}^{b}l_{i}(x)dx

Значения интегралов от $l_{i}$ $l_{i}$ не зависят от $f(x)$ $f(x)$ и их можно вычислить заранее с использованием последовательности $x_{j}$ $x_{j}$ .

Литература

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

Ссылки

М. А. Тынкевич. Глава 7.6.1. Интерполяционный многочлен Лагранжа // . — Кемерово, 2002. — ISBN 5-89070-042-1 . (недоступная ссылка)
А. Г. Хованский. . Видео-лекция. VI Летняя школа «Современная математика», Дубна, 2006.