Кратные интегралы.

ЛЕКЦИЯ №6. КРАТНЫЕ ИНТЕГРАЛЫ 6.1. Метод ячеек Приближенная оценка двукратного интеграла по прямоугольной области может быть дана следующим образом d b   f ( x, y)dxdy  S f ( x , y ) , (6.1) c a где S  (b  a )( d  c), x  ab cd , y 2 2 - координаты центра прямоугольника площади S. Для повышения точности разобьем область интегрирования на прямоугольные ячейки системой линий, параллельных осям координат, и построим обобщенную квадратурную формулу d b   f ( x , y )dxdy   S где i f ( xi , y i ) , (6.2) i c a x i , y i - координаты центра прямоугольной ячейки. Оценим погрешность интегрирования. Пусть в обобщенной квадратурной формуле (6.2) стороны прямоугольника разбиты соответственно на N и M равных частей. Тогда погрешность этой формулы для единичной ячейки Ri  1 ba 2 d c 2 S i [( ) f xx  ( ) f yy ] 24 N M . Суммируя погрешность по всем ячейкам, получим суммарную погрешность R  1 ba 2 d c 2 [( )  f xx dxdy  ( )  f yy dxdy]  O( N 2  M 2 ) . 24 N M G G т.е. формула имеет второй порядок точности. Если граница области интегрирования G криволинейная, то формулу (6.2) применяют, накладывая на область G прямоугольную сетку. Все ячейки разделяют на внутренние и граничные. Площадь внутренней ячейки равна произведению ее сторон. Площадь граничной ячейки вычисляют приближенно, заменяя в пределах данной ячейки истинную границу области хордой. Значения этих площадей подставляют в (6.2). Погрешность формулы (6.2) при этом будет такой. В каждой внутренней ячейке O( N 2 ) , в каждой граничной ячейке относительная ошибка есть ошибка составляет O( N 1 ) , так как центр тяжести прямоугольной ячейки не совпадает с центром тяжести входящей в интеграл части. Но граничных ячеек примерно в N раз меньше чем внутренних. Поэтому общая погрешность будет O( N 2 ) , т.е. имеем второй порядок точности. Можно граничные ячейки вообще не включать в сумму. Погрешность при этом будет O( N 1 ) . 6.2. Последовательное интегрирование Рассмотрим интеграл по прямоугольной области d b b c a a I    f ( x, y) dx dy   F ( x) dx , (6.3) где d F ( x)   f ( x, y)dy . c По каждой координате введем сетку узлов. Каждый однократный интеграл вычисляют на данной сетке по квадратурным формулам: I   Ai F ( x i ) , где F ( xi )   Aj f ( xi , y j ) . i j Тогда I   Ai Aj f ( xi , y j )   g i j f ( xi , y j ) . i j i j Для разных направлений можно использовать квадратурные формулы разных порядков точности (трапеций, прямоугольников, средних, Симпсона и т.д.). Можно также использовать формулы Гаусса, тогда gi j  xi  1 (b  a )( d  c) i  j , 4 ab ba  i , 2 2 yj  cd d c  j , 2 2 1  i, j  n , где i , j ,  i ,  j - нули многочленов Лежандра и веса формулы Гаусса. Теперь пусть область интегрирования ограничена непрерывными однозначными кривыми ( ( x)   ( x)) y   ( x), y   ( x) и двумя вертикалями x  a и x  b . Имеем b  ( x) b a  ( x) a I   f ( x, y )dxdy   dx G  f ( x, y)dy   F ( x)dx , где  ( x) F ( x)   f ( x, y)dy  . ( x) Отсюда n  ( xi ) m i 1  ( xi ) j 1 I   Ai F ( xi ) , F ( xi )  n  f ( xi , y )dy   Bi j f ( xi , y j ) , m I   f ( x, y )dxdy   Ai Bi j f ( xi , y j ) , G i 1 j 1 где Ai ,Bi j - известные постоянные. Пример. Получить кубатурную формулу для вычисления двукратного интеграла d b   f ( x, y)dxdy , используя по каждому направлению формулу Симпсона. c a Пусть областью интегрирования является прямоугольник со сторонами (b-a) и (dc). Разобьем каждую сторону данного прямоугольника на N и M интервалов, т.е. введем на каждой границе области сетку с N+1 и M+1 узлами. Соответствующие шаги сетки hx  ba d c , hy  . N M Согласно (6.3) получим d b  b d b a c a f ( x, y )dxdy   dx  f ( x, y )dy   F ( x) dx  c a h  x 3 N 1 2  F i 0 2i (6.4)  4 F2i1  F2i2 , где каждое слагаемое под знаком суммы d Fk   f ( xk , y ) dy  c hy 3 M 1 2  f (x , y k j 0 2j  )  4 f ( xk , y2 j 1 )  f ( xk , y2 j 2 ) , (6.5) здесь k  2i, 2i  1, 2i  2 . Чтобы получить представление о том, как выглядит формула численного вычисления двукратного интеграла в развернутом виде, рассмотрим частный случай, когда N=M=2. В данном случае имеем 9 узлов с координатами ( x0 , y0 ), ( x0 , y1 ), ( x0 , y2 ), ( x1, y0 ), ( x1, y1 ), ( x1, y2 ), ( x2 , y0 ), ( x2 , y1 ), ( x2 , y2 ), Тогда с учетом (6.4) и (6.5) кубатурной формуле можно придать вид d b hx   f ( x, y)dxdy  3 F  4 F1  F2   hx h y c a 9 [ f ( x0 , y 0 )  f ( x 2 , y 0 )  f ( x0 , y 2 )   f ( x 2 , y 2 )  4[ f ( x1 , y 0 )  f ( x0 , y1 )  f ( x 2 , y1 )  f ( x1 , y 2 )]  16 f ( x1 , y1 )]. При выводе принято во внимание, что d F0   f ( x0 , y ) dy  c d F1   f ( x1 , y ) dy  c d F2   f ( x2 , y ) dy  c hy 3 hy 3 hy 3 [ f ( x0 , y0 )  4 f ( x0 , y1 )  f ( x0 , y2 )], [ f ( x1 , y0 )  4 f ( x1 , y1 )  f ( x1 , y2 )], [ f ( x2 , y0 )  4 f ( x2 , y1 )  f ( x2 , y2 )] (6.6) При построении алгоритма программы конструировать общие формулы типа (6.6) нет нужды. Рабочими формулами являются (6.4), (6.5), т.е. алгоритм расчета двукратного интеграла строится согласно (6.4), где каждое слагаемое определяется из (6.5), как это и показано в предыдущем примере. Если область интегрирования криволинейная, то в самом простом варианте можно построить прямоугольник R , стороны которого параллельны осям координат и ввести вспомогательную функцию ( x, y )  G  f ( x, y ), f * ( x, y )   . ( x, y )  R  G  0, В таком случае, очевидно  f ( x , y )dxdy   f G * ( x , y )dxdy . G Последний интеграл может быть вычислен по общей кубатурной формуле. Подчеркнем еще раз, что вычисление интегралов по каждому направлению x и y выполняется независимо, т.е. по разным направлениям могут быть применены различные методы численного интегрирования, например, по одному направлению метод Симпсона, а по другому метод Гаусса. Аналогично формулами алгоритму вычисления двукратного интеграла, задаваемого (6.4), (6.5), строятся алгоритмы численных расчетов интегралов более высокой кратности.