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

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

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

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

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

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

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

на тему

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

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

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

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

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

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

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

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

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

PAGE 9

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
NURBIZ.KZ - каталог компаний и предприятий Казахстана и Алматы

Фараби Холл / Farabi Hall

Скидка 100%

Специальная подарочная цена на проведение банкетов от 50 до 300 человек!

Что мешает запоминать информацию – учимся на ошибках

Декретный отпуск не помеха для получения образования