Численное интегрирование методом Симпсона (реализация на C++)

C++   26 Февраль 2012  Автор статьи:  

Приведем пример реализации численного интегрирования методом Симпсона. На вход алгоритму подается набор точек, по которым требуется найти приближенное значение интеграла неизвестной функции. На выходе алгоритм выдает найденное приближенное значение с 8-ю знаками точности. В данном методе интерполирование неизвестной функции осуществляется многочленом второй степени, т.е. на каждом отрезке функция приближается параболой.

Сложность алгоритма — O(n), где n — число точек, по которым осуществляется приближенное интегрирование. Все вычислительные операции производятся с типом данных long double в целях повышения точности вычислений.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#include <iostream>
#include <vector>
using namespace std;

// Для удобства хранения заданных точек
// создадим соответствующую структуру
struct point
{
    long double x, y;  
};

int main()
{
    // Объявляем и считываем число точек,
    // по которым будем вычислять приближенное
    // значение интеграла функции
    int pointsCount;
    cin >> pointsCount;

    // Точки будем хранить в векторе структур;
    // его размер, очевидно, равен pointsCount
    vector <point> points;
    points.resize (pointsCount);

    // Считываем абсциссы и ординаты точек
    for (int i = 0; i < pointsCount; i++)
    {
        cin >> points[i].x >> points[i].y;
    }  

    // Изначально приравниваем приближенное
    // значение интеграла к нулю
    long double integralValue = 0.0;

    // Для троек соседних точек считаем
    // площадь под графиком параболы,
    // которую они образуют, по
    // соответствующим теоретическим формулам;
    // заметим, что в этом случае указатель i
    // на каждой итерации сдвигается на 2 ед.
    for (int i = 2; i < pointsCount; i += 2)
    {
        long double h = points[i].x - points[i - 2].x;
        integralValue += (points[i - 2].y + 4 * points[i - 1].y + points[i].y) * h;
    }

    // Для небольшого ускорения работы алгоритма
    // деление на 6 выносят за пределы цикла 
    integralValue /= 6.0;
   
    // Выводим приближенное значение интеграла
    // c восемью знаками точности
    printf ("%.8llf", integralValue);

    return 0;
}

Научиться программировать

  • на Delphi

  • на Java

  • на C++