Интегралы: полный курс / Численное интегрирование / урок 43 из 64

Формула Симпсона и правило Рунге

⏱ 45 мин · 🧩 заданий: 7

Идея: параболы вместо отрезков

Трапеции заменяют функцию на каждом отрезке прямой. Возьмём пары отрезков и через три точки $(x_0, f_0)$ , $(x_1, f_1)$ , $(x_2, f_2)$ проведём параболу. Интеграл от неё:

\int_{x_0}^{x_2}P(x)\,dx = \frac h3\left(f_0 + 4f_1 + f_2\right).

(Проверьте на $x_0 = -h$ , $x_1 = 0$ , $x_2 = h$ для $P = 1, x, x^2$ .)

Составная формула Симпсона ( $n$ чётное):

S_n = \frac h3\left(f_0 + 4f_1 + 2f_2 + 4f_3 + 2f_4 + \dots + 4f_{n-1} + f_n\right).

Коэффициенты: $1, 4, 2, 4, 2, \dots, 4, 1$ .

Точность

R = -\frac{(b - a)h^4}{180}f^{(4)}(\xi), \qquad |R| \le \frac{(b - a)h^4}{180}\max|f^{(4)}|.

Четвёртый порядок: удвоение $n$ уменьшает ошибку в $16$ раз. Неожиданный бонус: формула точна не только для парабол, но и для кубических многочленов ( $f^{(4)} = 0$ ) — ошибки на $x^3$ по симметрии сокращаются.

Связь с предыдущим уроком: $S_{2n} = \frac{2M_n + T_n}{3}$ — взвешенное среднее средних прямоугольников и трапеций, при котором главные ошибки ( $\frac{h^2}{24}f''$ и $-\frac{h^2}{12}f''$ ) взаимно уничтожаются.

import numpy as np
import vis

f, exact = (lambda x: 1 / (1 + x**2)), np.pi / 4

def trap(n):
    x = np.linspace(0, 1, n + 1); h = 1 / n
    return h * (np.sum(f(x)) - (f(0) + f(1)) / 2)

def simpson(n):
    x = np.linspace(0, 1, n + 1); h = 1 / n
    return h / 3 * (f(x[0]) + f(x[-1]) + 4 * np.sum(f(x[1:-1:2])) + 2 * np.sum(f(x[2:-1:2])))

ns = [2, 4, 8, 16, 32, 64]
et = [abs(trap(n) - exact) for n in ns]
es = [abs(simpson(n) - exact) for n in ns]
vis.line(np.log2(ns), {"трапеции": et, "Симпсон": es}, logy=True, title="Симпсон: наклон −4 против −2 у трапеций", xlabel="log₂ n")
for n, a, b in zip(ns, et, es):
    print(f"n = {n:3}: трапеции {a:.2e}, Симпсон {b:.2e}")

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

На практике производные неизвестны, и оценка через $\max|f^{(4)}|$ бесполезна. Выход — посчитать дважды, с шагом $h$ и $\frac h2$ . Если метод имеет порядок $p$ , то $I - I_h \approx Ch^p$ , и

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

Для трапеций делим разность на $3$ , для Симпсона — на $15$ . Это апостериорная оценка: её дают сами вычисления.

Экстраполяция Ричардсона. Прибавив оценку ошибки, получим более точное значение:

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

Для трапеций это в точности формула Симпсона! Повторяя приём (трапеции → Симпсон → формула 6-го порядка → …), получаем метод Ромберга.

Адаптивное интегрирование и Гаусс

Современные программы не дробят отрезок равномерно: где функция «спокойная», шаг крупный, где резко меняется — мелкий. Ошибку на каждом куске оценивают по Рунге.

Квадратуры Гаусса выбирают не только веса, но и узлы: $n$ точек дают точность для многочленов степени $2n - 1$ . Формулы Гаусса и их адаптивные варианты лежат в основе scipy.integrate.quad и mpmath.quad.

import numpy as np

# Квадратура Гаусса–Лежандра: n узлов на [−1, 1], переводим на [0, 1]
for n in [2, 3, 5, 8]:
    nodes, weights = np.polynomial.legendre.leggauss(n)
    x = (nodes + 1) / 2
    approx = np.sum(weights * np.exp(-x**2)) / 2
    print(f"Гаусс, {n} узлов: ∫₀¹ e^(−x²) dx ≈ {approx:.12f}")
print("эталон (erf):         ", f"{0.746824132812427:.12f}")

Восемь значений функции — двенадцать верных знаков.

Итог

Симпсон: $\frac h3(f_0 + 4f_1 + 2f_2 + \dots + 4f_{n-1} + f_n)$ , порядок 4, точен для кубических многочленов.
Правило Рунге: ошибка $\approx\frac{I_{h/2} - I_h}{2^p - 1}$ ; экстраполяция Ричардсона повышает порядок.
Квадратуры Гаусса: $n$ узлов — точность для степени $2n - 1$ .

Практика решено 0 из 7

Вы не вошли — ответы проверяются, но прогресс не сохранится. Войти или зарегистрироваться.

Задание 1 Ввод формулы

Формула Симпсона для $\int_0^2 x^2\,dx$ при $n = 2$ ( $h = 1$ ): $\frac h3(f_0 + 4f_1 + f_2)$ .

Как вводить формулы

Степень — ^ (x^2, e^(3x)), умножение можно не писать (2x, x sin(x)), дробь — / (1/(x+1)), корень — sqrt(x), модуль — |x| или abs(x). Функции: sin cos tg ctg arcsin arccos arctg arcctg sh ch ln lg exp, константы pi и e. Аргументы функций лучше брать в скобки: sin(2x). Для неопределённого интеграла допишите + C. Индексы пишутся слитно: m1, v0; греческие буквы — словом или символом: alpha или α. Десятичная дробь — через точку или запятую, большие числа — через 10^: 6,67*10^(-11).

Эталонное решение (ваше может отличаться):

8/3

Задание 2 Ввод формулы

Симпсон для $\int_0^1 \frac{dx}{1 + x^2}$ при $n = 2$ . Введите до пяти знаков (точное значение $\frac\pi4 \approx 0{,}78540$ ).

Как вводить формулы

Эталонное решение (ваше может отличаться):

0.78333

Задание 3 Ввод формулы

Во сколько раз уменьшится ошибка формулы Симпсона при увеличении $n$ с 10 до 20?

Как вводить формулы

Эталонное решение (ваше может отличаться):

Задание 4 Ввод формулы

Правило Рунге: трапеции дали $T_{10} = 1{,}7197$ и $T_{20} = 1{,}7186$ (порядок $p = 2$ ). Уточнённое значение $T_{20} + \frac{T_{20} - T_{10}}{2^p - 1}$ = ?

Как вводить формулы

Эталонное решение (ваше может отличаться):

1.718233

Задание 5 Ввод формулы

Оценка погрешности Симпсона $|R| \le \frac{(b - a)h^4}{180}\max|f^{(4)}|$ . Для $\int_0^1 e^x\,dx$ и $h = 0{,}1$ найдите эту оценку (до двух значащих цифр).

Как вводить формулы

Эталонное решение (ваше может отличаться):

e/180 * 0.0001

Задание 6 Один вариант ответа

Двухточечная формула Гаусса $\int_{-1}^1 f\,dx \approx f\left(-\frac{1}{\sqrt3}\right) + f\left(\frac{1}{\sqrt3}\right)$ точна для многочленов степени до…

1 2 3 4

Задание 7 Программа на Python

Симпсон и Рунге. Напишите simpson(f, a, b, n) — составную формулу Симпсона (n чётное; если нечётное — выбросьте ValueError) и runge(f, a, b, n), которая возвращает пару (S_2n, err): значение по Симпсону с 2n отрезками и оценку его ошибки по Рунге $\frac{|S_{2n} - S_n|}{15}$ .

def simpson(f, a, b, n):
    return 0.0

def runge(f, a, b, n):
    return 0.0, 0.0

Какие тесты проверяют код (1)

import math
assert abs(simpson(lambda x: x**3 - 2 * x, 0, 2, 2) - 0) < 1e-12, "Симпсон точен для кубических многочленов"
assert abs(simpson(lambda x: 1 / (1 + x * x), 0, 1, 2) - 4.7 / 6) < 1e-12
try:
    simpson(math.sin, 0, 1, 3)
    assert False, "Для нечётного n нужно выбросить ValueError"
except ValueError:
    pass
val, err = runge(math.exp, 0, 1, 4)
true_err = abs(val - (math.e - 1))
assert true_err < 1e-5 and 0.5 < err / true_err < 2, "Оценка Рунге близка к настоящей ошибке"

Эталонное решение (ваше может отличаться):

def simpson(f, a, b, n):
    if n % 2:
        raise ValueError("n должно быть чётным")
    h = (b - a) / n
    s = f(a) + f(b)
    s += 4 * sum(f(a + k * h) for k in range(1, n, 2))
    s += 2 * sum(f(a + k * h) for k in range(2, n, 2))
    return s * h / 3

def runge(f, a, b, n):
    s1, s2 = simpson(f, a, b, n), simpson(f, a, b, 2 * n)
    return s2, abs(s2 - s1) / 15

Войдите, чтобы отмечать пройденные уроки.