Численное интегрирование

Метод прямоугольников
Одним из простейших методов численного интегрирования является метод прямоугольников. На частичном отрезке заменяют подынтегральную функцию полиномом Лагранжа нулевого порядка, построенным в одной точке. Естественно в качестве этой точки выбрать среднюю: . Тогда значение интеграла на частичном отрезке:
(1.1)
Подставив это выражение, получим составную формулу средних прямоугольников:
(1.2)
Графическая иллюстрация метода средних прямоугольников представлена на рис.1(a). Из рисунка видно, что площадь криволинейной трапеции приближенно заменяется площадью многоугольника, составленного из N прямоугольников. Таким образом, вычисление определенного интеграла сводится к нахождению суммы N элементарных прямоугольников.
Формулу (1.2) можно представить в ином виде:
или (1.3)
Эти формулы называются формулой левых и правых прямоугольников соответственно. Графически метод левых и правых прямоугольников представлен на рис.1 (б, в). Однако из-за нарушения симметрии в формулах правых и левых прямоугольников, их погрешность значительно больше, чем в методе средних прямоугольников
а) средние прямоугольники
б) левые прямоугольники
в) правые прямоугольники
Рис.1. Интегрирование методом прямоугольников
Метод трапеций
Если на частичном отрезке подынтегральную функцию заменить полиномом Лагранжа первой степени, то есть:
(1.4)
то искомый интеграл на частичном отрезке запишется следующим образом:
(1.5)
И, составная формула трапеций на всем отрезке интегрирования примет вид:
(1.6)
Графически метод трапеций представлен на рис.2.3. Площадь криволинейной трапеции заменяется площадью многоугольника, составленного из Nтрапеций, при этом кривая заменяется вписанной в нее ломаной. На каждом из частичных отрезков функция аппроксимируется прямой, проходящей через конечные значения, при этом площадь трапеции на каждом отрезке определяется по формуле 1.5.
Погрешность метода трапеций выше, чем у метода средних прямоугольников. Однако на практике найти среднее значение на элементарном интервале можно только у функций, заданных аналитически (а не таблично), поэтому использовать метод средних прямоугольников удается далеко не всегда.

Рис.2. Интегрирование методом методом трапеций
Метод Симпсона
В этом методе подынтегральная функция на частичном отрезке аппроксимируется параболой, проходящей через три точки ,,, то есть интерполяционным многочленом Лагранжа второй степени:
(1.7)
Проведя интегрирование, получим:
(1.8)
Это и есть формула Симпсона или формула парабол. На отрезке формула Симпсона примет вид:
(1.9)
Можно разбить отрезок интегрирования на четное количество 2N равных частей с шагом . Тогда можно построить параболу на каждом сдвоенном частичном отрезке и переписать выражения (1.7-1.9) без дробных индексов. Тогда формула Симпсона примет вид:
(1.10)
Графическое представление метода Симпсона показано на рис.3. На каждом из сдвоенных частичных отрезков заменяем дугу данной кривой параболой.

Рис.3. Метод Симпсона

Реализация
Приведем пример реализации численного интегрирования методом левых и правых прямоугольников. На вход алгоритму подается набор точек, по которым требуется найти приближенное значение интеграла неизвестной функции. На выходе алгоритм выдает найденное приближенное значение с 8-ю знаками точности. Заметим, что эти два метода практически идентичны, единственное отличие заключается в том, ординату какой точки взять в качестве высоты прямоугольника. Все вычислительные операции производятся с типом данных long double в целях повышения точности вычислений.
1 #include 2 #include 3 using namespace std;
4 // Для удобства хранения заданных точек
5 // создадим соответствующую структуру
6 struct point
7 {
8 long double x, y;
9 };
10 int main()
11 {
12 // Объявляем и считываем число точек,
13 // по которым будем вычислять приближенное
14 // значение интеграла функции
15 int pointsCount;
16 cin >> pointsCount;
17 // Точки будем хранить в векторе структур;
18 // его размер, очевидно, равен pointsCount
19 vector points;
20 points.resize (pointsCount);
21 // Считываем абсциссы и ординаты точек
22 for (int i = 0; i < pointsCount; i++)
23 {
24 cin >> points[i].x >> points[i].y;
25 }
26 // Изначально приравниваем приближенное
27 // значение интеграла к нулю
28 long double integralValue = 0.0;
29 // Для каждой пары соседних точек считаем
30 // площадь левого или правого прямоугольника,
31 // который они образуют вместе с осью абсцисс,
32 // по соответствующим теоретическим формулам
33 for (int i = 1; i < pointsCount; i++)
34 {
35 // Воспользуемся формулой левых прямоугольников
36 integralValue += (points[i].x — points[i — 1].x) * points[i — 1].y;
37 // Формула же для правых прямоугольников имеет вид
38 // integralValue += (points[i].x — points[i — 1].x) * points[i].y;
39 }
40 // Выводим приближенное значение интеграла
41 // c восемью знаками точности
42 printf ("%.8llf", integralValue);
43 return 0;
44 }
Приведем пример реализации численного интегрирования методом трапеций. На вход алгоритму подается набор точек, по которым требуется найти приближенное значение интеграла неизвестной функции. На выходе алгоритм выдает найденное приближенное значение с 8-ю знаками точности. Все вычислительные операции производятся с типом данных long double в целях повышения точности вычислений.
1 #include 2 #include 3 using namespace std;
4 // Для удобства хранения заданных точек
5 // создадим соответствующую структуру
6 struct point
7 {
8 long double x, y;
9 };
10 int main()
11 {
12 // Объявляем и считываем число точек,
13 // по которым будем вычислять приближенное
14 // значение интеграла функции
15 int pointsCount;
16 cin >> pointsCount;
17 // Точки будем хранить в векторе структур;
18 // его размер, очевидно, равен pointsCount
19 vector points;
20 points.resize (pointsCount);
21 // Считываем абсциссы и ординаты точек
22 for (int i = 0; i < pointsCount; i++)
23 {
24 cin >> points[i].x >> points[i].y;
25 }
26 // Изначально приравниваем приближенное
27 // значение интеграла к нулю
28 long double integralValue = 0.0;
29 // Для каждой пары соседних точек считаем
30 // площадь трапеции, которую они образуют
31 // вместе с осью абсцисс, по соответствующим
32 // теоретическим формулам
33 for (int i = 1; i < pointsCount; i++)
34 {
35 integralValue += (points[i].x — points[i — 1].x) * (points[i].y + points[i — 1].y);
36 }
37 // Для небольшого ускорения работы алгоритма
38 // деление пополам выносят за пределы цикла
39 integralValue /= 2.0;
40 // Выводим приближенное значение интеграла
41 // c восемью знаками точности
42 printf ("%.8llf", integralValue);
43 return 0;
44 }

Программа на С++ для интегрирования методом Симпсона приведена ниже:

1 #include 2 #include 3 using namespace std;
5 // f(x) — функция для численного интегрирования
6 double f(double x)
7 {
8 return x; // Сюда можно написать любую функцию
9 }
10 /* Интеграл f(x) от a до b:
11 N — шаг интегрирования
12 double integrate(double a, double b, int N)
13 {
14 double h = (b — a) / N;
15 double sum = 0.5 * h * (f(a) + f(b));
16 for (int k = 0; k < N; k++)
17 sum += h * f(a + h*k);
18 return sum;
19 }
20 // Рекурсивный вызов функции для интегрирования.
21 // a, b — пределы интегрирования
22 double area_under_curve(double a,double b)
23 {
24 double area = integrate(a,b,1000);
25 double check = integrate(a,b,5);
26 if (abs(area — check) > 0.00001)
27 {
28 double m = (a + b) / 2;
28 area = area_under_curve(a, m) + area_under_curve(m, b);
29 }
30 return area;
31 }
32 int main()
33 {
34 cout<<«Интеграл x.dx от 1 до 2 равен „<<integrate(1,2,1000)<<endl;
35 cout<<“Площадь под кривой y=x от 1 до 2 равна „<<area_under_curve(1,2);
36 return 0;
37 }

0 комментариев

Только зарегистрированные и авторизованные пользователи могут оставлять комментарии.