Ctrl+K

Практическая оценка погрешности

Практическая оценка погрешности#

Не всегда на практике удается воспользоваться априорной оценкой вида

\[ \varepsilon \leqslant C M_p h^p, \]

так как не всегда удается оценить \(M_p\) для подынтегральной функции. В этом случае, обычно, пользуются правилом Рунге для определения подходящего шага \(h\).

Правило Рунге#

Предположим, что нам известен порядок метода \(p\), которым мы хотим найти значение интеграла с заданной точностью \(\varepsilon\). Тогда результат, вычисленный этим методом с шагом \(h\) будет иметь вид

\[ I_h = I^* + \underbrace{C_1 h^p + C_2 h^{p+1} + \dots}_\text{ошибка интегрирования}. \]

Здесь \(I^*\) — точное значение интеграла. Для достаточно малых \(h\) можно записать

\[ I_h = I^* + C_h h^p, \]

где \(C_h\) — почти константа, то есть слабо зависит от \(h\).

Вычислим интеграл с шагом \(h\) и с шагом \(h/2\):

\[\begin{split} I_h = I^* + C_1 h^p + C_2 h^{p+1} + \dots\\ I_{h/2} = I^* + C_1 \left(\frac{h}{2}\right)^p + C_2 \left(\frac{h}{2}\right)^{p+1} + \dots \end{split}\]

Тогда

\[\begin{split} I_h - I_{h/2} = C_h h^p - C_{h/2} \left(\frac{h}{2}\right)^p \approx (2^p - 1) C_{h/2} \left(\frac{h}{2}\right)^p = (2^p - 1) (I_{h/2} - I^*)\\ \end{split}\]

Таким образом, погрешность интегрирования на сетке \(h/2\) может быть оценена как

\[ I_{h/2} - I^* \approx \frac{I_h - I_{h/2}}{2^p - 1}. \]

Алгоритм оценки погрешности#

Алгоритм применения правила Рунге следующий:

  1. Задать небольшое число интевалов, скажем \(n = 10\)

  2. Вычислить \(I_h, I_{h/2}\) и оценить ошибку

\[ \Delta_{h/2} = \frac{I_h - I_{h/2}}{2^p - 1}. \]

  1. Если \(|\Delta_{h/2}| < \varepsilon\), прекратить вычисления. Иначе измельчить сетку вдвое \(h := h/2\) и перейти к шагу 2.

  2. Ответом служит \(I_{h/2}\).

import numpy as np

def simpson(f, h):
    assert (len(f)-1) % 2 == 0
    return h/3. * (f[0] + 4 * sum(f[1:-1:2]) + \
                   2 * sum(f[2:-2:2]) + f[-1])

def runge(f, a, b, eps=1e-12):    
    n = 4
    Ih = simpson(f(np.linspace(a, b, n+1)), (b-a) / n)
    while True:
        Ih_2 = simpson(f(np.linspace(a, b, 2*n+1)), (b-a) / (2*n))
        Dh_2 = (Ih - Ih_2) / (2**4 - 1) # p = 4 для Симпсона
        print('I(h) = %.16f, err(h) = %.6e' % (Ih_2, Dh_2))
        if abs(Dh_2) < eps:
            break
        n *= 2; Ih = Ih_2
        if n > 10000: print('Too large n'); break
    return Ih_2

def f(x): return 1 / (1 + x**2)
a = 0; b = 0.5; exact = np.arctan(b)

runge(f, a, b) - exact

I(h) = 0.4636479223346336, err(h) = 3.157185e-07
I(h) = 0.4636476285453064, err(h) = 1.958596e-08
I(h) = 0.4636476102217171, err(h) = 1.221573e-09
I(h) = 0.4636476090771032, err(h) = 7.630759e-11
I(h) = 0.4636476090055746, err(h) = 4.768578e-12
I(h) = 0.4636476090011042, err(h) = 2.980246e-13

2.9809488211185453e-13

Условия корректности применения правила Рунге#

Для проверки корректности использования правила Рунге достаточно проверять условия

  • что метод фактически сходится с порядком \(p\),

\[ \frac{\Delta_{h}}{\Delta_{h/2}} \approx 2^p \]

  • что константа \(C_h\) слабо зависит от \(h\):

\[ C_h = \frac{\Delta_h}{h^p} \to C = \mathrm{const} \]

  • что выполнено эмпирическое правило

\[ \left|2^p \cdot \frac{I_{h/2} - I_h}{I_h - I_{2h}}-1\right|<0.1 \]
def runge_checks(f, a, b, eps=1e-12):
    n = 4
    Ih = simpson(f(np.linspace(a, b, n+1)), (b-a) / n)
    Dh = None;
    while True:
        h_2 = (b-a) / (2*n)
        Ih_2 = simpson(f(np.linspace(a, b, 2*n+1)), h_2)
        Dh_2 = (Ih - Ih_2) / (2**4 - 1) # p = 4 для Симпсона
        Ch_2 = Dh_2 / h_2**4; ps = np.log2(Dh / Dh_2) if Dh != None else np.nan
        print('I(h) = %.16f, err(h) = %.6e, p* = %4.2f, C = %.6e' % \
              (Ih_2, Dh_2, ps, Ch_2))
        if abs(Dh_2) < eps:
            break
        n *= 2; Ih = Ih_2; Dh = Dh_2
        if n > 10000: print('Too large n'); break
    return Ih_2

def f(x): return 1 / (1 + x**2)
a = 0; b = 0.5; exact = np.arctan(b)

runge_checks(f, a, b) - exact

I(h) = 0.4636479223346336, err(h) = 3.157185e-07, p* =  nan, C = 2.069093e-02
I(h) = 0.4636476285453064, err(h) = 1.958596e-08, p* = 4.01, C = 2.053736e-02
I(h) = 0.4636476102217171, err(h) = 1.221573e-09, p* = 4.00, C = 2.049459e-02
I(h) = 0.4636476090771032, err(h) = 7.630759e-11, p* = 4.00, C = 2.048366e-02
I(h) = 0.4636476090055746, err(h) = 4.768578e-12, p* = 4.00, C = 2.048089e-02
I(h) = 0.4636476090011042, err(h) = 2.980246e-13, p* = 4.00, C = 2.048009e-02

2.9809488211185453e-13

Экстраполяция Ричардсона и алгоритм Ронберга#

Вспомним равенство, которое ранее получили:

\[ I_{h/2} - I^* \approx \frac{I_h - I_{h/2}}{2^p - 1}. \]

Переставляя слагаемые, понимаем, что не изменяя сетку и не делая никаких новых вычислений функции, мы получили метод с бОльшим порядком сходимости! Такой метод уточнения точности называется экстраполяцией Ричардсона.

\[ I^* \approx I_{h/2} + \frac{I_{h/2} - I_h}{2^p - 1} + O(h^{p+1}) \]

В частности, если мы подставим в выражение выше квадратуру для формулы трапеций, мы получим формулу Симпсона.

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

Содержание