Вычисления определенного интеграла с помощью ф. – лы Симпсона на компьютере

МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ СТРОИТЕЛЬНЫЙ УНИВЕРСИТЕТ

































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

«Программа приближенного вычисления определенного интеграла с помощью ф – лы Симпсона на компьютере»





Выполнил:

 студент ф – та ЭОУС – 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), в итоге получаем малую фор – лу Симпсона 



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, листинг -  исходный код работающей программы с комментариями, а ручной счет предоставляет возможность проанализировать результаты выполнения программы. 

Блок – схема программы



Ввод a, b, e, f(x)









n = 4, h = (b – a) / n

s_ab = f(a) + f(b)

s_even = 0, s_res = 0













i = 2, n – 1, 2





s_even = s_even + f(a + i * h)









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                                                                                                          ДА







                                                                          НЕТ



Вывод s_res



		

	

	

	

	

	

Спецификации



Имя переменной

Тип

Назначение

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 <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); /* Возвращаем результат */

}

		

	

Ручной счет



Таблица константных значений для 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