МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ СТРОИТЕЛЬНЫЙ УНИВЕРСИТЕТ
КУРСОВАЯ РАБОТА
«Программа приближенного вычисления определенного интеграла с помощью ф – лы Симпсона на компьютере»
Выполнил:
студент ф – та ЭОУС – 1 – 12
Валюгин А. С.
Принял:
Зоткин С. П.
Москва 2001
1. Введение
Определенный интеграл от функции, имеющей неэлементарную первообразную, можно вычислить с помощью той или иной приближенной формулы. Для решения этой задачи на компьютере, среди
прочих, можно воспользоваться формулами прямоугольников, трапеций или формулой Симпсона. В данной работе рассматривается
именно последняя.
Рассмотрим функцию y = f(x). Будем считать, что на отрезке [a, b] она положительна и непрерывна.
Найдем площадь криволинейной трапеции aABb (рис. 1).
рис. 1
Для этого разделим отрезок [a, b] точкой c = (a + b) / 2 пополам и в точке C(c, f(c)) проведем
касательную к линии y = f(x). После этого разделим [a, b] точками p и q на 3 равные части и проведем через них прямые x = p и x = q. Пусть P и Q – точки пересечения этих прямых с касательной.
Соединив A с P и B с Q, получим 3 прямолинейные трапеции aAPp, pPQq, qQBb. Тогда площадь трапеции aABb можно
приближенно посчитать по следующей формуле
I » (aA + pP) / 2 * h + (pP + qQ) / 2 * h + (qQ +
bB) / 2 * h, где h = (b – a) / 3.
Откуда получаем
I » (b – a) / 6 * (aA + 2 * (pP + qQ) + bB)
заметим, что aA = f(a), bB = f(b), а pP + qQ = 2 * f(c), в итоге получаем малую
фор – лу Симпсона
I » (b – a) / 6 * (f(a) + 4 * f(c) + f(b)) (1)
Малая формула Симпсона дает интеграл с хорошей точностью, когда график
подинтегральной функции мало изогнут, в случаях же, когда дана более сложная функция малая формула Симпсона непригодна. Тогда, чтобы посчитать интеграл
заданной функции нужно разбить отрезок [a, b] на n частей и к каждому из отрезков применить формулу
(1). После указанных выше действий получится “большая” формула Симпсона, которая имеет вид,
ДА
НЕТ
3. Спецификации
Имя переменной
Тип
Назначение
n
int
Число разбиений отрезка [a, b]
i
int
Счетчик циклов
a
float
Нижний предел интегрирования
b
float
Верхний предел интегрирования
h
float
Шаг разбиения отрезка
e
float
Допустимая относительная ошибка
f
float (*)
Указатель на интегрируемую фун - цию
s_ab
float
Сумма значений фун – ции в точках a и b
s_even
float
Сумма значений фун – ции в нечетных точках
s_odd
float
Сумма значений фун – ции в четных точках
s_res
float
Текущий результат интегрирования
s_pres
float
Предыдущий результат интегрирования
4. Листинг программы
#include
#include
/* Прототип фун – ции, вычисляющей интеграл */
float integral(float, float, float, float (*)(float));
/* Прототип фун – ции, задающей интегрируемую фун – цию */
float f(float);
main()
{
float result;
result = integral(0, 6, .1, f);
printf("%f", result);
return 0;
}
/* Реализация фун – ции, задающей интегрируемую фун – цию */
float f(float x)
{
/* Функция f(x) = x³(x-5)² */
return pow(x, 3) * pow(x - 5, 2);
}
/* Реализация фун – ции, вычисляющей интеграл */
float integral(float a, float b, float e, float (*f)(float))
{
int n = 4, i; /* Начальное число разбиений 4 */
float s_ab = f(a) + f(b); /* Сумма значений фун – ции в a и b */
float h = (b – a) / n; /* Вычисляем шаг */
float s_even = 0, s_odd;
float s_res = 0, s_pres;
/* Сумма значений фун – ции в нечетных точках */
for (i = 2; i
s_even += f(a + i * h);
}
do {
s_odd = 0;
s_pres = s_res;
/* Сумма значений фун – ции в четных точках */
for (i = 1; i
s_odd += f(a + i * h);
}
/* Подсчет результата */
s_res = h / 3 * (s_ab + 2 * s_even + 4 * s_odd);
/* Избегаем деления на ноль */
if (s_res == 0) s_res = e;
s_even += s_odd;
n *= 2;
h /= 2;
} while (fabs((s_pres - s_res) / s_res) > e);/* Выполнять до тех пор, пока результат не будет удовлетворять
допустимой ошибке */
return fabs(s_res); /* Возвращаем результат */
}
5. Ручной счет
Таблица константных значений для n = 8
Имя переменной
Значение
a
0
b
6
e
.1
s_ab
216
h
.75
Подсчет s_even
i
a + i * h
f(a + i * h)
s_even
2
1.5
41.34375
41.34375
4
3
108
149.34375
6
4.5
22.78125
172.125
Подсчет s_odd
i
a + i * h
f(a + i * h)
s_odd
1
.75
7.62012
7.62012
3
2.25
86.14158
93.7617
5
3.75
82.3973
176.159
7
5.25
9.044
185.203
Подсчет s_res
ò f(x) dx
s_res = h / 3 * (s_ab + 2 * s_even + 4 * s_odd)
Абсолютная ошибка
324
325.266
1.266