Численное интегрирование

Основная идея большинства методов численного интегрирования состоит в замене подынтегральной функции на более простую, интеграл от которой легко вычисляется аналитически. При этом для оценки значения интеграла получаются формулы вида I≈∑i=1nwif(xi),I \approx \sum\limits_{i = 1}^n {{w_i}} f\left( {{x_i}} \right),

где nn — число точек, в которых вычисляется значение подынтегральной функции. Точки xi{x_i}называются узлами метода, числа wi{w_i} — весами узлов. При замене подынтегральной функции на полином нулевой, первой и второй степени получаются соответственно методы прямоугольников, трапеций и парабол (Симпсона). Часто формулы для оценки значения интеграла называют квадратурными формулами.

Частным случаем является метод построения интегральных квадратурных формул для равномерных сеток, известный как формулы Котеса.  Основной идеей метода является замена подынтегральной функции каким-либо интерполяционным многочленом.  После взятия интеграла можно написать ∫abf(x)dx=∑i=0nHif(xi)+rn(f),\int\limits_a^b {f\left( x \right)} dx = \sum\limits_{i = 0}^n {{H_i}} f\left( {{x_i}} \right) + {r_n}\left( f \right),

где числа Hi{{H_i}} называются коэффициентами Котеса и вычисляются как интегралы от соответствующих многочленов, стоящих в исходном интерполяционном многочлене для подынтегральной функции при значении функции в узле xi=a+ih{x_i} = a + ih, где h=(b−a)÷n{h = \left( {b - a} \right) \div n} — шаг сетки, nn  — число узлов сетки, а индекс узлов i=0,1,…,ni = 0,1,…,n. Слагаемое rn(f){r_n}\left( f \right) — погрешность метода, которая может быть найдена разными способами. Для нечетных n⩾1n \geqslant 1 погрешность может быть найдена интегрированием погрешности интерполяционного полинома подынтегральной функции.

Частными случаями формул Котеса являются: формулы прямоугольников n=0n = 0, формулы трапеций n=1n =1, формула Симпсона n=2n =2.

Метод прямоугольников

Пусть требуется определить значение интеграла функции на отрезке [a,b]\left[ {a,b} \right]. Этот отрезок делится точками x0,x1,…,xn−1,xn{x_0},{x_1},…,{x_{n - 1}},{x_n} на nn равных отрезков длиной Δx=b−an\Delta x = \frac{{b - a}}{n}. Обозначим через y0,y1,…,yn−1,yn{y_0},{y_1},…,{y_{n - 1}},{y_n} значение функции f(x)f\left( x \right) в точках x0,x1,…,xn−1,xn{x_0},{x_1},…,{x_{n - 1}},{x_n}. Далее составляем суммы y0Δx+y1Δx+…+yn−1Δx{y_0}\Delta x + {y_1}\Delta x + … + {y_{n - 1}}\Delta x. Каждая из сумм — интегральная сумма для f(x)f\left( x \right) на [a,b]\left[ {a,b} \right] и поэтому приближённо выражает интеграл∫abf(x)dx=b−an(y0+y1+…+yn−1).\int\limits_a^b {f\left( x \right)} dx = \frac{{b - a}}{n}\left( {{y_0} + {y_1} + … + {y_{n - 1}}} \right).

Если заданная функция — положительная и возрастающая, то эта формула выражает площадь ступенчатой фигуры, составленной из «входящих» прямоугольников, также называемая формулой левых прямоугольников, а формула∫abf(x)dx=b−an(y1+…+yn−1+yn)\int\limits_a^b {f\left( x \right)} dx = \frac{{b - a}}{n}\left( {{y_1} + … + {y_{n - 1}} + {y_n}} \right)

выражает площадь ступенчатой фигуры, состоящей из «выходящих» прямоугольников, также называемая формулой правых прямоугольников. Чем меньше длина отрезков, на которые делится отрезок [a,b]\left[ {a,b} \right], тем точнее значение, вычисляемое по этой формуле, искомого интеграла.

Очевидно, стоит рассчитывать на бо́льшую точность если брать, в качестве опорной точки для нахождения высоты, точку посередине промежутка. В результате получаем формулу средних прямоугольников:∫abf(x)dx=h∑i=1nf(xi−1+h2)=h∑i=1nf(xi−h2),\int\limits_a^b {f\left( x \right)} dx = h\sum\limits_{i = 1}^n {f\left( {{x_{i - 1}} + \frac{h}{2}} \right)} = h\sum\limits_{i = 1}^n {f\left( {{x_i} - \frac{h}{2}} \right)} ,

где h=b−an.h = \frac{{b - a}}{n}.

Учитывая априорно бо́льшую точность последней формулы при том же объёме и характере вычислений её называют формулой прямоугольников.

Метод трапеций

Если функцию на каждом из частичных отрезков аппроксимировать прямой, проходящей через конечные значения, то получим метод трапеций.

Площадь трапеции на каждом отрезке: Ii≈f(xi−1)+f(xi)2(xi−xi−1){I_i} \approx \frac{{f({x_{i - 1}}) + f({x_i})}}{2}({x_i} - {x_{i - 1}})

Погрешность аппроксимации на каждом отрезке: ∣Ri∣⩽(b−a)312n2M2,i ,\left| {{R_i}} \right| \leqslant \frac{{{{\left( {b - a} \right)}^3}}}{{12{n^2}}}{M_{2,i}}{\mkern 1mu} ,

где M2,i=max⁡x2[xi−1,xi]∣f»(x)∣{M_{2,i}} = \mathop {\max }\limits_{x2[{x_{i - 1}},{x_i}]} \left| {f»(x)} \right|

Полная формула трапеций в случае деления всего промежутка интегрирования на отрезки одинаковой длины hh :I≈h(f(x0)+f(xn)2+∑i=1n−1f(xi)),I \approx h\left( {\frac{{f({x_0}) + f({x_n})}}{2} + \sum\limits_{i = 1}^{n - 1} f ({x_i})} \right),

где h=b−an.h = \frac{{b - a}}{n}.

Погрешность формулы трапеций:∣R∣⩽(b−a)312n2M2 ,\left| R \right| \leqslant \frac{{{{\left( {b - a} \right)}^3}}}{{12{n^2}}}{M_2}{\mkern 1mu} ,

где M2=max⁡x2[a,b]∣f»(x)∣{M_2} = \mathop {\max }\limits_{x2[a,b]} \left| {f»(x)} \right|.

Метод парабол (метод Симпсона)

Использовав три точки отрезка интегрирования, можно заменить подынтегральную функцию параболой. Обычно в качестве таких точек используют концы отрезка и его среднюю точку. В этом случае формула имеет очень простой видI≈b−a6(f(a)+4f(a+b2)+f(b)).I \approx \frac{{b - a}}{6}\left( {f(a) + 4f\left( {\frac{{a + b}}{2}} \right) + f(b)} \right).

Если разбить интервал интегрирования на 2N2N равных частей, то имеемI≈b−a6N(f0+4(f1+f3+…+f2N−1)+2(f2+f4+…+f2N−2)+f2N),I \approx \frac{{b - a}}{{6N}}\left( {{f_0} + 4\left( {{f_1} + {f_3} + \ldots + {f_{2N - 1}}} \right) + 2\left( {{f_2} + {f_4} + … + {f_{2N - 2}}} \right) + {f_{2N}}} \right),где fi=f(a+(b−a)i2N).{f_i} = f\left( {a + \frac{{(b - a)i}}{{2N}}} \right).

Когда первообразная не выражается в элементарных функциях, как для e−x2e^{-x^2}, или значения функции известны только из измерений, интеграл вычисляют приближённо. Отрезок [a,b][a, b] делят на nn равных частей длиной h=b−anh = \frac{b - a}{n} с узлами xi=a+ihx_i = a + ih. Формула трапеций: Tn=h2[f(x0)+2f(x1)+⋯+2f(xn−1)+f(xn)]T_n = \frac{h}{2}\left[f(x_0) + 2f(x_1) + \dots + 2f(x_{n-1}) + f(x_n)\right]. Формула Симпсона (парабол) для чётного числа частей nn: Sn=h3[f(x0)+4f(x1)+2f(x2)+⋯+4f(xn−1)+f(xn)]S_n = \frac{h}{3}\left[f(x_0) + 4f(x_1) + 2f(x_2) + \dots + 4f(x_{n-1}) + f(x_n)\right].

Погрешность формулы трапеций убывает как h2h^2, а формулы Симпсона — как h4h^4: если уменьшить hh вдвое, первая уменьшится примерно в 4 раза, а вторая примерно в 16 раз. Обеим формулам нужны только значения функции в узлах, поэтому они работают и для функции, заданной таблицей. Следите, чтобы в формуле Симпсона nn было чётным.

0,40,81,21,620,20,40,60,811,20xyy = 1/x
2
Формула трапеций
0,7083
Точно, ln 2
0,6931
Погрешность
0,0152
Зелёная ломаная соединяет точки графика y=1xy = \frac{1}{x} в узлах, а площадь под ней равна TnT_n. Если удвоить nn, погрешность уменьшится примерно в четыре раза; при n=4n = 4 получается T4T_4 из задачи для тренировки.

Гидрологи по глубинам, измеренным через каждые несколько метров, вычисляют площадь поперечного сечения реки; тот же приём нужен всегда, когда функция задана только таблицей измерений, например скорость из записи поездки.

Сечение реки

Ширина небольшой реки 40 м. Глубины, измеренные через каждые 10 м от берега до берега, равны 0; 2,1; 3,4; 2,6 и 0 м. Вычислите площадь поперечного сечения реки по формуле трапеций и по формуле Симпсона.

Решение

  1. Узлы отстоят друг от друга на h=10h = 10 м, поэтому n=4n = 4, а значения равны f0=0f_0 = 0, f1=2,1f_1 = 2{,}1, f2=3,4f_2 = 3{,}4, f3=2,6f_3 = 2{,}6 и f4=0f_4 = 0.
  2. Формула трапеций: T4=102[0+2(2,1+3,4+2,6)+0]=5⋅16,2=81T_4 = \frac{10}{2}\left[0 + 2(2{,}1 + 3{,}4 + 2{,}6) + 0\right] = 5 \cdot 16{,}2 = 81 м².
  3. Формула Симпсона, так как n=4n = 4 чётное: S4=103[0+4⋅2,1+2⋅3,4+4⋅2,6+0]=103⋅25,6≈85,33S_4 = \frac{10}{3}\left[0 + 4 \cdot 2{,}1 + 2 \cdot 3{,}4 + 4 \cdot 2{,}6 + 0\right] = \frac{10}{3} \cdot 25{,}6 \approx 85{,}33 м².
  4. Формула трапеций заменяет дно ломаной, которая при таком закруглённом дне проходит выше него и даёт меньшую площадь. Формула Симпсона повторяет дно параболами и обычно точнее.

Ответ: По формуле Симпсона площадь сечения реки около 85 м², а формула трапеций даёт 81 м².

Задача для тренировки

Вычислите ∫12dxx\int_1^2 \frac{dx}{x} по формуле трапеций и по формуле Симпсона при n=4n = 4 и сравните результаты с точным значением ln⁡2≈0,6931\ln 2 \approx 0{,}6931.

Показать решение

h=0,25h = 0{,}25, значения в узлах: 11; 0,80{,}8; 0,66670{,}6667; 0,57140{,}5714 и 0,50{,}5. T4=0,125⋅5,5762≈0,6970T_4 = 0{,}125 \cdot 5{,}5762 \approx 0{,}6970 и S4=0,253⋅8,3190≈0,6933S_4 = \frac{0{,}25}{3} \cdot 8{,}3190 \approx 0{,}6933, погрешности около 0,00390{,}0039 и 0,00010{,}0001.