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

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

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

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

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

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

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

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

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

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

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

${\displaystyle x_{j}=x_{0}+jh.}$

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

${\displaystyle x_{j}-x_{i}=(j-i)h.}$

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

${\displaystyle 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)}}$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

${\displaystyle 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)}$

${\displaystyle 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)}$

${\displaystyle 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)}$

${\displaystyle 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)}$

${\displaystyle 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).}$

Получим

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

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)foriinrange(n):ifi==_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=0foriinrange(n):ifcoefficients[i][0]!=float('inf'):# изменение последовательности, степень увеличиваетсяfiltered_coefficients[j][0]=coefficients[i][1]filtered_coefficients[j][1]=coefficients[i][0]j+=1returnfiltered_coefficientsdef get_polynomial_l(_xi:np.ndarray):    '''    Определение базисных полиномов    :param _xi: массив значений x    :return:    '''n=int(_xi.shape[0])pli=np.zeros((n,n),dtype=float)forplinrange(n):coefficients=get_coefficients(pl,_xi)foriinrange(1,n-1):# проходим по массиву coefficientsifi==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)forvalinrange(n):clone_pli[val]=pli[pl][val]zeros_pli=np.zeros(n,dtype=float)forjinrange(n-1):# проходим по строке pl массива pliproduct_1=clone_pli[j]*coefficients[i][0]product_2=clone_pli[j]*coefficients[i][1]zeros_pli[j]+=product_1zeros_pli[j+1]+=product_2forvalinrange(n):pli[pl][val]=zeros_pli[val]returnplidef 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)foriinrange(n):forjinrange(n):polynomial_l[i][j]*=_yi[i]L=np.zeros(n,dtype=float)foriinrange(n):forjinrange(n):L[i]+=polynomial_l[j][i]returnL# результат в виде массива коэффициентов многочлена при x в порядке увеличения степени# [ 0.         -1.47747378  0.          4.8348476   0.        ]# т.е. результирующий многочлен имеет вид: y(x) = -1.47747378*x + 4.8348476*x^3

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

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

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

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

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

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

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

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