|
Курсовая: Вычисление кратных интегралов методом ячеек с автоматическим выбором шага |
| |
|
|
|
Курсовая: Вычисление кратных интегралов методом ячеек с автоматическим выбором шага
Курсовая: Вычисление кратных интегралов методом ячеек с автоматическим выбором шага
Министерство образования Украины
Днепропетровский государственный университет
–––––––––––––––––––––––––––––––––––––––––––––
Факультет прикладной математики
Кафедра вычислительной механики и прочности конструкций
КУРСОВАЯ РАБОТА
по численным методам в механике
на тему
Вычисление кратных интегралов
методом ячеек
с автоматическим выбором шага
Исполнитель: студент группы ПД-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
|
 |
|
