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

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

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

Выполнил:

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

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

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

                                                                                                         
ДА

                                                                         
НЕТ

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

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

Подсчет s_even

Подсчет
s_odd

Подсчет s_res