Метод Симпсона на компьютере

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

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

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

лы Симпсона на компьютере»

Выполнил:

студент ф – та ЭОУС – 1 – 12

Валюгин А. С.

Принял:

Зоткин С. П.

Москва 2001

1. Введение

Определенный интеграл от функции, имеющей неэлементарную первообразную,

можно вычислить с помощью той или иной приближенной формулы. Для решения

этой задачи на компьютере, среди прочих, можно воспользоваться формулами

прямоугольников, трапеций или формулой Симпсона. В данной работе

рассматривается именно последняя.

Рассмотрим функцию y = f(x). Будем считать, что на отрезке [a, b] она

положительна и непрерывна. Найдем площадь криволинейной трапеции aABb (рис.

1).

[pic]

рис. 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). На этом отрезке функция непрерывна и принимает

только неотрицательные значения, т. е. знакопостоянна.

[pic]

рис. 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 |float |Сумма значений фун – ции в нечетных |

| | |точках |

|s_odd |float |Сумма значений фун – ции в четных |

| | |точках |

|s_res |float |Текущий результат интегрирования |

|s_pres |float |Предыдущий результат интегрирования |

4. Листинг программы

#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 < 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 |.1 |

|s_ab |216 |

|h |.75 |

Подсчет s_even

|i |a + i * |f(a + i *|s_even |

| |h |h) | |

|2 |1.5 |41.34375 |41.34375 |

|4 |3 |108 |149.34375|

|6 |4.5 |22.78125 |172.125 |

Подсчет s_odd

|i|a + i *|f(a + i |s_odd |

| |h |* h) | |

|1|.75 |7.62012 |7.6201|

| | | |2 |

|3|2.25 |86.14158|93.761|

| | | |7 |

|5|3.75 |82.3973 |176.15|

| | | |9 |

|7|5.25 |9.044 |185.20|

| | | |3 |

Подсчет s_res

|( f(x)|s_res = h / 3 * (s_ab + 2 * s_even +|Абсолютная ошибка |

|dx |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)