Лабораторна робота № 9
Програмування чисельного інтегрування функції
1. Порядок виконання роботи
1.1. Скласти алгоритмічною мовою Фортран програму чисельного інтегрування функції.
1.2. Відлагодити на комп’ютері складену програму.
1.3. Ввести числові дані та отримати результат.
1.4. Скласти звіт про роботу й захистити його.
2. Вказівки до виконання роботи
2.1. Чисельне інтегрування функції
В задачах електромеханіки досить часто виникає необхідність обчислення означених інтегралів.
Якщо функція неперервна на якомусь інтервалі і відома її первісна, то означений інтеграл від цієї функції в границях інтервалу можна обчислити за формулою Ньютона-Лейбніца. Однак на практиці у багатьох випадках інтегрування функції внаслідок її складності пов’язане із значними труднощами, а в окремих випадках нездійсненне. Крім того, дуже часто підінтегральна функція задається таблично (іноді графічно) і тоді можливість аналітичного інтегрування відпадає. У таких випадках застосовуються методи чисельного інтегрування, тобто знаходження означеного інтеграла чисельними методами.
Чисельне інтегрування є наближеним і полягає в обчисленні означеного інтеграла на підставі використання протабульованої функції, тобто значень підінтегральної функції у вузлах.
Для чисельного інтегрування функції, заданої аналітично, її необхідно перед інтегруванням протабулювати.
Геометрична інтерпретація означеного інтеграла – це площа, яка обмежується графіком функції , віссю абсцис і вертикальними прямими що відповідають абсцисам а, b границь інтервалу інтегрування (рис. 19).
Серед чисельних методів інтегрування широке застосування отримали методи прямокутників, трапецій і Сімпсона (парабол).
2.1.1. Метод прямокутників
Метод прямокутників передбачає обчислення означеного інтеграла за формулою (варіант 1)

де N – кількість рівновіддалених вузлів на інтервалі інтегрування; – значення підінтегральної функції у її вузлах, – крок аргументу функції. За цим методом підінтегральна функція замінюється ступінчастою лінією, а сам означений інтеграл – сумою площ прямокутників, як зображено на рис. 19.
Заміна підінтегральної функції ступінчастою лінією у методі прямокутників приводить до обчислення інтеграла з деякою похибкою. Так, як видно з рис. 19, для заданої функції отримане значення інтеграла буде меншим, бо при його обчисленні не враховуються площі криволінійних трикутників на ділянках зростання функції. Очевидно, що чим менша величина кроку, тим вища точність обчислення інтеграла.
Рис. 19. Геометрична інтерпретація означеного інтеграла методом прямокутників (варіант 1)
Означений інтеграл також може бути обчислений методом прямокутників за формулою (варіант 2)

Геометрична інтерпретація означеного інтеграла для цього випадку зображена на рис. 20.
Рис. 20. Геометрична інтерпретація означеного інтеграла методом прямокутників (варіант 2)
Приклад. Для заданої функції

скласти програму її чисельного інтегрування методом прямокутників у діапазоні від a до b при заданій у цьому діапазоні кількості N рівновіддалених вузлів.
Один із можливих варіантів програми:
С
С
С
С
С
Лабораторна робота № 9
Обчислення означеного інтеграла методом прямокутників
Виконав ст. гр. ЕМ–11 Вдалий Б.

C
Описання змінної дійсного типу



REAL INT

C
Оголошення розміру масиву



DIMENSION Y(100)

C
Введення вхідних даних у діалоговому режимі



1
2
3
WRITE (*, 1)
FORMAT (2X, ‘A=’,$)
READ (*, *) A
WRITE (*, 2)
FORMAT (2X, ‘B=’,$)
READ (*, *) B
WRITE (*, 3)
FORMAT (2X, ‘N=’,$)
READ (*, *) N

C
Обчислення величини кроку



H = (B–A)/(N–1)

C
Задавання біжучого значення аргументу вузла



X=A

С
“Очищення“ комірки пам’яті, де буде вестися підрахунок суми



S = 0.

C
Обчислення суми значень функції у її вузлах від 1 до N–1



DO 4 I = 1, N–1

C
Обчислення значення функції у n-му вузлі



Y(I) = COS(ALOG(X)**2/3.)

C
Обчислення суми значень функції у її вузлах від 1 до n



S = S + Y(I)

C
Обчислення значення аргументу в (n+1)-му вузлі


4
X = X + H

C
Обчислення означеного інтеграла



INT = S*H

C
С
Друкування обчисленого означеного інтеграла
на екран монітора



WRITE (*, *) ‘ INT =’, INT

C
Закінчення роботи програми



STOP
END


Обчислення означеного інтеграла може здійснювати і підпрограма. Один із можливих варіантів такої програми з використанням підпрограми обчислення методом прямокутників означеного інтеграла, який передбачає використання для обчислення підінтегральної функції підпрограми-функції, має вигляд:
С
С
C
C
C
C
Обчислення означеного інтеграла методом прямокутників
з використанням підпрограми для обчислення інтеграла
та підпрограми-функції для обчислення підінтегральної функції
Основна програма

C
Введення вхідних даних у діалоговому режимі



1
2
3
WRITE (*, 1)
FORMAT (2X, ‘A=’,$)
READ (*, *) A
WRITE (*, 2)
FORMAT (2X, ‘B=’,$)
READ (*, *) B
WRITE (*, 3)
FORMAT (2X, ‘N=’,$)
READ (*, *) N

C
C
Звертання до підпрограми PR
обчислення методом прямокутників означеного інтеграла



CALL PR(A, B, N, Spr)

C
С
Друкування обчисленого означеного інтеграла
на екран монітора



WRITE (*, *) ‘Spr =’, Spr

C
Закінчення роботи програми



STOP
END



C
C
Підпрограма PR
обчислення методом прямокутників означеного інтеграла



SUBROUTINE PR(A, B, N, Spr)

C
**********************************************************

C

Опис формальних параметрів.

C
Вхідні величини:

C
C
C

A – значення лівої границі інтервалу інтегрування;
B – значення правої границі інтервалу інтегрування;
N – кількість рівновіддалених вузлів на інтервалі інтегрування.

C
Вихідні величини:

C

Spr – величина означеного інтегралу.

C
**********************************************************

C
Обчислення величини кроку



H = (B–A)/(N–1)

С
“Очищення“ комірки пам’яті, де буде вестися підрахунок суми



S = 0.

C
Обчислення суми значень функції у її вузлах від 2 до N



4
DO 4 I = 2, N
X = A + (I–1)*H
S = S + FUN(X)

C
Обчислення означеного інтеграла



Spr = S*H

C
Закінчення роботи підпрограми



RETURN
END




C
C
Підпрограма-функція FUN(X)
обчислення підінтегральної функції



FUNCTION FUN(X)
FUN = COS(ALOG(X)**2/3.)
RETURN
END

2.1.2. Метод трапецій
При обчисленні означеного інтеграла методом трапецій він наближено замінюється площею трапецій, утворених ламаною лінією, яка апроксимує задану підінтегральну функцію (рис. 21). У цьому випадку означений інтеграл обчислюється за формулою

З порівняння рис. 19, 20 і 21 видно, що метод трапецій дає точніші результати обчислення означеного інтеграла, ніж метод прямокутників.
Рис. 21. Геометрична інтерпретація означеного інтеграла методом трапецій
Один із можливих варіантів фрагменту програми:
C
C
C
Зміна біжучого значення аргументу вузла,
обчислення біжучого значення функції,
обчислення суми значень функції в усіх її вузлах



4
DO 4 I = 1, N
X = A + (I–1)*H
Y(I) = COS(ALOG(X)**2/3.)
S = S + Y(I)

C
Обчислення означеного інтеграла



INT = (S – (Y(1)+Y(N))/2.) * H


Тут слід звернути увагу на те, що змінна S містить суму значень функції у її всіх вузлах, а формула обчислення означеного інтеграла методом трапецій передбачає, що значення функції для першого та останнього вузла повинно множитися на коефіцієнт 1/2.
2.1.3. Метод Сімпсона
Ще точніше, у порівнянні з методом трапецій, можна обчислити означений інтеграл методом Сімпсона, за яким трапеції розглядаються як криволінійні. У межах кожної такої трапеції підінтегральна функція замінюється інтерполяційним поліномом другого порядку, побудованим за значеннями функції в трьох вузлах трапеції: на її початку, в середині та кінці. Формула Сімпсона обчислення означеного інтеграла має вигляд:

Увага! У методі Сімпсона кількість N вузлів повинна бути непарною. Тут слід звернути увагу на коефіцієнт (1, 4 чи 2), на який множиться значення функції у n-му вузлі: для усіх парних вузлів це коефіцієнт 4, для першого та останнього вузла коефіцієнт 1, а для решти непарних вузлів – коефіцієнт 2.
Один із можливих варіантів фрагменту програми, де передбачено вибір коефіцієнта, на який множиться значення функції у вузлі:
C
C
C
C
Зміна біжучого значення аргументу вузла,
обчислення біжучого значення функції,
вибір коефіцієнта,
обчислення суми значень функції в усіх її вузлах з врахуванням вибраного коефіцієнта



4
DO 4 I = 1, N
X = A + (I–1)*H
Y(I) = COS(ALOG(X)**2/3.)
AK = 2.
IF (I/2*I .EQ. I) AK = 4.
S = S + Y(I)*AK

C
Обчислення означеного інтеграла



INT = (S – (Y(1)+Y(N)) * H/3.


У цьому фрагменті програми оператор
IF (I/2*I .EQ. I) AK = 4.
присвоює значення коефіцієнту AK = 4. лише для парних І. Для непарних I (див. попередній оператор) він залишається рівним 2.
2.1.4. Обчислення методом Сімпсона означеного інтеграла із заданою точністю
Вибір більшої кількості N вузлів на інтервалі інтегрування дає можливість підвищити точність обчислення значення інтеграла, однак вимагає збільшення об’єму пам’яті комп’ютера і часу на обчислення. Тому доцільно при розрахунку означеного інтеграла задаватися точністю обчислень ?.
Для обчислення означеного інтеграла із заданою точністю ? організовують ітераційний цикл, збільшуючи на кожній ітерації кількість М = N–1 ділянок розбиття заданого інтервалу у два рази. Так, якщо на першій ітерації кількість ділянок розбиття є М1, то на другій ітерації ця кількість збільшується вдвічі, тобто М2 = 2М1 і т.д. На і-тій ітерації

Ітераційний процес закінчується, коли модуль різниці наближених значень означеного інтеграла і-ої та -ої ітерації не перевищує заданої точності ?.
Далі наведено текст підпрограми SIMP обчислення методом Сімпсона означеного інтеграла із заданою точністю. Ця підпрограма міститься у програмному забезпеченні даної дисципліни.
C
C
C
Підпрограма SIMP
обчислення методом Сімпсона означеного інтеграла
із заданою точністю



SUBROUTINE SIMP(A, B, EPS, IP, Simt, N)

C
**********************************************************

C

Опис формальних параметрів.

C
Вхідні величини:

C
C
C
C

A – значення лівої границі інтервалу інтегрування;
B – значення правої границі інтервалу інтегрування;
EPS – задана точність;
IP – ознака друку:

C
C


IP = 0 – не виводити хід ітераційного процесу;
IP = 1 – виводити хід ітераційного процесу.

C
Вихідні величини:

C
C

Simt – величина означеного інтегралу;
N – кількість рівновіддалених вузлів на інтервалі інтегрування,

C


яка забезпечує задану точність.

C
**********************************************************

C
Кількість ділянок розбиття інтервалу інтегрування



M = 2

C
Величина кроку аргументу функції



H = (B–A)/M

C
Обчислення суми значень підінтегральної функції на границях інтервалу



AB = FUN(A) + FUN(B)

C
Обчислення означеного інтегралу при М = 2



S0 = (AB + 4*FUN(A+H))*H/3.

C
C
Виведення кількості вузлів та обчисленого значення інтегралу
на екран монітора при умові, що IP = 1



1
IF (IP .EQ. 1) WRITE (*, 1) M+1, S0
FORMAT(2X, ‘N = ’, I3, 3X, ‘INTEGRAL = ’, E12.5)

C
Обчислення нового значення кількості ділянок розбиття


2
M = 2*M

C
Обчислення нового значення величини кроку



H = 0.5*H

C
Обчислення нового значення означеного інтегралу



3
X = A + H
S = AB +4.*FUN(X)
MM3 = M – 3
I = 1
I = I + 1
X = X + H
S = S + 2*FUN(X)
I = I + 1
X = X + H
S = S +4.*FUN(X)
IF (I .LE. MM3) GOTO 3
Simt = S*H/3.
N = M+1

C
C
Виведення кількості вузлів та обчисленого значення інтегралу
на екран монітора при умові, що IP = 1



IF (IP .EQ. 1) WRITE (*, 1) N, Simt

C
Перевіряння умови досягнення заданої точності



IF (ABS(Simt – S0) .LE. EPS) RETURN

C
Присвоєння змінній S0 нового значення означеного інтегралу



S0 = Simt

C
Перехід на виконання нової ітерації



GOTO 2

C
Кінець програмної одиниці



STOP
END


2.2. Завдання
Скласти програму обчислення означеного інтеграла чисельними методами прямокутників, трапецій та Сімпсона. Ця програма повинна складатися із основної програми та підпрограм. Основна програма має бути призначена для введення вхідної інформації і виведення результатів обчислень, а кожна із підпрограм – для обчислення означеного інтеграла одним із заданих методів. У цій програмі слід використати наведену підпрограму SIMP обчислення означеного інтеграла методом Сімпсона із заданою точністю. Для обчислення підінтегральної функції – використати підпрограму-функцію.
Порівняти отримані різними методами значення означеного інтегралу.
Дані для розрахунку наведені в табл. 10. Точність ? = 0,0001.
Таблиця 10
Дані для обчислення означеного інтеграла
№ варіанту

Підінтегральна функція

a

b

N

1
2
3
4
5

1

2,4
5,1
28

2

0,3
2,2
100

3

1,2
4,1
30

4

4,6
7,6
31

5

0,0
2,0
21

6

0,4
5,3
50

7

–0,5
2,5
31

8

–1,5
1,4
30

9

2,6
5,5
30

10

1,5
4,5
31

11

–0,4
1,6
41

12

1,4
7,2
59

1
2
3
4
5

13

–0,6
2,0
27

14

0,0
2,5
51

15

–1,5
1,5
31

16

0,0
2,0
21

17

0,0
6,5
27

18

0,5
5,0
91

19

0,0
1,3
27

20

5,5
15,5
41

21

1,6
4,8
63

22

3,5
6,5
31

23

1,0
3,0
41

24

0,0
1,5
76

25

1,5
5,5
81

26

2,6
6,3
38

27

1,5
4,5
31

28

4,5
8,5
41

29

0,4
1,2
41

30

3,5
6,5
31


Список літератури
1. Алексеев В. Е., Ваулин А. С., Петрова Г. Б. Вычислительная техника в инженерных и экономических расчетах. Сборник задач и упражнений: Учеб. пособие для вузов / Под ред. А. В. Петрова. – М.: Высш. шк., 1984. – 136 с., ил.
2. Бартеньев О. В. Современный Фортран. – 2-е изд., испр. – М.: «Диалог-МИФИ», 1998. – 397 с.
3. Бартеньев О. В. Фортран для студентов. – М.: «Диалог-МИФИ», 1999. – 400 с.
4. Бартеньев О. В. Visual Fortran: новые возможности. – М.: «Диалог-МИФИ», 1999. – 304 с.
5. Белецки Я. Фортран 77: Пер. с польск. О. И. Гуськовой / Под ред. В. Р. Носова; Предисл. В. Р. Носова. – М.: Высш. шк., 1991. – 207 с.: ил.
6. Бронштейн И. Н., Семендяев К. А. Справочник по математике. – 13-е изд., исправленное. – М.: Наука, 1986. – 544 с.
7. Демидович Б. П., Марон И. А. Основы вычислительной математики. – М.: Наука, 1970. – 664 с.
8. Корн Г., Корн Т. Справочник по математике. – 4-е издание. – М.: Наука, 1978. – 832 с.
9. Программирование на Фортране 77: Пер. с англ./ Дж. Ашкрофт, Р. Элдридж, Р. Полсон, Г. Уилсон. – М.: Радио и связь, 1990. – 272 с.: ил.
10. Рудавський Ю. К., Костробій П. П., Луник Х. П., Уханська Д. В. Лінійна алгебра та аналітична геометрія: Навч. посібник. – Львів: Видавництво Державного Університету “Львівська політехніка”, 1999. – 262 с.
11. Соловьев П. В. Fortran для персонального компьютера. – М.: Арист, 1991. – 223 с.
12. Уорд Т., Бромхед Э. Фортран и искусство программирования персональных ЭВМ: Пер. с англ. – М.: Радио и связь, 1993. – 352 с.: ил.
13. Фигурнов В. Э. IBM для пользователя. – М.: Финансы и статистика, 1990. 240 с.: ил.
14. Штыков В. В. FORTRAN & WIN32 API: Создание программного интерфейса для Windows средствами современного Фортрана. – М.: «Диалог-МИФИ», 2001. – 304 с.