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

1 ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ 14.1. Постановка задачи Если функция f(х) непрерывна на отрезке [a,b] и известна ее первообразная F(x), то определенный интеграл от этой функции в пределах от a до b может быть вычислен по формуле Ньютона-Лейбница b  f(x)dx = F(b) F(a) a где F' (х) = f(x). Но часто возникают ситуации, когда вычислить интеграл можно только с помощью численных методов: 1) F(x) не выражается через элементарные функции. b  x dx; е 2 a 2) F(x) существует и выражается через элементарные функции, но ее сложно найти b 30 29  x + x +...+ 1 dx;  40 39  x + x +...+ 1 a 3) найдена F(x), но сложно вычислить ее значение; 4) f(х) задана таблично или графиком. b f(x)dx. Итак, как вычислить  a Обычный прием состоит в том, что данную функцию f(х) на рассматриваемом отрезке [a,b] заменяют интерполирующей функцией Pn(x) простого вида, а затем приближенно полагают: b b a a  f(x)dx   Pn(x)dx . b Функция Pn(x) должна быть такова, чтобы интеграл  Pn(x)dx выa числялся непосредственно. Можно использовать интерполяционный многочлен Pn(x) различной степени n, n =  14.2. Формулы прямоугольников При n=f(x)  P0(x) 2 Для построения Р0(х) требуется одна точка (х0, f(х0) ). P0(x) = f(x0) Формула левых прямоугольников: ( а, f(а) ) b b b b a a a a  P(x)dx =  f(x)dx   f(а)dx = f(а) x | = f(а) (b  a) Геометрическая иллюстрация Формула правых прямоугольников: ( b, f(b) ) b b b b a a a a P(x)dx =   f(x)dx   f(b)dx = f(b) x | = f(b) (b  a) Геометрическая иллюстрация Формула центральных прямоугольников: ( a+b a+b ;f( ) ) 2 2 b b b b a+b a+b a+b   ) dx = f ( )x |=f( ) (b  a) P(x)dx = f ( f(x)dx   2 2 2  a a a a Геометрическая иллюстрация 3 14.3. Формула трапеций При n= f x   P1 x  Для построения Р1(х) требуется две точки: х у y0  f a  x0  a y1  f b x1  b   1 xx  j  P1 x     yi    i 0  j 0 xi  x j  j i   1 Р1(x) = y0 xb xa 1 + f(b) = (f(b)(xa) f(a)(xb)) ab ba ba Р1(x) = f(a) b xx1 xx0 + y1 x0x1 x1x0 b b 1   P(x)dx = (f(b)(xa) f(a)(xb)) dx = f(x)dx   ba a a a = b 1  f(b) (xa) dx  a ba  = b f(a) a   2 b 2 b   f(b) (xa)  f(a) (xb)   =  2 a 2 a    (ab)2  f(a) + f(b) f(a) (ba) = 2 2  (xb) dx  = 1  (ba)2  f(b) 2 ba  1 ba Геометрическая иллюстрация 4 14.4. Формула Симпсона При n= f x   P2 x  Для построения Р2(х) требуется три точки: х x0  a x1  ab 2 x2  b у y0  f a   ab y1  f    2  y2  f b  2 xx  j P2 x     yi   i 0  j 0 xi  x j j i  2      a+b a+b )(xb) (xa) (x ) 2 2 a+b (xa) (xb) Р2(x) = f(a) + f( ) + f(b) a+b 2 a+b a+b a+b (a ) (ab) ( a) ( b) (ba) (b ) 2 2 2 2 ( х b b ba a+b f(x)dx  P(x)dx = f(a) + f(b) + 4 f( )) (   6 2 a a 5 Геометрическая иллюстрация 14.5. Обобщенные формулы На практике обычно пользуются обобщенными формулами, т.к. [a,b] может быть большим и, следовательно, большой и погрешность вычисления интеграла по формулам прямоугольников, трапеции и Симпсона. Как же добиться повышения точности вычисления? [a,b] разбивают на n равных частей точками a=x0x1 ...  xn=b, и на каждом отрезке [xi, xi+1] применяется конкретный метод прямоугольников, трапеции или Симпсона, результаты суммируются, пользуясь условием аддитивности определения интеграла. ba Величина h = - шаг интегрирования, xi=x0+ih , где i= 0, n . n Обобщенная формула левых прямоугольников Геометрическая иллюстрация b n1 f(x)dx   h f(xi)  a i=0 Введем обозначение f(xi) = fi: b n1 f(x)dx  h  fi  a i=0 (14.1) 6 Обобщенная формула правых прямоугольников Геометрическая иллюстрация b n f(x)dx   h f(xi)  a i=1 Введем обозначение f(xi) = fi: b n f(x)dx  h  fi  a i=1 (14.2) Обобщенная формула центральных прямоугольников Геометрическая иллюстрация b a Введем обозначение f( b n1 f(x)dx   h f(  i=0 xi + xi+1 ) 2 xi + xi+1 ) = fi + 12 : 2 n1 f(x)dx  h  fi + a 1 2 i=0 Обобщенная формула трапеции (14.3) 7 Геометрическая иллюстрация b  n1  n1 f + f n1 f f f0 + fn  i i+1  i + i+1=h f0 + f1 + f1 + f2 +...+ fn-1 + fn =h f(x)dx  h =h f +  i    2 2   2 2 2 2  2 2 2   2  a   i=1 i=0 i=0  n1 f(a) + f(b) f(x)dx  h   fi +  2 a  i=1 b     (14.4) Обобщенная формула Симпсона b n1 h f(x)dx    fi + fi+1 + 4fi + 6  a i=0 b a f(x)dx  h n1  f + f + 4fi + 6   i i+1 i=0 1 2    1 2    (14.5) 14.6. Оценка погрешности b J = f(x)dx  a Если подынтегральная функция f имеет на отрезке [a,b] непрерывную производную f , то оценка погрешностей формул (14.1) и (14.2) дается неравенством 2  b  a J  J n  M1  , (14.6) 2n где M 1  max f ( x ) . a ,b 8 Если подынтегральная функция f имеет на отрезке [a,b] непрерывную вторую производную f , то оценка погрешностей формулы (14.3) дается неравенством 3  b  a J  Jn  M2  , (14.7) 24n 2 f ( x ) . где M 2  max a ,b  Если подынтегральная функция f имеет на отрезке [a,b] непрерывную вторую производную f , то оценка погрешностей формулы (14.4) дается неравенством 3  b  a J  Jn  M2  , (14.8) 12n 2 f ( x ) . где M 2  max a ,b Если подынтегральная функция f имеет на отрезке [a,b] непрерывную четвертую производную f IV, то оценка погрешностей формулы (14.5) дается неравенством J  Jn  M4 5  b  a  180n 4 , (14.9) f IV ( x ) . где M 4  max a ,b  Эмпирический критерий оценки точности вычисления интеграла На практике широко применяется следующий прием, пригодный для каждого из рассматриваемых методов. Искомый интеграл вычисляется дважды: при делении отрезка [a,b] на n частей и на 2n частей. Полученные интегралы Jn J2n сравниваются, и совпадающие первые десятичные знаки считаются верными. 14.7. Квадратурные формулы Гаусса При переходе от формулы трапеций к формуле Симпсона точность вычисления интеграла при фиксированном значении n значительно возрастает, так как возрастает порядок аппроксимации подынтегральной функции f(x): линейная аппроксимация сменяется на параболическую. При этом метод трапеций приводит к точному значению, если f(x) линейная функция на интервале интегрирования, а метод Симпсона – если f(x) квадратичная функция. Разбиение отрезка интегрирования может быть произвольным, хотя на практике чаще предпочитают использовать равноотстоящее. 9 Существует иной метод построения квадратурных формул, в котором центральное место играет выбор узлов интерполирования подынтегральной функции. Традиционно при получении квадратурных формул Гаусса в исходном интеграле выполняется замена переменных, переводящая интеграл по отрезку a, b в интеграл по отрезку  1,1: a  1 2 x  (b  a) t , при этом , b  1 ba (14.10) 1 1 x  (b  a)t  (b  a) . 2 2 Тогда b 1 1 1 1 J   f ( x)dx  (b  a)  f (b  a)t  (b  a) dt . (14.11) 2 2 2 a 1 Далее, не теряя общности, будем рассматривать метод Гаусса при-   1 менительно к интегралу вида  (t )dt . 1 y y A2 S3 S2 S A1 S1 –1 a) 1 t –1 t1 t2 1 t б) На рис. a) в качестве узлов линейной интерполяции выбраны концы отрезка  1,1, что соответствует применению малой формулы трапеций. Различие в площадях S криволинейной трапеции, ограниченной сверху кривой y  (x) , и обычной трапеции, ограниченной сверху прямой, проведенной через концы графика функции на отрезке интегрирования, фиксировано видом заданной подынтегральной функции и определяет погрешность подсчета значения определенного интеграла. Сделаем узлы интерполяции подвижными. Тогда их можно выбрать таким образом, чтобы разность между площадями криволинейной и обычной трапеций была минимальна. В идеале: S1+S3=S2, тогда интеграл аппроксимирован точно. Но нахождение значений t1 и t2 настолько же трудо- 10 емкая задача, как и вычисление точного значения определенного интеграла. Сформулируем задачу следующим образом: подобрать t1 и t2 таким образом, чтобы площадь трапеции, ограниченной сверху прямой, проходящей через точки A1 t1 , (t1 ) и A2 t2 , (t2 ) , была равна точному значению определенного интеграла от многочлена наивысшей возможной степени. Так как положение точек A1 и A2 определено четырьмя координатами, то и этот многочлен может задаваться четырьмя коэффициентами, то есть являться многочленом третьей степени: P3 (t )  a0  a1t  a2t 2  a3t 3 . (14.12) Уравнение прямой (A1, A2) имеет вид t  t1   1  , t2  t1 2  1 (14.13) 2  1 t2 1  t12 y t , t2  t1 t2  t1 где 1  (t1 ) , 2  (t2 ) . Таким образом, необходимо отыскать значения t1 и t2 так, чтобы выполнялось равенство 1  2  1 t2 1  t12  2 3 1 t  t t  t  t dt  1 a0  a1t  a2t  a3t dt . (14.14)  2 1  2 1 1 Найдем интегралы в левой и правой части равенства: 2  1 t 2 t 2  t1 2 1 1 t2 t 2 1  t12 1  t  a0t 1  a1 2 t  t 2 1 1 1 0 0 1 t3  a2 3 1 1 t4  a3 4 1 1 , 1 0 t21  t12 2  2a0  a2 , t2  t1 3 t21  t12 1  a0  a2 . t2  t1 3 Подставим значения 1  P3 (t1 ) , 2  P3 (t2 ) : t2 a0  a1t1  a2t12  a3t13   t1 a0  a1t2  a2t22  a3t23   a0  1 a2 . t2  t1 t2  t1 3 Сгруппируем в левой части слагаемые относительно коэффициентов многочлена: 1 a0  0  a1  t1t2  a2  t1t2 t1  t2   a3  a0  a2 . 3 Данное равенство выполняется тогда и только тогда, когда 1 t1t2   , t1t2 t1  t2   0 . 3 2 11 1 1 , t2  . 3 3 Вычислим, учитывая эти значения, площадь трапеции, ограниченной сверху прямой (A1, A2): 1 1 1  2 1  2  1 t 2 1  t12  t 2 1  t12 3 3 1 t  t t  t  t dt  2 t  t  2 1 1  1  2   2 1  2 1 2 1  3 3  1   1        . 3  3  Формула 1  1   1  (14.15) 1 (t )dt    3    3  называется квадратурной формулой Гаусса. Данная формула дает точные значения для подынтегральной функции, являющейся многочленом не выше третьей степени. Для иных функций формула Гаусса даст лишь приближенное значение, но из геометрической интерпретации ясно, что точность его будет выше, нежели значение полученной с помощью формулы трапеций. Следовательно t1   b Применительно к исходному интегралу  f ( x)dx формула Гаусса a примет вид: J ba ba ba  b  a b  a     f  f  . 2   2 2 3 2 3   2 (14.16) Для повышения точности вычисления интеграла разобьем отрезок a, b на n равных частей длины h  b  a . Получим обобщенную формулу n Гаусса h n1   h h  h h   (14.17) J n    f  xi     f  xi    . 2 i 0   2 2 3 2 2 3   Если подынтегральная функция f имеет на отрезке [a,b] непрерывную четвертую производную f IV, то оценка погрешностей формулы (14.17) дается неравенством 5  b  a , (14.18) J  Jn  M 4  4320n 4 f IV ( x ) . где M 4  max a ,b 