Численное интегрирование: формулы Гаусса
Общие соображения
При использовании многочлена Лагранжа получается квадратурная формула вида
которая точна для многочленов степени . Но может быть, можно подобрать такие узлы и веса, при которых формула будет точно для многочленов более высоких степеней?
Без уменьшения общности можно рассматривать только интегралы на . Для начала попробуем “сделать” формулу точной для всех многочленов степени . Очевидно, что если формула точна хотя бы для одного такого многочлена, то она точна и для всех остальных. Поэтому рассмотрим многочлен
для которого, согласно квадратурной формуле, получаем
Отсюда получается, что для точности формулы для многочленов порядка на узлы нужно наложить условие:
Теперь попробуем ещё улучшить формулу и сделать её точной для многочленов степени . В качестве испытуемого многочлена выберем такой, у которого правая часть формулы наверняка будет нулевой. Например, можно использовать :
Продолжая подбирать многочлены и выписывать нулевые интегралы, мы можем существенно конкретизировать узлы и увеличить точность интегрирования.
Однако, для многочлена точность формулы уже гарантировать нельзя. Нетрудно построить многочлен степени , для которого формула будет неточна при любом выборе узлов:
Поэтому максимальная степень многочлена, для которой формула может быть верна равна . Причём, для того, чтобы она была верна, должны выполняться следующие условия:
Это равенство означает, что многочлен ортогонален всем многочленам меньших степеней с весом 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))
Определение весов
Теперь нужно определить веса. Это можно сделать несколькими способами.
Способ первый: очевидный
Взять функцию для нахождения весов в методе Лагранжа и посчитать его на полученных узлах.
Способ второй: менее очевидный
Так как формула точна для многочленов степеней :
Решение этой системы – искомые веса.
Способ третий: совсем нетривиальный
Внимание! В этом способе используются свойства многочленов Лежандра! Ни в коем случае не пытайтесь повторить это в домашних условиях!
Из статьи про метод полинома Лагранжа мы уже знаем, что
где
Следовательно,
Остановимся подробнее на интеграле. Для его взятия рассмотрим более общий интеграл
который сводится к искомому с точностью до постоянного множителя подстановкой . Воспользуемся рекуррентным соотношением для многочленов Лежандра
для преобразования числителя и лучшего понимания происходящего:
Ясное дело, что первый интеграл равен нулю, откуда
Ну и напоследок
откуда окончательно