Рефераты. Курсовая: Вычисление кратных интегралов методом ячеек с автоматическим выбором шага

Курсовая: Вычисление кратных интегралов методом ячеек с автоматическим выбором шага

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

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

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

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

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

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

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

на тему

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

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

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

Исполнитель: студент группы ПД-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)

, будем иметь:

(4)

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

получаем:

2.2 Метод ячеек

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

(5)

соответственно площадь ячейки и координаты её центра, получим:

(6)

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

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

(7)

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

(8)

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

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

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

(9)

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

–координаты центра тяжести, вычисляемые по обычным формулам:

(10)

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

; этот объём вычислим приближённо. Эти площади подставим в (6) и
вычислим интеграл.

, если функция дважды непрерывно дифференцируема; это означает второй
порядок точности.

, и для хорошей точности потребуется более подробная сетка.

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

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

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

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

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

(11)

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

. Тогда главный член погрешности имеет вид:

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

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

(12)

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

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

; здесь узлами будут служить проекции хорд на ось ординат.

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

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

(13)

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

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

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

разобьём пополам точками:

.

точек сетки. Имеем:

(14)

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

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

или

(15)

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

(15()

. Кратности этих значений обозначены на рис. 5.

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

кубатурной формулы.

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

и

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

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

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

(16)

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

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

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

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

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

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

. Начальные условия для применения процедуры:

.

. Некоторым недостатком этой процедуры является необходимость
запоминания отрезков разбиения, интегрирование по которым на данный
момент не произведено.

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

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

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

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

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

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

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

.

.

.

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

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

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

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

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

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

O

x

y

Рис. 1

M1

MN

(()

O

b1

a1

a2

b2

x1

x2

Рис. 2

Рис. 3

(2(x2)

(1(x2)

O

a2

b2

x1

x2

Рис. 4

=A1

Рис. 5

O

x1

x2

(R)

1

1

1

1

4

4

4

4

16

h1

h1

k1

k1

y1

y1

y1

=A1

Рис. 6

O

x1

x2

R

O

x

y

Рис. 7

R

(

1

2

3

x1

x2

O

G

0.5

1

0.75

1

0.5