Численное вычисление определенных интегралов

Численное вычисление определенных интегралов Постановка задачи. Классификация методов b Пусть требуется найти значение I интеграла  f ( x )dx для некоторой a заданной на отрезке a; b функции f(x). В математическом анализе обосновывается аналитический способ нахождения значения I с помощью формулы Ньютона-Лейбница I  F (b)  F (a) , где F(x) – первообразная для данной функции f(x). Эта формула не всегда применима [1]: - не существует первообразная среди элементарных функций для многих b b b 2 sin x dx функций f(x) (например, интегралы  ,  e  x dx ); dx ,  ln x a x a a - для многих реальных приложений определенного интеграла характерна дискретность задания подынтегральной функции f(x), например, функция представлена в виде таблицы (массива) значений в узлах некоторой расчётной сетки. Также возможно, что первообразная существует, но имеет слишком громоздкое выражение. В этих случаях применяются методы численного интегрирования. b Задача численного интегрирования заключается в вычислении  f ( x )dx a на основе значений подынтегральной функции y  f (x) с помощью специальных приближенных формул, называемыми квадратурными формулами (механическими квадратурами) или формулами численного интегрирования. Название «квадратурные формулы» связано с геометрическим смыслом определенного интеграла: вычисление b I   f ( x )dx a при f ( x )  0 равносильно построению квадрата, равновеликого криволинейной трапеции с основанием a; b и «крышей» f (x ) . В многомерном случае (размерность интеграла больше единицы) формулы для приближенного вычисления называют кубатурными. Методы численного интегрирования можно сгруппировать по способу аппроксимации подынтегральной функции [2]. 1. Методы Ньютона-Котеса основаны на полиномиальной аппроксимации функции f (x ) . Алгоритмы этого класса отличаются только степенью полинома. Узлы аппроксимирующего полинома, как правило, равноотстоящие. 1 2. Методы сплайн-интерполяции основаны на аппроксимации f (x ) сплайн-кусочным полиномом. Методы различаются по типу выбранных сплайнов. Такие методы используются в задачах, где алгоритмы сплайновой аппроксимации применяются для обработки данных. 3. В методах наивысшей алгебраической точности (метод Гаусса) используются специально выбранные неравноотстоящие узлы, обеспечивающие минимальную погрешность интегрирования при заданном количестве узлов. 4. Методы Монте-Карло используются чаще всего при вычислении кратных интегралов, и узлы выбираются специальным образом с помощью датчика случайных чисел. Ответ носит вероятностный характер. 5. В класс специальных группируются методы, алгоритмы которых разрабатываются на основе учёта особенностей конкретных подынтегральных функций, что позволяет существенно сократить время и уменьшить погрешность вычисления интегралов. Независимо от выбранного метода в процессе численного интегрирования необходимо вычислить приближенное значение интеграла и оценить погрешность. Погрешность уменьшается при увеличении n – количества разбиений отрезка a; b , но при этом возрастает погрешность округления за счёт суммирования значений интегралов, вычисленных на частичных отрезках (рисунок 1). Рисунок 1. Зависимость полной погрешности R от количества разбиений n интервала интегрирования Погрешность зависит от свойств подынтегральной функции и длины h частичного отрезка. Методы Ньютона-Котеса Пусть y  f (x) – непрерывная функция на отрезке a; b . Требуется вычислить b  f ( x )dx . a 2 Разобьём отрезок интегрирования на n равных частей точками a  x0  x1  x2   xn  b так, что ba xi 1  xi   h, i  0,n. n Пусть yi  f ( xi ) – значения f (x ) в точках деления. Приближенное значение интеграла ищут в виде линейной комбинации значений функции f(x) на отрезке a; b . Формулы прямоугольников На отрезке xi ; xi 1  возьмём произвольную точку ci . Интерполяционным многочленом в случае одного узла является f (ci ) . x i 1  f ( x)dx  h  f (ci ), h  xi 1  xi (1) xi Геометрическая интерпретация (1) – площадь, ограниченная графиком функции y  f (x) на отрезке xi ; xi 1  , осью x, ординатами f ( xi ) и f ( xi 1 ) , заменяется площадью прямоугольника с высотой, равной f (ci ) . Общий вид формулы прямоугольников: n 1 b f (ci ),  f ( x)dx  h  i 0 ci  xi ; xi 1 . (2) a Частные случаи (2) 1. Если ci совпадает с левым концом частичного промежутка интегрирования (рисунок 2), формулу называют формулой левых прямоугольников. ci  xi : b  n 1 f ( x )dx  I п   h  f ( xi ) . i 0 a Рисунок 2. Метод левых прямоугольников 3 2. Если ci совпадает с правым концом частичного промежутка интегрирования (рисунок 3), формулу называют формулой правых прямоугольников. ci  xi 1 : b  f ( x)dx  I a п n  h  f ( xi ) . i 1 Рисунок 3. Метод правых прямоугольников 3. Если ci середина частичного промежутка интегрирования (рисунок 4), формулу называют формулой средних прямоугольников. b n 1 xi  xi 1 h п (3) ci  :  f ( x )dx  I  h  f ( xi  ) . 2 2 i  a Рисунок 4. Метод средних прямоугольников Очевидно, что формулы левых и правых прямоугольников дают двухстороннее приближение к значению I интеграла от монотонной функции. Для возрастающей на a; b функции I п   I  I п  , для убывающей – I п  I  I п . 4 а) б) Рисунок 5. Геометрическое оценивание интеграла от монотонной функции с помощью а) I п  , б) I п  Погрешность составных прямоугольников [1]: методов левых, правых и средних ba f ()  h , 2 ba R п  (h)  f ()  h , 2 ba R п (h)  f ()  h 2 ,   (a; b) . (4) 24 Как видно из формул, при увеличении числа n элементарных отрезков, на которые разбивается промежуток интегрирования a; b , ошибка численного интегрирования по формуле средних прямоугольников убывает пропорционально квадрату шага h (метод второго порядка). Погрешность численного интегрирования непрерывной дифференцируемой функции по формулам левых и правых прямоугольников убывает по линейному закону. R п  (h)  Формула трапеций Выведем формулу трапеций, используя интерполяционный полином Ньютона (см. лекцию по интерполяции: y 2 y0 x  x0 x  x1   N n ( x )  y0  0 x  x0   1! h 2! h 2 3 y0 n y0 x  x0 x  x1 x  x2     n x  x0 x  x1 x  x2 x  xn 1  ,  3! h 3 n! h где n – степень полинома, (x0, y0), (x1, y1), (x2, y2), …, (xn, yn) – узлы интерполяции). xi ; xi 1  подынтегральную функцию Заменим на отрезке интерполяционным многочленом первой степени и проинтегрируем: 5 x i 1 x i 1 xi xi  f ( x)dx   yi yi  x  xi      y  x  x dx  y x  i i i   h h 2 2 xi 1  xi yi h 2 y h h  yi h   yi h  i   yi 1  yi  , h 2 2 2 yi  yi 1  yi Геометрическая интерпретация формулы трапеций (рисунок 6): на элементарном отрезке xi ; xi 1  площадь под кривой y  f (x) полагают равной площадью трапеции с основаниями yi, yi+1 и высотой h  xi 1  xi . Рисунок 6. Геометрическая интерпретация формулы трапеций Общая (составная) формула в случае линейной интерполяции называется формулой трапеции и имеет вид: b  y0  yn n 1  h (5) yi  .  f ( x)dx  2  y0  2 y1    yn 1   yn   h 2   i  1  a Погрешность1 составной формулы трапеций (5): b  a  f ()  h 2   (a; b) R( h )   , (6) 12 1 dx Пример 1. Вычислить интеграл I   [3]. Этот интеграл 2 1  x вычисляется по формуле Ньютона-Лейбница:  1 I  arctg x 0   0.785398 . 4 Используем теперь для вычисления данного интеграла формулы прямоугольников и трапеций. Разобьём отрезок интегрирования 0;1 на десять равных частей: n = 10, h = 0,1. Вычислим значения подынтегральной 1 Примечание. На первый взгляд, неожиданно, что метод трапеций имеет погрешность в два раза больше по абсолютной величине, чем метод средних прямоугольников, хотя аппроксимация подынтегральной функции проводилась полиномом первой, а не нулевой степени. То есть выбранный вариант аппроксимации подынтегральной функции прямой, проходящей через её значения на границах, не является оптимальным. Задача выбора аппроксимации полиномом заданной степенью с наименьшей возможной погрешностью была решена Гауссом, что привело к развитию целого класса методов [2]. 6 функции yi в точках разбиения xi и в полуцелых точках xi+1/2(серединах полученных отрезков) – таблица 1. Таблица 1. Значения подынтегральной функции и численное значение интеграла по формулам средних прямоугольников и трапеции к примеру 1 xi 0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 yi 1,000000 0,990099 0,961538 0,917431 0,862069 0,800000 0,735294 0,671141 0,609756 0,552486 0,500000 xi+1/2 0,05 0,15 0,25 0,35 0,45 0,55 0,65 0,75 0,85 0,95 yi+1/2 0,997506 0,977995 0,941176 0,890869 0,831601 0,767754 0,702988 0,640000 0,580552 0,525624 Iтр = 0,784981 0,000417 0,785398 Iп = 0,785606 ΔIп = 0,000208 ΔIтр = I= По формуле средних прямоугольников получим Iп = 0,784981, по формуле трапеций Iтр = 0,785606. Погрешность в вычислении интеграла составляет I п  I п  I  0.785606  0.785398  0.00021 (около 0,027 %). По формуле трапеций погрешность равна I тр  I тр  I  0,784981  0.785398  0.00042 (около 0,054 %). Таким образом, в рассмотренном примере лучшую точность вычисления даёт формула прямоугольников. Повышение точности здесь объясняется способом вычисления элементарных площадей, использующим значения в центральной точке на отрезке xi ; xi 1 . Использование формул левых и правых прямоугольников в этом примере приведёт к погрешности более 3 %. Формула Симпсона2 Если для каждой пары отрезков xi ; xi  2  построить многочлен второй степени, проинтегрировать его и воспользоваться свойством аддитивности интеграла, то получим формулу Симпсона. xi 2 xi 2   yi 2 yi  f ( x)dx    yi  h x  xi   2h 2 x  xi x  xi 1 dx ,  x x  i i 2 Симпсон Томас (Simpson Thomas); (1710–1761) – английский математик, ряд работ которого посвящён элементарной геометрии, тригонометрическому анализу и математическому анализу. В 1734 году Симпсон вывел формулу приближенного интегрирования (формула Симпсона). 7 xi 1  xi  h , xi 2  xi xi 2  xi yi yi  x  xi      y  x  x dx  y x  i i i   h h 2 2 2 yi x  xi x  xi 1 dx  2h 2 xi 2  xi x i 1 , xi 2 yi x  xi x  xi  h dx  2h 2 xi 2 3 2 2 yi   x  xi   x  xi   ,    h  2 2h 2  3  xi yi 2 yi  8h 3 4h 3  2  f ( x)dx  yi  2h  2h  4h  2h 2  3  2     xi h h  6 yi  6yi  2 yi  6 yi 1  yi  2 yi 1  yi  2   3 3 h   yi  4 yi 1  yi  2  . 3 h f ( x )dx   y0  4 y1  2 y2  4 y3  2 y4    2 yn  2  4 yn 1  yn   3 xi 2  b  a  n n   1 2 2   h   y0  yn  4 y2i 1  2  y2i  . (7) 3 i 1 i 1    Погрешность составной формулы Симпсона (7): b  a  f IV  ()  h 4   (a; b) R( h )   , . (8) 180 Из выражений остаточных членов (4), (6), (8) видно, что формулы средних прямоугольников и трапеций точны для многочленов первой степени, то есть для линейных функций, а формула Симпсона точна для многочленов третьей степени (для них остаточный член равен нулю). Погрешность формул средних прямоугольников и трапеций имеет второй порядок относительно h, а формула Симпсона является формулой четвёртого порядка точности. Формулы прямоугольников и трапеций в отдельности уступают при интегрировании гладких функций формуле Симпсона. Однако в паре они обладают следующим свойством: если f () не изменяют знака на a; b , то эти формулы дают двусторонние приближения для интеграла I, так как согласно (6) и (8) их остаточные члены имеют противоположные знаки [4]. На основании формул прямоугольников и трапеций можно получить уточнённые значения интегралов, если учесть характер погрешности этих формул [3]: I  2 I п  I тр 3 .   8 Для рассмотренного выше примера 1 получено Iп = 0,784981, Iтр = 0,785606. (с точность до I  2  0,785606  0,784981 3  0,785398 погрешностей округления), т.е. все шесть разрядов равны точным значениям. Метод Симпсона позволяет получить высокую точность, если f IV  ( x ) не слишком велика. В противном случае метод второго порядка может дать большую точность. Например, для функции f ( x)  25x 4  45x 2  7 формула трапеций при n  2 для 1  f ( x )dx  6 1 – даёт точный результат, тогда как по формуле Симпсона получается 1 2  f ( x )dx   3 . 1 Геометрическая интерпретация формулы Симпсона (рисунок 7). На отрезке xi ; xi  2  длиной 2h строится парабола, проходящая через точки xi , yi  , xi 1 , yi 1  , xi  2 , yi  2  . Площадь под параболой, заключённая между осью x и прямыми x  xi , x  xi  2 , принимают равной значению x i 2  f ( x )dx . xi Рисунок 7. Геометрическая интерпретация формулы Симпсона Для формулы Симпсона число разбиений отрезка интегрирования – чётное. 1 dx Пример. Вычислить по методу Симпсона интеграл I   . Значения 2 1  x функции при n = 10, h = 0,1 приведены в таблице 1. По формуле Симпсона: 9 0.1  y0  4 y1  y3  y5  y7  y9   2 y2  y4  y6  y8   y10   0.785398 , 3 что совпадает с точным значением (шесть значащих цифр). I Формулы Ньютона-Котеса высших порядков При замене подынтегральной функции интерполяционным многочленом третьей степени, получим квадратурную формулу Ньютона: x i 3 3h  f ( x)dx  8  yi  3 yi 1  3 yi  2  yi  3  , xi правило «трёх восьмых». На отрезке a; b (составная формула): b 3h  f ( x)dx  8  y0  yn   2 y3  y6    yn  3   a  3 y1  y2  y4  y5    yn  2  yn 1  Погрешность формулы равна: b  a  f IV  ()  h 4   (a; b) R( h )   , , 80 то есть при одинаковом шаге формула Ньютона менее точна, чем формула Симпсона. В смысле порядка точности квадратурные формулы с нечётным числом ординат являются более выигрышными. Фиксируя степень k интерполяционного многочлена равной 4, 5 и т.д., придём к частным формулам Ньютона-Котеса, подобным полученным простейшим формулам трапеций и Симпсона. Представим все эти частные случаи, включая уже рассмотренные, формулой вида xk k Ai yi ,  f ( x)dx   i 0 Ai  xk  x0 H i  hkHi , i  0, 1,k , x0 Hi – постоянные, называемые коэффициентами Котеса. Коэффициенты Котеса Коэффициенты Котеса для квадратурных формул со степенью интерполяционного полинома k представлены в таблице 2. Коэффициенты Котеса для каждого k представлены в виде дробей Hˆ Hi  i N с общим знаменателем N. Для контроля заметим, что k  Hˆ i  N . i 0 10 Таблица 2. Коэффициенты Котеса Степень полинома, k 1 2 3 4 5 6 Ĥ0 1 1 1 7 19 41 Ĥ1 1 4 3 32 75 216 Ĥ2 1 3 12 50 27 Ĥ3 Ĥ4 1 32 50 272 7 75 27 Ĥ5 19 216 Ĥ6 Общий знаменатель 41 N 2 6 8 90 288 840 При больших k (например, 8) коэффициенты Котеса могут быть отрицательными. Правило Рунге практического оценивания погрешностей На практике воспользоваться оценками квадратурных формул (4), (6), (8) – оценить максимальное значение производной – удаётся редко. Поэтому вычисление интегралов с нужной точностью обычно производят посредством последовательного дробления шага (как правило, делением пополам) до выполнения некоторых критериев точности. Обозначим через I h приближенное значение интеграла I, найденное по одной из трёх формул (3), (5), (7). Вычислив I 2 h , I h , I h 2 и убедившись, что выполняется неравенство [4]: Ih  Ih 2 2p   1  0.1 , I 2h  I h где p  2 для формул прямоугольников и трапеций, p  4 для формулы Симпсона, можно приблизительно оценить погрешность I  I h 2 по правилу Рунге: Ih 2  Ih I  Ih 2  p . 2 1 Кроме того, возможно найти уточнённое по Ричардсону значение интеграла I: 2 p Ih 2  Ih * Ih  2 p 1 с погрешностью I  I h*  O h p  2 . Формулы трапеций и Симпсона удобны тем, что при переходе от h к h/2 все вычисленные значения функции f могут быть использованы. Примеры вычисления уточнённых по Ричардсону значений интегралов для формул трапеций и Симпсона приведены в таблицах 3–4. Обозначим I2, I4 – приближенное значение для n = 2 и n = 4; I2,4 – уточнённое значение, I – точное значение интеграла. Тогда для формулы трапеций   11 22 I 4  I 2 I 2, 4   I 4  I 4  I 2  3 , 22  1 для формулы Симпсона 24 I 4  I 2 I 2, 4   I 4  I 4  I 2  15 . 24  1 Таблица 3. Экстраполирование для случая формулы трапеций I2 I4 I2,4 I 1,571 1,896 2,004 2,000 0,877 0,881 0,8823 0,8821 185,7090 179,5385 177,4817 177,4836 0,9695 0,9389 0,9287 0,9267  I   sin xdx 2 2 I   e  x dx 7 I   x 2 ln xdx 3 4 dx I  5 x 2 Таблица 4. Экстраполирование для случая формулы Симпсона I2 I4 I2,4 I 2,094 2,004 1,998 2,000 0,7833 0,7853 0,7854 0,7854 177,454 177,481 177,483 177,4836 0,0577 0,0541 0,0539 0,0533  I   sin xdx 1 dx 2 0 1 x I  7 I   x 2 ln xdx 3 4 I  dx 25  x  2 32 Примечание. Для применения метода Рунге необходимо знать, каков порядок точности исходной формулы. Фактический порядок точности p квадратурной формулы может отличаться от теоретического, например, в силу недостаточной гладкости функции. Определить фактический порядок точности p квадратурной формулы для заданной подынтегральной функции можно на основании процесса Эйткена [5]: p  log 2 I h  I 2h . Ih 2  Ih 12 Простейший алгоритм численного интегрирования, основанный на правиле Рунге 1. Задаём на отрезке интегрирования n узлов и вычисляем по выбранной квадратурной формуле приближенное значение In интеграла I по этим узлам. 2. Увеличиваем число узлов вдвое и вычисляем I2n. 3. По правилу Рунге находим оценку погрешности ошибки: I I   2 np n . 2 1 4. Если    , где  – заданная точность, то полагаем I  I 2n   и заканчиваем вычисления, иначе повторяем шаги 2–4. СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ 1. Вержбицкий, В. М. Основы численных методов: Учебник для вузов / В. М. Вержбицкий. – М.: Высш. шк., 2002. – 840 с. 2. Мудров А. Е. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль. – Томск: МП «Раско», 1991. – 272 с. 3. Турчак Л. И. Основы численных методов: Учеб. пособие. – М.: Наука. Гл. ред. физ.-мат. лит., 1987. – 320 с. 4. Волков Е. А. Численные методы: Учеб. пособие для вузов. – М.: Наука., 1987. – 248 с. 5. Калиткин Н. Н. Численные методы. – М.: Наука, 1978. – 512 с. 13