Алгоритм Ромберга

Если подынтегральная функция "хорошая", т. е. не меняется на интервале интегрирования слишком быстро и не обращается на нем в бесконечность, то численное вычисление интеграла не принесет никаких неприятных сюрпризов. Приведем основные идеи итерационного алгоритма Ромберга, который применяется для интегрирования таких функций с параллельным апостериорным контролем ошибки интегрирования.

• Сначала строится несколько интерполирующих полиномов, которые заменяют на интервале интегрирования подынтегральную функцию f(x). В качестве первой итерации полиномы вычисляются по 1, 2 и 4 интервалам. Например, первый полином, построенный по 1 интервалу - это просто прямая линия, проведенная через две граничные точки интервала интегрирования, второй - квадратичная парабола и т. д.
• Интеграл от каждого полинома с известными коэффициентами легко вычисляется аналитически. Таким образом, определяется последовательность интегралов от интерполирующих полиномов: I₁, I₂, I₄,… Например, по правилу трапеций I₁=(b-a)*(f(a)+f(b))/2 и т. д.
• Из-за интерполяции по разному числу точек вычисленные интегралы I₁, I₂,… несколько отличаются друг от друга. Причем, чем больше точек используется для интерполяции, тем интеграл от интерполяционного полинома ближе к искомому интегралу, стремясь к нему в пределе бесконечного числа точек. Поэтому определенным образом осуществляется экстраполяция последовательности I₁, I₂, I₄,… до нулевой ширины элементарного интервала. Результат этой экстраполяции J₀ принимается за нулевое приближение к вычисляемому интегралу.
• Осуществляется переход к новой итерации с помощью еще более частого разбиения интервала интегрирования, добавления нового члена последовательности интерполирующих полиномов и вычисления нового (N-го) приближения Ромберга J_N.
• Чем больше количество точек интерполяции, тем ближе очередное приближение Ромберга к вычисляемому интегралу и, соответственно, тем меньше оно отличается от приближения предыдущей итерации. Как только разница между двумя последними итерациями |J_N-J_N-1| становится меньше погрешности ( абсолютной, или в единицах |J_N|), итерации прерываются, и J_N принимается в качестве результата интегрирования.