Общие соображения

При использовании многочлена Лагранжа получается квадратурная формула вида

которая точна для многочленов степени . Но может быть, можно подобрать такие узлы и веса, при которых формула будет точно для многочленов более высоких степеней?

Без уменьшения общности можно рассматривать только интегралы на . Для начала попробуем “сделать” формулу точной для всех многочленов степени . Очевидно, что если формула точна хотя бы для одного такого многочлена, то она точна и для всех остальных. Поэтому рассмотрим многочлен

для которого, согласно квадратурной формуле, получаем

Отсюда получается, что для точности формулы для многочленов порядка на узлы нужно наложить условие:

Теперь попробуем ещё улучшить формулу и сделать её точной для многочленов степени . В качестве испытуемого многочлена выберем такой, у которого правая часть формулы наверняка будет нулевой. Например, можно использовать :

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

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

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

Это равенство означает, что многочлен ортогонален всем многочленам меньших степеней с весом 1 на отрезке . Семейство таких многочленов называется ортогональными многочленами Лежандра .

Это означает, что

откуда следует, что узлами разбиения являются корни многочлена Лежандра .

Определение узлов разбиения

Для определения корней можно воспользоваться следующим алгоритмом:

def lejendre_roots(n):
    p = lejendre(n)
    f = to_func(p)
    d = to_func(derivative(p))
    xs = []
    for i in range(n):
        x = cos(pi * (4 * i + 3) / (4 * n + 2)) # начальное значение
        delta = 1 # произвольное значение, чтобы зайти в цикл
        while abs(delta) > 1e-8:
            delta = f(x) / d(x)
            x -= delta
        xs.append(x) # добавляем корень к ответу и ищем следующий
    return xs

где legendre(n) возвращает массив коэффициентов полинома Лежандра, derivative(p) – массив коэффициентов производной многочлена, а to_func(p) превращает массив коэффициентов многочлена в функцию:

def lejendre(n):
    if n == 0:
        return [1]      # P₀(x) = 1
    elif n == 1:
        return [0, 1]   # P₁(x) = x

    # далее обрабатываем остальные случаи с помощью рекуррентной формулы
    a, b = [1], [0, 1]

    for i in range(1, n):
        a, b = b, a
        b = b + [0.0, 0.0]

        b[0] *= -i / (i + 1.0)
        for j in range(1, i + 2):
            b[j] = (2 * i + 1.0) / (i + 1.0) * a[j - 1] - i / (i + 1.0) * b[j]
    return b


def derivative(p):
    return [a * (i + 1) for i, a in enumerate(p[1:])]


def to_func(p):
    return lambda x: sum(a * x ** i for i, a in enumerate(p))

Определение весов

Теперь нужно определить веса. Это можно сделать несколькими способами.

Способ первый: очевидный

Взять функцию для нахождения весов в методе Лагранжа и посчитать его на полученных узлах.

Способ второй: менее очевидный

Так как формула точна для многочленов степеней :

Решение этой системы – искомые веса.

Способ третий: совсем нетривиальный

Внимание! В этом способе используются свойства многочленов Лежандра! Ни в коем случае не пытайтесь повторить это в домашних условиях!

Из статьи про метод полинома Лагранжа мы уже знаем, что

где

Следовательно,

Остановимся подробнее на интеграле. Для его взятия рассмотрим более общий интеграл

который сводится к искомому с точностью до постоянного множителя подстановкой . Воспользуемся рекуррентным соотношением для многочленов Лежандра

для преобразования числителя и лучшего понимания происходящего:

Ясное дело, что первый интеграл равен нулю, откуда

Ну и напоследок

откуда окончательно