Численное интегрирование: многочлен Лагранжа

Рассмотрим многочлен вида

$P_n^0(x) = \frac{x - x_1}{x_0 - x_1} \cdot \frac{x - x_2}{x_0 - x_2} \cdot \cdots \cdot \frac{x - x_n}{x_0 - x_n} = \prod_{i=1}^n \frac{x - x_i}{x_0 - x_i}.$

Очевидно, что $P_n^0(x_i) = \delta_{0i}$ . Сложив $n+1$ подобных многочленов можно получить

$L_n(x) = \sum_{i=0}^n f(x_i) \prod_{j=0}^n' \frac{x - x_j}{x_i - x_j},$

где штрих у произведения означает, что произведение берётся по всем $j \neq i$ . Этот многочлен называется многочленом Лагранжа. Нетрудно видеть, что

$L_n(x_j) = \sum_{i=0}^n f(x_i) \delta_{ij} = f(x_j),$

то есть этот многочлен совпадает с интерполяционным $P_n(x)$ . Однако, при его интегрировании обнаруживается интересная особенность: для интегрирования нужно знать лишь интегралы от полиномов вроде рассмотренного выше $P_n^0(x)$ :

$\int\limits_a^b L_n(x)\,dx = \int\limits_a^b \sum_{i=0}^n f(x_i) \prod_{j=0}^n' \frac{x - x_j}{x_i - x_j}\,dx = \sum_{i=0}^n f(x_i) \int\limits_a^b\prod_{j=0}^n' \frac{x - x_j}{x_i - x_j}\,dx = \sum_{i=0}^n c_i f(x_i),$

где коэффициенты

$c_i = \int\limits_a^b\prod_{j=0}^n' \frac{x - x_j}{x_i - x_j}\,dx$

зависят только от разбиения и не зависят от интегрируемой функции. Следовательно, они общие для всех функций и для расчёта интеграла больше не нужно решать систему.

Для большего обобщения можно рассмотреть отображение

$y : [a, b] \to [-1, 1],\ y(x) = \frac{2x - b - a}{b - a}$

Тогда

$c_i = \frac{b - a}{2}\int\limits_{-1}^{1} \prod_{j=0}^n' \frac{y - y_j}{y_i - y_j}\,dy.$

Точки $y_i$ называются узлами разбиения, а соответствующие им значения $c_i$ – весами.

Теперь обратимся к задаче поиска весов. Подынтегральное произведение является многочленом степени $n$ и может быть представлено в виде

$\prod_{j=0}^n' \frac{y - y_j}{y_i - y_j} = \frac{P_n(y)}{P_n(y_i)}.$

Для определения коэффициентов $P_n(y)$ рассмотрим вспомогательную задачу: определим коэффициенты произведения многочлена степени $n$ на двучлен $(x - a)$ . Расписывая, получим:

$P_n(x)\cdot(x-a) = x P_n(x) - a P_n(x) = \sum_{i=1}^{n+1} p_{i-1} x^i - \sum_{i=0}^n ap_i x^i = p_n x^{n + 1} + \sum_{i=1}^n (p_{i-1} - ap_i) x^i - ap_0.$

Этот процесс можно описать в программе, представив многочлен в виде списка его коэффициентов:

def polynomial_multiply(p, a):
    '''
        Вычисляет коэффициенты многочлена, который получается при умножении
        многочлена с коэффициентами p на (x-a)

        p -- список коэффициентов [a0, a1, … , an]
    '''
    n = len(p)
    result = []

    # начинаем формировать результат:
    # 1. добавляем свободный член
    result.append(-a * p[0])

    # 2. добавляем все коэффициенты, кроме самого старшего
    for i in range(1, n):
        result.append(p[i-1] - a * p[i])

    # 3. добавляем старший коэффициент
    result.append(p[-1])
    return result

# пример использования:
polynomial_multiply([1], 2)     # [-2, 1]    : 1 * (x - 2) = x - 2
polynomial_multiply([-2, 1], 2) # [4, -4, 1] : (x - 2) * (x - 2) = x^2 - 4x + 4

Используя полученное выражение для коэффициентов, можно определить коэффициенты $P_n(y)$ , последовательно домножая $1$ на $y - y_j\ (j \ne i)$ :

def polynomial_from_roots(ys):
    '''
        Вычисляет список коэффициентов многочлена, имеющего корни ys
    '''
    result = [1]
    for y in ys:
        result = polynomial_multiply(result, y)
    return result

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

def polynomial_value(p, x):
    '''
        Вычисляет значение многочлена с коэффициентами p в точке x
    '''
    return sum(a * x ** i for i, a in enumerate(p))

Для определения веса нужно проинтегрировать полученный многочлен на промежутке $[-1, 1]$ :

$\int\limits_{-1}^1 P_n(x)\,dx = \int\limits_{-1}^1\sum_{i=0}^n a_i x^i\,dx = \sum_{i=0}^n a_i \int\limits_{-1}^1 x^i\,dx = \sum_{i=0}^n a_i \frac{1 + (-1)^i}{i+1}.$

Видно, что слагаемые с нечётными $i$ обращаются в $0$ , что позволяет написать вот такой вот код:

def polynomial_integrate(p):
    '''
        Интегрирует многочлен с коэффициентами p на [-1, 1]
    '''
    return sum(2 * a / (i + 1) for i, a in enumerate(p) if i % 2 == 0)

Ну а теперь, собрав всё вместе, можно определить веса:

def weights(ys):
    '''
        Рассчитывает веса
    '''
    result = []
    for i, y in enumerate(ys):
        # сначала формируем список [yj], j ≠ i
        yj = ys[:i] + ys[i + 1:]

        # определяем коэффициенты числителя
        p = polynomial_from_roots(yj)

        # определяем знаменатель
        value = polynomial_value(p, y)

        # рассчитываем вес и добавляем его в список весов
        result.append(polynomial_integrate(p) / value / 2)
    return result

Внимательный читатель заметит, что произведения отличаются друг от друга только одним множителем, поэтому их можно не пересчитывать несколько раз. Вместо этого можно посчитать произведение всех двучленов, а затем получать каждый многочлен после однократного деления полного произведения на $y - y_i$ . Для этого можно воспользоваться схемой Горнера.

Ну а теперь можно коротко и ясно описать интегрирование:

def integrate_lagrange(f, a, b, n):
    '''
        Интегрирование полиномом Лагранжа степени n
    '''
    # создаём разбиение отрезка [-1, 1]
    ys = [-1 + 2 * float(i) / n for i in range(n + 1)]

    # считаем веса:
    w = weights(ys)

    # считаем ответ
    x = lambda y: (b + a) / 2 + y * (b - a) / 2 # отображает разбиение ys на
                                                # промежуток интегрирования
    s = sum(w[i] * f(x(y)) for i, y in enumerate(ys))
    return s * (b - a)

Обратите внимание, что пределы интегрирования нужно передавать как числа с плавающей точкой, то есть не 0 или 1, а 0.0 и 1.0. Это особенность python2.