Блюмин А.Г., Федотов А.А., Храпов П.В. - Численные методы (1078410), страница 4
Текст из файла (страница 4)
метод ячеек имеет второй порядок точности относительно шагов сетки h1и h2 .29Глава 2. Приближенное вычисление двойного интегралаЗаметим, что поскольку в оценке (2.2.6) отброшены более высокие степениh1 и h2 , то соотношение для погрешности (2.2.7) является асимптотическим,т.е. выполняется при h1 → 0 и h2 → 0 с точностью до членов более высокогопорядка малости по h1 и h2 .Для вычисления интеграла (2.1.1) с заданной точностью можно воспользоваться, как это следует из (2.2.8), правилом Рунге практической оценкипогрешности. С помощью разложения f (x, y) в окрестности центра каждойячейки по формуле Тейлора до членов с производными четвертого порядкаможно получить оценку не только главного члена погрешности (2.2.8), но иоценить следующие по порядку малости h1 и h2 члены погрешности.
Заметим, что члены в разложении f (x, y) по формуле Тейлора, содержащие производные третьего порядка, в силу симметрии области интегрирования ∆Dijотносительно точки разложения не вносят вклад в погрешность интегрирования. Поэтому, для того чтобы учесть следующие после главного члена попорядку малости h1 и h2 члены погрешности, необходимо разлагать f (x, y) дочленов, содержащих производные четвертого порядка. В результате интеграл(2.1.1) можно представить в видеI = Ih + c1 · h21 + c2 · h22 + O(h41 + h21 · h22 + h42 ).(2.2.9)Выражения для c1 и c2 были приведены выше (см. формулу (2.2.7)).
Здесьважно подчеркнуть, что c1 и c2 — не зависящие от h1 и h2 постоянные величины, причем они не должны одновременно обращаться в 0.Таким образом, если известен главный член погрешности, то можно увеличить точность вычисления интеграла (2.1.1)I ≈ Ih + c1 · h21 + c2 · h22 .(2.2.10)Однако постоянные c1 и c2 являются неизвестными величинами. Для то30Глава 2. Приближенное вычисление двойного интегралаго, чтобы вычислить интеграл (2.1.1) с учетом главного члена погрешности,можно поступить следующим образом. Сначала находим значение Ih , затем— значение Ih/2 .
Здесь Ih/2 — значение интеграла (2.1.1), вычисленное поформуле (2.2.5) с шагами h1 /2 и h2 /2. Теперь, наряду с выражением (2.2.10),можно написать соотношениеI ≈ Ih/2 + c1 ·h122+ c2 ·h222.(2.2.11)Тот факт, что сетка по каждой переменной x и y сгущается в одинаковоечисло раз, позволяет в выражении (2.2.10) выделить главный член погрешности [c1 · (h1 /2)2 + c2 · (h2 /2)2 ] формулы (2.2.11):" 2 # 2h2h1+ c2 ·.I ≈ Ih + 4 c1 ·22(2.2.12)Из выражений (2.2.11) и (2.2.12) следует, что"" 2 2 2 # 2 #h1h1h2h2Ih + 4 c1 ·≈ Ih/2 + c1 ·.+ c2 ·+ c2 ·2222Из этого соотношения получаем выражение для главного члена погрешности формулы (2.2.11)c1 ·h122+ c2 ·h222≈Ih/2 − Ih.3(2.2.13)Теперь, согласно формуле (2.2.11), имеем приближенную оценку погрешности по правилу РунгеI − Ih/2 ≈Ih/2 − Ih.3(2.2.14)Наконец, подставляя выражение (2.2.13) в (2.2.11), получаем значение интеграла (2.1.1) с учетом главного члена погрешности, т.е.I ≈ Ih∗ =4 · Ih/2 − Ih,331Глава 2.
Приближенное вычисление двойного интегралагде Ih∗ — уточненное по Ричардсону значение интеграла I.Замечания. 1. Подчеркнем, что для практической оценки погрешности поправилу Рунге сетка по каждой переменной сгущается в одинаковое числораз, т.е. отношение m/n при сгущении сетки должно оставаться постоянным.В противном случае не удается в результате двойного пересчета интеграла(2.1.1) по двум сеткам с разными размерами ячеек составить формулы типа(2.2.10) и (2.2.11), из которых можно найти главный член погрешности.2. Понятно, что если одновременно c1 = 0 и c2 = 0, то для оценки погрешности вычисления интеграла (2.1.1) правило Рунге в виде (2.2.14) неприменимо.Мы получили формулу (2.2.5) для вычисления интеграла в простейшемслучае — для прямоугольной области.
Если область не прямоугольная, тов ряде случаев исходный интеграл по такой области соответствующей заменой переменных удобно преобразовать к двойному интегралу по прямоугольной области. Например, если область задана в виде криволинейногочетырехугольника D = {a 6 x 6 b, ϕ1 (x) 6 y 6 ϕ2 (x)} (рис. 2.2,a), тос помощью замены переменных x = x(u) = (b − a)u + a, y = ϕ1 (x(u))++ v · (ϕ2 (x(u)) − ϕ1 (x(u))) исходная область D преобразуется в квадратнуюобласть D0 = {0 6 u 6 1, 0 6 v 6 1} (рис.
2.2,б )Напомним правило замены переменных в двойном интеграле. Если ограниченная замкнутая область D в плоскости Oxy взаимно однозначно отображается на область D0 на плоскости Ouv с помощью непрерывно дифференцируемых функций x = x(u, v), y = y(u, v), причем якобиан преобразования ∂x ∂x ∂u ∂v 6= 0,J = ∂y∂y ∂u ∂v 32ôû é ôГлава 2.
Приближенноеìv вычисление двойного интегралаayávn2(x)1DD’n1(x)0ab0x1uРис. 2.2.то справедлива формулаZZZZf (x, y) dx dy =f (x(u, v), y(u, v)) | J | du dv.DD0Методом ячеек можно вычислить интеграл и по области сложной формы,например, с криволинейной границей (рис. 2.3).Интеграл в этом случае будем вычислять следующим образом. Нало-yжим на область D прямоугольную сетDку, и в интегральную сумму (2.2.5) будем включать только те ячейки, всеточки которых принадлежат областиD. В итоге на порядок понижаетсяточность формулы (2.2.5), поэтому длявычисления интеграла с достаточнойx0Рис.
2.3.точностью требуется сетка с более мелкими ячейками.33Глава 2. Приближенное вычисление двойного интегралаСледует отметить, что метод ячеек (2.2.5) легко переносится на большеечисло измерений (для вычисления тройных и большей кратности интегралов).
В случае однократного интеграла аналогом метода ячеек является метод средних прямоугольников (1.2.6), рассмотренный в гл. 1.2.2.2. Последовательное интегрирование с использованием формулы трапецийДругой метод вычисления двойных интегралов — их сведение к последовательному вычислению однократных интегралов. Снова рассмотрим интегралпо прямоугольной области D = {a 6 x 6 b, c 6 x 6 d} (рис. 2.4). Интеграл(2.1.1) можно вычислить последовательным интегрированиемZZZd ZbI=f (x, y) dx dy = f (x, y) dx dy.cD(2.2.15)aЭто выражение перепишем в видеZdI=ZbF (y) dy,F (y) =cf (x, y) dx.(2.2.16)aДля вычисления этих интегралов могут быть использованы известныеформулы из гл.
1. Например, пусть и по направлению x, и по направлениюy для приближенного вычисления применяется формула трапеций (1.2.7).ТогдаF (yj ) ≈ h1 ·mXq1,i · f (xi , yj ),i=0гдеq1,i 1/2 при i = 0 и i = m,= 1 при i = 1, 2, ..., m − 1.34(2.2.17)Глава 2. Приближенное вычисление двойного интегралаиI ≈ h2 ·nXq2,j · F (yj ),(2.2.18)j=0гдеq2,j 1/2 при j = 0 и j = n,= 1 при j = 1, 2, ..., n − 1.Подставляя выражение (2.2.17) в (2.2.18), получаем формулу последовательного интегрированияI ≈ Ih = h1 · h2 ·m XnXqij · f (xi , yj ),(2.2.19)i=0 j=0гдеqij = q1,i · q2,j1/4 при i = 0, m; j = 0, n; i = 0 и i = m, j = 1, ..., n − 1,1/2 при= j = 0 и j = n, i = 1, ..., m − 1; 1 при i = 1, ..., m − 1, j = 1, ..., n − 1.На рис.
2.4 приведена сетка, которая используется при приближенномвычислении интеграла (2.1.1) по формуле (2.2.19). Точками, кружочкамии квадратиками показаны узловые точки, в которых коэффициенты qij =1, 1/2 и 1/4, соответственно.Легко убедиться в том, что для дважды непрерывно дифференцируемойфункции f (x, y) формула (2.2.19) имеет второй порядок точности относительно шагов h1 и h2 и что можно применить правило Рунге практическойоценки погрешности.Замечание. Если в методе последовательного интегрирования воспользуемся формулой средних прямоугольников при интегрировании по каждомуиз направлений x и y, то в результате получим расчетную формулу методаячеек (2.2.5).35Глава 2.
Приближенное вычисление двойного интегралаyh1yn=dyjh2yj -1y0=c0x0=axi -1xixm=bxРис. 2.4.Случай сложной области. Метод последовательного интегрирования можно непосредственно применять и к области произвольной формы, например скриволинейной границей (см. рис. 2.3). Однако для получения простых расчетных формул на практике всегда стараются свести исходный интеграл ксумме интегралов по прямоугольным областям.2.2.3. Последовательное интегрирование с использованием квадратурных формул ГауссаДля получения квадратурной формулы более высокой точности, чем формулы (2.2.19) можно воспользоваться квадратурными формулами Гаусса.При этом предварительно заменой переменных x(u) = (a + b)/2 + u · (b − a)/2,y(v) = (c + d)/2 + v · (d − c)/2 прямоугольная область {a 6 x 6 b, c 6 x 6 d}преобразуется в квадратную область D = {−1 6 u 6 1, −1 6 v 6 1}.
Поэтому будем считать, что с самого начала требуется вычислить интеграл по36Глава 2. Приближенное вычисление двойного интегралаобласти D = {−1 6 x 6 1, −1 6 y 6 1}:Z1 Z1ZZf (x, y) dx dy =I=f (x, y) dxdy.(2.2.20)−1 −1DПрименяя для интегрирования (2.2.20) и по направлению x, и по направлению y квадратурную формулу Гаусса с одинаковым числом узлов, получаем следующую формулу последовательного интегрированияIn =n XnXqi · qj · f (xi , yj ).(2.2.21)i=1 j=1Значения координат узловых точек и весовых коэффициентов по направлениям x и y берутся из табл. 1.1.
Расположение узловых точек для n = 3 иn = 4 проиллюстрировано на рис. 2.5,a и рис. 2.5,б, соответственно:y1-1ya11 x-11 x-1-1Рис. 2.5.37ɛГлава 2. Приближенное вычисление двойного интеграла2.3. ЗаданиеДля предложенного варианта лабораторной работы вычислите двойнойинтеграл по области D, где D — криволинейный четырехугольник {a 6 x 6 b,ϕ1 (x) 6 y 6 ϕ2 (x)}:ZZf (x, y) dx dy.I=DИнтегралы вычислите:1) аналитически,2) численно с точностью до ε = 0.0001:• методом ячеек,• последовательным интегрированием с использованием формулытрапеций для интегрирования по направлениям x и y.При численном решении область D предварительно отобразите в квадратD0 (см. рис. 2.2).