задачи численного интегрирования приходится решать для
вместо функции, которую требуется проинтегрировать, интегрируется интерполяционный полином
пусть требуется определить значение интеграла функции f(x) на отрезке [a,b]. этот отрезок делится точками на n равных отрезков, каждый длиной
d = (b - a) / n
обозначим через y₁,y₂,...,yₙ значение функции f(x) в точках x₁,x₂,...,xₙ, каждая из которых принадлежит одному из элементарных отрезков:
далее составляем сумму:
S = y₁ * d + y₂ * d + .. + yₙ * d
эта сумма приближённо выражает интеграл для f(x) на [a,b]
S = (y₁ + y₂ + .. + yₙ) * d
если функцию на каждом из частичных отрезков аппроксимировать прямой, параллельной оси абсцисс и проходящей через начало, середину или конец или через любую вообще точку отрезка, то получим метод прямоугольников
площадь прямоугольника на каждом отрезке:
S = yi * d
полная формула в случае деления всего промежутка интегрирования на n отрезков одинаковой длины d:
S = (y₁ + y₂ + .. + yₙ) * dочевидно, стоит рассчитывать на бо́льшую точность если брать в качестве опорной точки для нахождения высоты точку посередине промежутка
ε ∝ y'(ζ) * n * d² / 4 ∝ y' * n * ((xn-x1)/n)²
ε ∝ Ο (1/n)если функцию на каждом из частичных отрезков аппроксимировать прямой, проходящей через конечные значения, то получим метод трапеций
площадь трапеции на каждом отрезке:
S = d * fi = d * yi + d * (yi+1 - yi) / 2
полная формула в случае деления всего промежутка интегрирования на отрезки одинаковой длины d:
S = (f₀ / 2 + f₁ + f₂ + .. + fn-1 + fₙ / 2) * d
ε ∝ y"(ζ) * n * d³ / 12
ε ∝ Ο (1/n²)множитель 1/12 появляется из
использовав три точки отрезка интегрирования, можно заменить подынтегральную функцию квадратичной параболой, проходящей через точки. обычно в качестве таких точек используют концы отрезка и его среднюю точку
выведем формулу Симпсона:
рассматриваем функцию
f(x) = a * x² + b * x + c (0)
рассмотрим отрезок [0, 2*h]. очевидно, что
f(0) = c (1)
f(h) = a * h² + b * h + c (2)
f(2*h) = 4 * a * h² + 2 * b * h + c (3)
2h
найдем ∫ f(x) dx
0
2h 2h
∫ f(x) dx = a * x³/3 + b * x²/2 + c * x | = (h / 3) * (8 * a * h² + 6 * b * h + 6 * c)
0 0
используем равенства (1-3), чтобы осуществить переход для выражения в скобках:
(1) 4 * (2) (3)
------------------------------------------------------------------------------------
с + 4 * ( a * h² + b * h + c ) + 4 * a * h² + 2 * b * h + c =
8 * a * h² + 6 * b * h + 6 * c =
f(0) + 4 * f(h) + f(2*h)
итак:
2h
∫ f(x) dx = (h / 3) * (f(0) + 4 * f(h) + f(2*h))
0
итак, формула для площади на отрезке d со средней точкой h, имеет вид (учитываем, что d = 2 * h):
S = d * fi = d * ( yi-1 + 4 * yi + yi+1 ) / 6
тогда формула для [a, b] с равномерной сеткой и четным числом шагов n (при длине шага d):
S = d * (ya + 4*f₁ + 2*f₂ + 4*f₃ + 2*f4 + ... + 2*fn-2 + 4*fn-1 + yb) / 6
ε ∝ Ο (1/n⁴) независимо от выбранного метода в процессе численного интегрирования необходимо оценить погрешность. погрешность уменьшается при увеличении количества разбиений. однако при этом возрастает погрешность округления
инженеры считают интеграл методом трапеций или прямоугольников. математики - симпсоном
--- математики - они оценивают точность приближения (после чего используют то, что лучше работает в конкретной ситуации)
--- нет, точность приближения оценивают физики (промежуточное звено между инженерами и математиками). математиков она вообще не интересует, а инженеры пользуются тем, что в справочнике написано
--- в общем-то верно: для инженерных задач (как правило) второго порядка точности считается достаточным, что трапеции и есть
есть N единичных интервалов по времени, N велико
есть зашумленный сигнал вида y(t) = const + помеха(t)
введем обозначения:
помеха в нечетный момент времени odd_fluct
помеха в четный момент времени even_fluct
нужно функцию проинтегрировать; чем болеe результат будет похож на const*N, тем лучше
интеграл методом трапеций ~ const*N + N/2 слагаемых вида (odd_fluct + even_fluct)
интеграл методом Симпсона ~ const*N + N/2 слагаемых вида (2/3 odd_fluct + 4/3 even_fluct)
теперь ключевое
предположим, что рапределение помех нормальное
предположим, что помехи не коррелированы между собой
обозначим дисперсию помехи как D
D(интеграл методом трапеций) = N*D
D(интеграл методом Симпсона) = 0.5N * (4/9+16/9) * D = 1.1*N*D
на 11% дисперсия возросла
а теперь пересчитайте эту дисперсию в предположении, что:
если подынтегральная функция имеет только первую "хорошую" производную - то трапеции симпсона в общей ситуации, разумеется, уделают. если и первую и вторую - симпсон уже начнёт уделывать трапеции (т.е. при одинаковой точности результата точек ему нужно будет меньше). что до остального - само собой, надо смотреть на конкретику ситуации (можно ли в этой конкретной таки прикинуть пресловутые вторые производные? или действительно нельзя?)
в физике никакой "данной точки" попросту вообще нету (измерение непрерывной величины всегда даёт не одно конкретное значение, а целый диапазон, причём попытка его сужения на микроуровень рано или поздно вообще даст заведомую фигню типа "температуры одной молекулы"), всё существенно упрощается: "числовая функция" оказывается на поверку чем-то вроде пары функций (в зазор между которыми мы укладываем вот те самые реально измеряемые прямоугольники "промежуток времени на диапазон значений параметра"), а вопрос о гладкости зависимости расшифровывается как "существует ли внутри намерянного зазора гладкая функция с требуемыми свойствами". погрешность в интеграле при этом складывается из двух частей:
1) принципиально неустранимая, отвечающая ИЗНАЧАЛЬНО существующему "физическому зазору"
2) математическая, связанная с ОБРАБОТКОЙ аппроксимирующей "физический зазор" математической функции по кубатурной формуле
привлекая вторые производные и т.д. мы просто предполагаем, что вот эта вышеприведенная вторая составляющая погрешности в нашей конкретной ситуации существенно больше, нежели первая составляющая - а таких ситуаций в приложениях действительно ведь куча
octave :
function s = matinteg (a, b, n, f)
A = tril (ones(n)) ;
t = linspace (a, b, n) ;
y = arrayfun (f, t) ;
S = (((b - a)/n) * A) * y' ;
S(n)
endoctave > matinteg (2,4,100, @ (x) x + 2) ans = 10