МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ СТРОИТЕЛЬНЫЙ УНИВЕРСИТЕТ
КУРСОВАЯ РАБОТА
«Программа приближенного вычисления определенного интеграла с помощью ф
– лы Симпсона на компьютере»
Выполнил:
студент ф – та ЭОУС – 1 – 12
Валюгин А. С.
Принял:
Зоткин С. П.
Москва 2001
Введение
Определенный интеграл от функции, имеющей неэлементарную первообразную,
можно вычислить с помощью той или иной приближенной формулы. Для решения
этой задачи на компьютере, среди прочих, можно воспользоваться формулами
прямоугольников, трапеций или формулой Симпсона. В данной работе
рассматривается именно последняя.
Рассмотрим функцию 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), в итоге
получаем малую фор – лу Симпсона
Малая формула Симпсона дает интеграл с хорошей точностью, когда график
подинтегральной функции мало изогнут, в случаях же, когда дана более
сложная функция малая формула Симпсона непригодна. Тогда, чтобы
посчитать интеграл заданной функции нужно разбить отрезок [a, b] на n
частей и к каждому из отрезков применить формулу (1). После указанных
выше действий получится “большая” формула Симпсона, которая имеет вид,
где 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, листинг – исходный код работающей программы с комментариями,
а ручной счет предоставляет возможность проанализировать результаты
выполнения программы.
Блок – схема программы
ДА
НЕТ
Спецификации
Имя переменной Тип Назначение
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 Предыдущий результат интегрирования
Листинг программы
#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 e);/* Выполнять до тех пор, return fabs(s_res); /* Возвращаем результат */ } Ручной счет Таблица константных значений для 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 Ввод a, b, e, f(x) n = 4, h = (b – a) / n s_ab = f(a) + f(b) s_even = 0, s_res = 0 s_even = s_even + f(a + i * h) i = 2, n – 1, 2 s_odd = 0, s_pres = s_res i = 1, n – 1, 2 s_odd = s_odd + f(a + i * h) s_res = h / 3 * (s_ab + 2 * s_even + 4 * s_odd) s_even = s_even + s_odd, n = n / 2, h = h / 2 | (s_pres – s_res) / s_res | > e I ( (b – a) / 6 * (f(a) + 4 * f(c) + f(b)) (1) Вывод s_res I ( h / 3 * (Yкр + 2 * Yнеч + 4 * Yчет) (2)
пока результат не будет удовлетворять допустимой ошибке */
ошибка
Нашли опечатку? Выделите и нажмите CTRL+Enter
