5.2. Метод Ромберга

Метод Ромберга (Romberg’s method) позволяет быстро интегрировать с заданной точностью. Метод основывается на применении экстраполяции Ричардсона к формуле трапеций при удвоении числа узлов.

5.2.1. Экстраполяция Ричардсона

Частный случай этого метода мы рассматривали при выводе Формулы Симпсона.

[89] Рассмотрим вычисление интеграла \(I_h\) на равномерной сетке с шагом \(h\). На каждом интервале отрезка применяется одна и та же квадратурная формула. Для некоторых формул (например, для формулы трапеций) удаётся получить вид погрешности \(I - I_h\). Предположим, что это разложение выглядит так

(5.10)\[I_h = I + a_1 h^{\alpha_1} + a_2 h^{\alpha_2} + ... + a_m h^{\alpha_m} + O(h^{\alpha_{m+1}}).\]

Здесь \(0 < \alpha_1 < \alpha_2 < ... < \alpha_m\) – степени разложения, а \(a_i\) его ненулевые коэффициенты.

Теперь рассмотрим (5.10) на последовательности сеток с шагами \(h_0 = h\), \(h_1 = qh\), \(h_2 = q^2 h\), …, \(h_k = q^k h\), \(k=0,...,m\), \(0 < q < 1\). Тогда, согласно (5.10), запишем

(5.11)\[I_{h_k} = I + a_1 h_k^{\alpha_1} + a_2 h_k^{\alpha_2} + ... + a_m h_k^{\alpha_m} + O(h_k^{\alpha_{m+1}}),\quad k=0,...,m.\]

Введём обозначение \(I^{(1)}_{h_k} = I_{h_k}\) и запишем для двух сеток

\[\begin{split}I^{(1)}_{h_{k-1}} = I + a_1 h_{k-1}^{\alpha_1} + O(h^{\alpha_2}_{k-1})\\ I^{(1)}_{h_k} = I + a_1 h_k^{\alpha_1} + O(h^{\alpha_2}_{k}).\end{split}\]

Скомбинируем выражения так, чтобы слагаемое при \(a_1\) занулить, получим

(5.12)\[I = I^{(2)}_{h_{k-1}} + O(h^{\alpha_2}_{k-1}),\]

где

\[I^{(2)}_{h_{k-1}} = \frac{I^{(1)}_{h_k} - q^{\alpha_1}I^{(1)}_{h_{k-1}}}{1 - q^{\alpha_1}}.\]

Таким образом, мы выполнили один шаг экстраполяции Ричардсона. Важно отметить, что согласно (5.12) \(I^{(2)}\) приближает \(I\) лучше, чем \(I^{(1)}\), у которых ошибка порядка \(h^{\alpha_1}\).

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

(5.13)\[\begin{split}\begin{align} I_{h_{k-1}}^{(j+1)} = \frac{I^{(j)}_{h_k} - q^{\alpha_j}I^{(j)}_{h_{k-1}}}{1 - q^{\alpha_j}},\\ j = 1,..., m, \quad k = 1, 2, ..., m - j + 1, \\ I_{h_k}^{(1)} = I_{h_k},\quad k =0,1,...,m. \end{align}\end{split}\]

Точность приближения для (5.13) составляет

\[I = I_{h_k}^{(j)} + O(h_k^{\alpha_j}).\]

Важно отметить структуру формулы (5.13): \((j+1)\)-ое приближение на сетке с шагом \(h_{k-1}\) можно получить по \((j)\)-приближениям на сетках с шагами \(h_k\) и \(h_{k-1}\).

\ j  1    2    3   
k --------------
0   01 → 02 → 03
       ↗    ↗
1   11 → 12
       ↗
2   21

Теперь рассмотрим частный случай экстраполяции Ричардсона, метод Ромберга.

Положим в основу экстраполяции (5.13) формулу трапеций \(I_{h_k}^{(1)}=T_f(h_k)\) и будем мельчить сетку вдвое (\(q=0.5\)). Для составной формулы трапеций разложение (5.11) содержит только чётные степени, т.е. \(\alpha_j = 2 j\). Тогда, из (5.13) получим

(5.14)\[I_{h_{k-1}}^{(j+1)} = \frac{4^j I^{(j)}_{h_k} - I^{(j)}_{h_{k-1}}}{4^j - 1}.\]

Прежде чем переходить к реализации, рассмотрим особенность формулы трапеции при измельчении сетки, это позволит избежать перевычислений в (5.14).

5.2.2. Удвоение числа узлов

Рассмотрим равномерную сетку \(a = x_0 < x_1 < x_2 < x_3 < x_4 = b\) из \(n=4\) интервалов с шагом \(h_4\) и сетку из \(n=2\) интервалов \(a = x_0 < x_2 < x_4 =b\) с шагом \(2h_4\).

На крупной сетке формула трапеций даёт

(5.15)\[T_f(n=2) = 2h_4 \Big[0.5 \big(f(x_0) + f(x_4)\big) + f(x_2)\Big].\]

А на сетке вдвое мельче

(5.16)\[T_f(n=4) = h_4 \Big[ 0.5 \big(f(x_0) + f(x_4)\big) + f(x_1) + f(x_2) + f(x_3) \Big].\]

Сравнивая (5.15) и (5.16), обнаруживаем их связь

\[T_f(n=4) = 0.5 T_f(n=2) + h_4 \big(f(x_1) + f(x_3)\big).\]

Таким образом, мы в частном случае получили, что вычислив интеграл на крупной сетке, мы можем вычислить его на сетке вдвое мельче, используя только значения в новых узлах (\(x_1\), \(x_3\)).

Это верно и в общем случае

(5.17)\[T_f(2n) = 0.5 T_f(n) + h_{2n} \sum_{i=1}^n f(x_{2i-1}), \quad h_{2n} = \frac{b-a}{2n},\quad x_{2i-1} = a + h_{2n} (2i-1).\]

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

(5.18)\[|T_f(2n) - T_f(n)| < \text{tol}.\]

Например, при \(\text{tol} = 10^{-6}\) мы можем добиться точности до 6 знаков после запятой. Конечно, ошибки округления и обусловленность задачи никуда не делись.

Также отметим, что выражение (5.18) является абсолютной, а не относительной точностью.

Реализация может быть такой

Функция 5.10 (trapezoid_tol)

Формула трапеций с контролем точности

"""
Вычисляет интеграл ∫`f`dx на [`a`, `b`] с точностью `atol` по формуле трапеций,
удваивая число разбиений интервала, но не более `maxstep` раз.
Возвращает значение интеграла.
"""
function trapezoid_tol(f, a, b; atol=1e-3, maxstep::Integer=100)
    nc, hc = 1, float(b - a)
    Tc = hc * (f(a) + f(b)) / 2
    for step in 1:maxstep
        Tp, np = Tc, nc
        hc /= 2
        nc *= 2
        Tc = Tp / 2 + hc * sum(f, (a + hc*(2i-1) for i in 1:np))
        abs(Tc - Tp) < atol && return Tc
    end
    error("Точность не удовлетворена.")
end

5.2.3. Реализация

Можно обойтись и двумя векторами.

Итак, хранить \(I_{h_{k-1}}^{(j)}\) будем в элементе матрицы I[k, j].

Первый элемент I[1, 1] заполним вручную: он соответствует формуле трапеций на сетке из одного интервала (2 узла). Затем, пока не достигнем необходимой точности, будем делать следующее

1. Измельчим сетку вдвое и подсчитаем интеграл по формуле трапеций на ней (5.17), так получим из I[k, 1] элемент I[k+1, 1];

2. Вычислим по доступной побочной диагонали слева направо все элементы от I[k, 2] вплоть до I[1, k+1] по формуле перехода (5.14);

3. Будем продолжать 1 и 2 до тех пор, пока приближения I[1, k+1] и I[2, k] станут неотличимы в заданной точности atol, или пока не исчерпаем количество измельчений maxstep.

Пример заполнения матрицы I. Допустим, максимальное измельчение maxstep = 4, т.е. от \(h=b-a\) до \(h = (b-a) / 2^4\). Ниже приведена схема вычислений и обозначения

• ⋅: элемент матрицы, который не пригодится в вычислениях;

• □: элемент, который, возможно, предстоит определить;

• ○ и ●: уже вычисленные элементы, причём ● сравниваются между собой в конце итерации;

• стрелка вниз показывает переход по умельчению сетки (шаг 1);

• остальные стрелки показывают переход по формуле (шаг 2).

Этот алгоритм является типичной задачей динамического программирования.

     1 итерация             2 итерация           3 итерация             4 итерация

○ → ●   □   □   □      ○   ○ → ●   □   □      ○   ○   ○ → ●   □      ○   ○   ○   ○ → ●
↓ ↗                          ↗                          ↗                          ↗
●   □   □   □   ⋅      ○ → ●   □   □   ⋅      ○   ○ → ●   □   ⋅      ○   ○   ○ → ●   ⋅
                       ↓ ↗                          ↗                          ↗
□   □   □   ⋅   ⋅      ○   □   □   ⋅   ⋅      ○ → ○   □   ⋅   ⋅      ○   ○ → ○   ⋅   ⋅
                                              ↓ ↗                          ↗
□   □   ⋅   ⋅   ⋅      □   □   ⋅   ⋅   ⋅      ○   □   ⋅   ⋅   ⋅      ○ → ○   ⋅   ⋅   ⋅
                                                                     ↓ ↗
□   ⋅   ⋅   ⋅   ⋅      □   ⋅   ⋅   ⋅   ⋅      □   ⋅   ⋅   ⋅   ⋅      ○   ⋅   ⋅   ⋅   ⋅

Функция 5.11 (romberg)

Метод Ромберга с контролем точности

"""
Вычисляет интеграл ∫`f`dx на отрезке [`a`, `b`] методом Ромберга.
Разбивает отрезок пополам не более `maxstep > 0` раз.
Возвращает значение интеграла, если приближения отличаются не более чем на `atol`.
"""
function romberg(f, a, b; atol=1e-6, maxstep::Integer=100)
    maxstep = max(1, maxstep)  # хотя бы одно разбиение
    I = Matrix{Float64}(undef, maxstep+1, maxstep+1)
    I[1, 1] = (b - a) * (f(a) + f(b)) / 2
    for i in 2:maxstep+1
        let hc = (b - a) / 2^(i-1), np = 2^(i-2)
            I[i, 1] = I[i-1, 1] / 2 + hc * sum(f, (a + hc * (2i-1) for i in 1:np))
        end
        for k in i-1:-1:1
            I[k, i-k+1] = (2^i*I[k+1, i-k] - I[k, i-k]) / (2^i - 1)
        end
        abs(I[1, i] - I[2, i-1]) < atol && return I[1, i]
    end
    error("Точность не удовлетворена.")
end

Пример 5.12

Сравним время вычисления интеграла

\[\int_0^3 x e^{\sin 2x}\]

нашими реализациями trapezoid_tol и romberg.

julia> using BenchmarkTools

julia> foo(x) = x * exp(sin(2x));

julia> a, b, atol = 0, 3, 1e-6;

julia> @benchmark trapezoid_tol($foo, a, b; atol=$atol)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  121.110 μs … 452.256 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     121.582 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   122.370 μs ±   6.715 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▇█▄▄▁▁                                                        ▁
  ████████▇▆▆▅▆▄▅▃▃▄▄▁▄▅▄▅▅▅▆▅▆▅▅▄▅▄▃▅▅▄▄▄▁▄▄▄▄▄▄▄▄▄▄▃▃▄▅▅▅▁▃▄▅ █
  121 μs        Histogram: log(frequency) by time        144 μs <

 Memory estimate: 272 bytes, allocs estimate: 17.

julia> @benchmark romberg($foo, a, b; atol=$atol)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):   8.434 μs …   3.711 ms  ┊ GC (min … max):  0.00% … 99.35%
 Time  (median):     13.210 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   20.341 μs ± 103.650 μs  ┊ GC (mean ± σ):  20.16% ±  4.08%

  ▃   ▃█▆▄▃▄▄▃▂▁                                               ▁
  █▅▆▅███████████▇█▇▇▇▇▇▇▆▆█▇▇▇▇▇▇▇▆▆▆▅▅▅▅▄▃▅▄▅▃▃▄▄▅▅▄▄▃▄▃▃▄▃▂ █
  8.43 μs       Histogram: log(frequency) by time      57.9 μs <

 Memory estimate: 79.89 KiB, allocs estimate: 6.

Получаем, что на данной машине реализация метода Ромберга (несмотря на аллоцирование матрицы), примерно в 6 раз (≈ 120 мкс / 20 мкс) быстрее вычислила интеграл.

Если немного модифицировать код и посмотреть, сколько разбиений сетки произошло при вычислениях, то получим, что trapezoid_tol потребовалось 12 разбиений сетки, а romberg всего 8. Это означает, что trapezoid_tol вычисляла функцию примерно на \(2^{12} - 2^8 \approx 3800\) раз больше.

В свою очередь, одно вычисление функции foo стоит около 14 нс.

julia> x = 1.5;

julia> @benchmark foo($x)
BenchmarkTools.Trial: 10000 samples with 998 evaluations.
 Range (min … max):  13.746 ns … 171.441 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     14.281 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   14.639 ns ±   3.666 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▂█▆▃█▄▂▁▁                                                    ▂
  ██████████▇▆▆▆▆▆▅▇▄▅▅▆▅▅▄▅▄▃▁▅▃▄▃▁▄▄▄▅▁▃▅▄▄▁▁▄▄▁▃▁▅▄▅▃▄▃▅▁▄▄ █
  13.7 ns       Histogram: log(frequency) by time      24.7 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

Значит, romberg только на вычислениях функции «сэкономил» примерно 3800 × 14 нс ≈ 53 мкс, что составляет около 40% времени вычисления trapezoid_tol.