Вычисление кратных интегралов

Министерство образования Украины

Днепропетровский государственный университет

–––––––––––––––––––––––––––––––––––––––––––––

Факультет прикладной математики

Кафедра вычислительной механики и прочности конструкций

КУРСОВАЯ РАБОТА

по численным методам в механике

на тему

Вычисление кратных интегралов

методом ячеек

с автоматическим выбором шага

Исполнитель: студент группы ПД-97-1 Коваленко А.В.

Руководитель: профессор Мусияка В.Г.

Днепропетровск 1999

Содержание

1 Постановка задачи 2

2 Теоретическая часть 2

2.1 Понятие о кубатурных формулах 2

2.2 Метод ячеек 3

2.3 Последовательное интегрирование 5

2.4 Кубатурная формула типа Симпсона 6

2.5 Принципы построения программ с автоматическим выбором шага 8

3 Список использованной литературы 9

4 Практическая часть 9

4.1 Решение задачи 9

4.2 Блок-схема программы 10

4.3 Листинг программы 12

4.4 Результаты решения 13

1 Постановка задачи

Найти при помощи метода ячеек значение интеграла , где – область, ограниченная функциями .

2 Теоретическая часть

Рассмотрим K-мерный интеграл вида:

(1)

где - некоторая K-мерная точка. Далее для простоты все рисунки будут сделаны для случая K=2.

2.1 Понятие о кубатурных формулах

Кубатурные формулы или, иначе формулы численных кубатур предназначены для численного вычисления кратных интегралов.

Пусть функция определена и непрерывна в некоторой ограниченной области . В этой области выбирается система точек (узлов) . Для вычисления интеграла приближённо полагают:

(2)

Чтобы найти коэффициенты , потребуем точного выполнения кубатурной формулы (2) для всех полиномов

(3)

степень которых не превышает заданного числа . Для этого необходимо и достаточно, чтобы формула (2) была точной для произведения степеней . Полагая в (1) , будем иметь:

(4)

Таким образом, коэффициенты формулы (2), вообще говоря, могут быть определены из системы линейных уравнений (4).

Для того чтобы система (4) была определённой, необходимо, чтобы число неизвестных было равно числу уравнений. В случае получаем:

2.2 Метод ячеек

Рассмотрим K-мерный интеграл по пространственному параллелепипеду . По аналогии с формулой средних можно приближённо заменить функцию на её значение в центральной точке параллелепипеда. Тогда интеграл легко вычисляется:

(5)

Для повышения точности можно разбить область на прямоугольные ячейки (рис. 2). Приближённо вы­числяя интеграл в каждой ячейке по формуле средних и обозначая через соответственно пло­щадь ячейки и координаты её центра, получим:

(6)

Справа стоит интегральная сумма; следовательно, для любой непрерывной она сходится к зна­чению интеграла, когда периметры всех ячеек стремятся к нулю.

Оценим погрешность интегрирования. Формула (5) по самому её выводу точна для . Но непосредственной подстановкой легко убедиться, что формула точна и для любой линейной функции. В са­мом деле, разложим функцию по формуле Тейлора:

(7)

где , а все производные берутся в центре ячейки. Подставляя это разложение в правую и левую части квадратурной формулы (5) и сравнивая их, аналогично одномерному случаю легко получим выражение погрешности этой формулы:

(8)

ибо все члены разложения, нечётные относительно центра симметрии ячейки, взаимно уничтожаются.

Пусть в обобщённой квадратурной формуле (6) стороны пространственного параллелепипеда разбиты соответст­венно на N1, N2, …, Nk равных частей. Тогда погрешность интегрирования (8) для единичной ячейки равна:

Суммируя это выражение по всем ячейкам, получим погрешность обобщённой формулы:

(9)

т.е. формула имеет второй порядок точности. При этом, как и для одного измерения, можно применять метод Рунге–Ромберга, но при одном дополнительном ограничении: сетки по каждой переменной сгу­щаются в одинаковое число раз.

Обобщим формулу ячеек на более сложные области. Рассмотрим случай K=2. Легко сообразить, что для линейной функции формула типа (5) будет точна в области произвольной формы, если под S подразуме­вать площадь области, а под ­–координаты центра тяжести, вычисляемые по обычным формулам:

(10)

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

Для области с произвольной границей формулу (6) применяют иным способом. Наложим на область сетку из K-мерных параллелепипедов (рис.3). Те ячейки сетки, все точки которых принадлежат области, на­зовём внутренними; если часть точек ячейки принадлежит области, а часть – нет, то назовём ячейку гра­ничной. Объём внутренней ячейки равен произведению её сторон. Объёмом граничной ячейки будем считать объем той её части, которая попадает внутрь ; этот объём вычислим приближённо. Эти площади подставим в (6) и вычислим интеграл.

Оценим погрешность формулы (6). В каждой внутренней ячейке ошибка составляет по отношению к значению интеграла по данной ячейке. В каждой граничной ячейке относительная ошибка есть , ибо центр ячейки не совпадает с центром тяжести входящей в интеграл части. Но самих граничных ячеек примерно в раз меньше, чем внутренних. Поэтому при суммировании по ячейкам общая погрешность будет , если функция дважды непрерывно диф­ференцируема; это означает второй порядок точности.

Вычисление объёма граничной ячейки довольно трудоёмко, ибо требует определения поло­жения границы внутри ячейки. Можно вычислять интегралы по граничным ячейкам более грубо или во­обще не включать их в сумму (6). Погрешность при этом будет , и для хорошей точности потре­буется более подробная сетка.

Мы видели, что к области произволь­ной формы метод ячеек трудно применять; поэтому всегда желательно заменой переменных преобразовать об­ласть интегрирования в прямоугольный параллелепипед (это относится практически ко всем методам вычисления кратных интегралов).

2.3 Последовательное интегрирование

Снова рассмотрим интеграл по K-мерной области, разбитой сеткой на ячейки (рис. 2). Его можно вычислить последовательным интегрированием:

Каждый однократный интеграл легко вычисляется на данной сетке по квадратурным формулам типа:

Последовательное интегрирование по всем направлениям приводит к кубатурным формулам, которые являются прямым произведением одномерных квадратурных формул:

(11)

Например, при K=2, если по каждому направлению выбрана обобщённая формула трапеций, а сетка рав­номерная, то веса кубатурной формулы равны соответственно для внутренних, граничных и угловых узлов сетки. Легко показать, что для дважды непрерывно дифференцируемых функций эта формула имеет второй порядок точности, и к ней применим метод Рунге–Ромберга.

Вообще говоря, для разных направлений можно использовать квадратурные формулы разных порядков точности . Тогда главный член погрешности имеет вид:

Желательно для всех направлений использовать квадратурные формулы одинакового порядка точности.

Можно подобрать веса и положение линий сетки так, чтобы одномерная квадратурная формула была точна для многочлена максимальной степени, т.е. была бы

формулой Гаусса, тогда, для случая K=2:

(12)

где –нули многочленов Лежандра и соответствующие веса. Эти формулы рассчитаны на функции высокой гладкости и дают для них большую экономию в числе узлов по сравнению с более простыми формулами.

Произвольная область. Метод последовательного интегрирования можно применять к области про­извольной формы, например, с криволинейной границей. Рассмотрим этот случай при K=2. Для этого проведём через область хорды, па­раллельные оси , и на них введём узлы, расположенные на каждой хорде так, как нам требуется (рис. 4). Представим интеграл в виде:

Сначала вычислим интеграл по вдоль каждой хорды по какой-нибудь одномерной квадратурной фор­муле, используя введённые узлы. Затем вычислим интеграл по ; здесь узлами будут служить проекции хорд на ось ординат.

При вычислении интеграла по имеется одна тонкость. Если область ограничена гладкой кривой, то при длина хорды стремится к нулю не линейно, а как ; значит, вблизи этой точки . То же будет при . Поэтому интегрировать непосредственно по формулам высокого порядка точности бессмысленно. Целесообразно выделить из основную осо­бенность в виде веса , которому соответствуют ортогональные многочлены Чебышева второго рода.

Тогда второе интегрирование выполняется по формулам Гаусса–Кристоффеля:

(13)

где , а и –нули и веса многочленов Чебышева второго рода.

Чтобы можно было применять эту формулу, надо ординаты хорд на рис. 4 заранее выбрать в соответствии с узлами (13). Если это не было сделано, то придётся ограничиться интегрированием по обобщённой формуле трапеций, причём её эффективный порядок точности в этом случае будет ниже второго.

2.4 Кубатурная формула типа Симпсона

Пусть сначала область интегрирования есть K-мерный пространственный параллелепипед (рис. 5), стороны которого параллельны осям координат. Каждый из промежутков разобьём по­полам точками:

, где .

Всего таким образом, получим точек сетки. Имеем:

(14)

Находим K-мерный интеграл, вычисляя каждый внутренний интеграл по квадратурной формуле Симпсона на соответствующем отрезке. Проведём полностью все вычисления для случая K=2:

Применяя к каждому интегралу снова формулу Симпсона, получим:

или

(15)

Формулу (15) будем называть кубатурной формулой Симпсона. Следовательно,

(15)

где – сумма значений подынтегральной функции в вершинах прямоугольника , – сумма значений в серединах сторон прямоугольника , – значение функ­ции в центре прямоугольника . Кратности этих значений обозначены на рис. 5.

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

Опять рассмотрим случай K=2. Положим, что стороны прямоугольника мы разделили соответственно на и равных частей; в результате получилась относительно крупная сеть прямоугольников (на рис. 6 вершины этих прямоугольников отмечены более крупными кружками). Каждый из этих прямоугольников в свою очередь разделим на четыре равные части. Вершины этой последней мелкой сети прямоугольников при­мем за узлы кубатурной формулы.

Пусть и . Тогда сеть узлов будет иметь следующие координаты:

и

Для сокращения введём обозначение

Применяя формулу (15) к каждому из прямоугольников крупной сети, будем иметь (рис.6):

Отсюда, делая приведение подобных членов, окончательно находим:

(16)

где коэффициенты являются соответствующими элементами матрицы

Если область интегрирования – произвольная, то строим параллелепипед , стороны которого параллельны осям координат (рис. 83). Рассмотрим вспомогательную функцию

В таком случае, очевидно, имеем:

Последний интеграл приближённо может быть вычислен по общей кубатурной формуле (16).

2.5 Принципы построения программ с автоматическим выбором шага

При написании программ численного интегрирования желательно, чтобы для любой функции распределение узлов являлось оптимальным или близким к нему. Однако в случае резко меняющихся функций возникают некоторые проблемы. Если первоначальная сетка, на которой исследуется подынте­гральная функция, частая, то сильно загружается память ЭВМ; если она редкая, то не удаётся хорошо аппроксимировать оптимальное распределение узлов на участках резкого изменения подынтегральной функции. Рассмотрим некоторые из процедур распределения узлов интегрирования, обеспечивающие лучшее приближение к оптимальному распределению узлов для функций с особенностями.

Пусть на элементарном отрезке интегрирования вычисляется приближённое значе­ние интеграла и мера погрешности . Требуется вычислить . Первая процедура, которую естественно назвать горизонтальной, определяется заданием параметров . Полагаем . Предположим, что каким-то образом уже вычислено приближён­ное значение интеграла . Программа располагает в каждый момент времени некоторым значе­нием , с которым надо начинать считать оставшуюся часть интеграла. Вычисляем величину , соответствующую отрезку . Если оказалось , то вычисляем приближённое значение и полагаем . Мы получили приближённое значение величины . В случае полагаем , в противном случае полагаем . Мы готовы к следующему шагу. Если оказалось , то принимаем за новое значение вели­чины и возвращаемся к исходной позиции: вычислено значение интеграла и задан шаг . Начальные условия для применения процедуры:

Процедура должна также иметь блок окончания работы: если оказалось, что , то следует положить . Установилась практика брать .

Другая процедура, которую можно назвать вертикальной, определяется заданием числа и заключается в следующем. Пусть на каком-то шаге возникает необходимость вычисления интеграла по отрезку разбиения : ; вычисляется величина , соответствующая этому отрезку. Если она оказалась меньше , то этот интеграл вычисляется по соответствующей формуле и программа переходит к следующему справа отрезку разбиения. В противном случае отрезки и объявляются отрезками разбиения, и программа обращается к вычислению интеграла по ле­вому из этих отрезков. В начале работы программа обращается к вычислению исходного интеграла . Некоторым недостатком этой процедуры является необходимость запоминания отрезков раз­биения, интегрирование по которым на данный момент не произведено.

3 Список использованной литературы.

1. Бахвалов Н.С. Численные методы. т.1 – М.: Наука. 1975.

2. Демидович Б.П., Марон И.А. Основы вычислительной математики. – М.: Наука, 1966.

3. Калиткин Н.Н Численные методы. – М.: Наука, 1978.

4. Мусіяка В.Г. Основи чисельних методів механіки. – Дніпропетровськ: Видавництво ДДУ, 1993.

4 Практическая часть

4.1 Решение задачи

Н

1

0.5

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

Н

0.5

1

0.75

е учитывая граничные ячейки, получаем: .

Дополнение от граничных ячеек: .

Окончательно получаем:

4.2 Блок-схема программы

За программой и блок-схемой по данной теме обращайтесь по адресу: shuric_1@mail.ru

4.4 Результаты решения

Расчёт проводился при точности eps=1E-6.

Интеграл равен: 0.221612

Количество ячеек равно 8525.