Программирование и комп-ры

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


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



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

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


                                  Выполнил:
                         студент ф – та ЭОУС – 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 
                 #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)




смотреть на рефераты похожие на "Метод Симпсона на компьютере "