Правило Рунге оценки погрешности при интегрировании

Теорию рекомендую почитать здесь (стр. 15-18).

Перейдем к практике.

Требуется рассчитать определенный интеграл функции f(x) = ln(x + (x² + 4)^0.5) методом трапеций на произвольном участке [a, b] с заданной точностью е = 10^-5.

import matplotlib.pyplot as plt
import math
import numpy as np

#метод трапеций
def trapeze(y, ab, h):
    a, b = ab[0], ab[1]
    s = 0
    while round(a, 8) < b:
        s += 0.5 * h * (y(a) + y(a + h))
        a += h
    return s

def f(x):
    return math.log(x + (x ** 2 + 4) ** 0.5) #функция
ab = [0, 1.5] #участок интегрирования

#оценка погрешности по правилу Рунге
e = 0.00001 #требуемая точность
k = 2 #порядок точности квадратурной формулы
N = 1
ei = 0
while True:
    N += 1 #число разбиений
    h = (ab[1] - ab[0]) / N #шаг
    ei = (trapeze(f, ab, h / 2) - trapeze(f, ab, h)) / (2 ** k - 1)
    if ei < e: break

print(N, ei, trapeze(f, ab, h))

#график функции
xi = np.linspace(ab[0], ab[1], N + 1)
yi = [f(i) for i in xi]
plt.plot(xi, yi)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid()
plt.show()

Результат: при числе разбиений N = 22 площадь под графиком S = 1.5794027981035632 с точностью е = 9.686138265981489e-06.

Теперь найдем интеграл аналитически и сравним точности рассчетов методами трапеций и Ньютона ("3/8").

Код:

import matplotlib.pyplot as plt
import math
import numpy as np

#метод трапеций
def trapeze(y, a, b, N):
    h = (b - a) / N
    s = 0
    while round(a, 8) < b:
        s += 0.5 * h * (y(a) + y(a + h))
        a += h
    return s

#метод Ньютона
def nyton(y, a, b, N):
    h = (b - a) / N #шаг
    s = 0
    while round(a, 8) < b:
        s += (h / 8) * (y(a) + 3 * y(a + h / 3) + 3 * y(a + h / 1.5) + y(a + h))
        a += h
    return s

def f(x):
    return math.log(x + (x ** 2 + 4) ** 0.5) #функция
a, b = 0, 1.5 #участок интегрирования
N = 22

#интеграл заданной функции, посчитанный вручную
def integral(x):
    return x * math.log(x + (x ** 2 + 4) ** 0.5) - (x ** 2 + 4) ** 0.5

real_s = integral(b) - integral(a)
about_s = trapeze(f, a, b, N)
e = real_s - about_s
print(real_s)
print(about_s)
print(e)

#график функции
xi = np.linspace(a, b, N + 1)
yi = [f(i) for i in xi]
plt.plot(xi, yi)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid()
plt.show()

#оценка погрешностей по правилу Рунге для метода трапеций
k = 2 #порядок точности квадратурной формулы
N = 8
while N < 130:
    ei = (trapeze(f, a, b, 2 * N) - trapeze(f, a, b, N)) / (2 ** k - 1)
    print(N, ei)
    N = N * 2

#оценка погрешностей по правилу Рунге для метода Ньютона
k = 4 #порядок точности квадратурной формулы
N = 8
while N < 130:
    ei = (nyton(f, a, b, 2 * N) - nyton(f, a, b, N)) / (2 ** k - 1)
    print(N, ei)
    N = N * 2

На выходе получим зависимость точности е от числа разбиений N.

Для метода трапеций она будет иметь вид:

8 7.331217001693702e-05
16 1.8314912301642394e-05
32 4.577909425984572e-06
64 1.1444262219365935e-06
128 2.8610335989220914e-07

А для метода Ньютона:

8 -1.5556134454660273e-09
16 -9.70164541295541e-11
32 -6.060278205192541e-12
64 -3.787192781601334e-13
128 -2.3655151911346668e-14

Если присмотреться, можно увидеть интересную картину: соседние значения погрешностей для метода трапеций отличаются примерно в 4 раза, а для метода Ньютона - в 16 раз. Почему?

Все дело в точности расчетов обоими методами! Если обратиться к теории, которую я рекомендовал почитать в самом начале, найдем одну занятную формулу:

В левой части данной формулы находится разность интегралов, вычесленных одним и тем же методом с числами разбиений, отличающимися в 2 раза. А в правой части следует обратить внимание на множитель (0.5h)^k. Так вот, для метода трапеций порядок точности квадратурной формулы k = 2, стало быть множитель приобретает вид h^k/ 4, а для метода Ньютона k = 4, значит получим h^k/ 16. Вот и весь секрет!

Статьи

Правило Рунге оценки погрешности при интегрировании