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

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

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

тАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУтАУ

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

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

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

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

на тему

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

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

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

Исполнитель: студент группы ПД-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. Мусiяка В.Г. Основи чисельних методiв механiки. тАУ Днiпропетровськ: Видавництво ДДУ, 1993.

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

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

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

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

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

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

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

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

4.4 Результаты решения Расчёт проводился при точности eps=1E-6. Интеграл равен: 0.221612 Количество ячеек равно 8525.

Вместе с этим смотрят:

Вычисление кратных интегралов методом ячеек с автоматическим выбором шага
Вычисление определенного интеграла
Вычисление определенного интеграла методами трапеций и средних прямоугольников
Вычисление определителя методом Гаусса с выбором главного элемента