Вычисления определенного интеграла с помощью ф. – лы Симпсона на компьютере
МОСКОВСКИЙ
ГОСУДАРСТВЕННЫЙ СТРОИТЕЛЬНЫЙ УНИВЕРСИТЕТ
КУРСОВАЯ РАБОТА
«Программа приближенного вычисления определенного
интеграла с помощью ф – лы Симпсона на компьютере»
Выполнил:
студент ф – та ЭОУС – 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). После указанных выше действий получится “большая” формула Симпсона, которая имеет вид,
I
»
h / 3 * (Yкр + 2 * Yнеч + 4 * Yчет) (2)
|
|
где Yкр = y1 + yn, Yнеч
= y3 + y5 + … +
yn – 1, Yчет = y2 + y4 + … + yn – 2, а h = (b – a) / n.
Задача. Пусть нужно проинтегрировать
функцию f(x) = x³(x - 5)² на отрезке [0, 6] (рис.
2). На этом отрезке функция непрерывна и принимает только неотрицательные
значения, т. е. знакопостоянна.
рис. 2
Для
выполнения поставленной задачи составлена нижеописанная программа, приближенно
вычисляющая определенный интеграл с помощью формулы Симпсона. Программа состоит
из трех функций main, f и integral. Функция main вызывает функцию integral для вычисления интеграла и
распечатывает на экране результат. Функция f принимает аргумент x типа float
и возвращает значение интегрируемой функции в этой точке. Integral – основная функция программы: она
выполняет все вычисления, связанные с нахождением определенного интеграла. Integral принимает четыре параметра: пределы
интегрирования типа float,
допустимую относительную ошибку типа float и указатель на интегрируемую функцию. Вычисления выполняются до тех
пор, пока относительная ошибка, вычисляемая по формуле
| (In/2 – In) / In | ,
где In интеграл при числе разбиений n, не будет меньше требуемой.
Например, допустимая относительная ошибка e = 0.02 это значит, что максимальная погрешность в
вычислениях будет не больше, чем In *
e = 0.02 * In. Функция реализована с
экономией вычислений, т. е. учитывается, что Yкр постоянная, а Yнеч = Yнеч + Yчет, поэтому эти значения вычисляются единожды. Высокая
точность и скорость вычисления делают использование программы на основе формулы
Симпсона более желательным при приближенном вычислении интегралов, чем
использование программ на основе формулы трапеции или метода прямоугольников.
Ниже предлагается блок – схема, спецификации, листинг и
ручной счет программы
на примере поставленной выше задачи. Блок
– схема позволяет
отследить и понять особенности алгоритма программы, спецификации дают представление о назначении каждой переменной в
основной функции integral, листинг - исходный код работающей программы с комментариями, а ручной счет предоставляет возможность проанализировать результаты
выполнения программы.
2.
Блок – схема
программы
ДА
НЕТ
|
|
|
3.
Спецификации
Имя
переменной
|
Тип
|
Назначение
|
n
|
int
|
Число
разбиений отрезка [a,
b]
|
i
|
int
|
Счетчик циклов
|
a
|
float
|
Нижний
предел интегрирования
|
b
|
float
|
Верхний
предел интегрирования
|
h
|
float
|
Шаг
разбиения отрезка
|
e
|
float
|
Допустимая
относительная ошибка
|
f
|
float (*)
|
Указатель
на интегрируемую фун - цию
|
s_ab
|
float
|
Сумма
значений фун – ции в точках a
и b
|
s_even
|
Сумма
значений фун – ции в нечетных точках
|
s_odd
|
float
|
Сумма
значений фун – ции в четных точках
|
s_res
|
float
|
Текущий
результат интегрирования
|
s_pres
|
float
|
Предыдущий
результат интегрирования
|
4.
Листинг
программы
#include
<stdio.h>
#include
<math.h>
/* Прототип фун – ции, вычисляющей
интеграл */
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 < n; i
+= 2) {
s_even += f(a + i * h);
}
do
{
s_odd
= 0;
s_pres
= s_res;
/* Сумма значений фун – ции в четных
точках */
for (i = 1; i < n; i += 2) {
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
|
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
|