Реализация некоторых численных методов
8.2 Численное интегрирование
Задача численного интегрирования состоит в замене исходной подинтегральной функции , для которой трудно или невозможно записать первообразную в аналитике, некоторой аппроксимирующей функцией
. Такой функцией обычно является полином (кусочный полином)
.
Таким образом



8.2.1 Обзор методов интегрирования
Методы вычисления однократных интегралов называются квадратурными (для кратных интегралов — кубатурными), и делятся на следующие группы:
- Методы Ньютона-Котеса. Здесь
— полином различных степеней. Сюда относятся метод прямоугольников, трапеций, Симпсона.
- Методы статистических испытаний (методы Монте-Карло). Здесь узлы сетки для квадратурного или кубатурного интегрирования выбираются с помощью датчика случайных чисел, ответ носит вероятностный характер. В основном применяются для вычисления кратных интегралов.
- Сплайновые методы. Здесь
— кусочный полином с условиями связи между отдельными полиномами посредством системы коэффициентов.
- Методы наивысшей алгебраической точности. Обеспечивают оптимальную расстановку узлов сетки интегрирования и выбор весовых коэффициентов
в задаче
(характерный пример — метод Гаусса).
8.2.2 Метод прямоугольников
Различают метод левых, правых и средних прямоугольников. Суть метода ясна из рисунка. На каждом шаге интегрирования функция аппроксимируется полиномом нулевой степени — отрезком, параллельным оси абсцисс.
Формулы метода прямоугольников можно получить из анализа разложения функции в ряд Тейлора вблизи некоторой точки
:

Рассмотрим диапазон интегрирования от до
, где
— шаг интегрирования. Вычислим интеграл от исследуемой функции на этом промежутке:






8.2.3 Метод средних прямоугольников
Здесь на каждом интервале значение функции считается в средней точке отрезка , то есть
.
Разложение функции в ряд Тейлора показывает, что в случае средних прямоугольников точность метода существенно выше:

Пример функции Maxima, реализующей метод средних прямоугольников, представлен ниже:
(%i1) intpr(f,n,a,b):=block([h,i,s,_x],h:(b-a)/n, _x:a+h/2, s:0, for i:1 thru n do (s:s+float(subst(_x,x,f)),_x:_x+h),s:s*h)$ (%i2) intpr(x^2,100,-1,1);

8.2.4 Метод трапеций
Аппроксимация в этом методе осуществляется полиномом первой степени. На единичном интервале

В случае равномерной сетки ( = const)

При этом , а
.
Погрешность метода трапеций в два раза выше, чем у метода средних прямоугольников. Однако на практике найти среднее значение на элементарном интервале можно только у функций, заданных аналитически (а не таблично), поэтому использовать метод средних прямоугольников удаётся далеко не всегда. В силу разных знаков погрешности в формулах трапеций и средних прямоугольников истинное значение интеграла обычно лежит между двумя этими оценками.
Пример реализации метода трапеций в виде функции приведен ниже:
(%i1) f:x^2;

(%i2) inttrap(f,n,a,b):=block([h,i,s],h:(b-a)/n, s:(float(subst(a,x,f))+float(subst(b,x,f)))/2, for i:1 thru n-1 do (s:s+float(subst(a+i*h,x,f))), s:s*h)$ (%i3) inttrap(f,100,-1,1);

Подинтегральная функция задаётся в виде выражения Maxima. Выражение, определяющее подинтегральную функцию, можно задавать и непосредственно при обращении к методу, как в следующем примере:
(%i4) inttrap(x*sin(x),100,-1,1);

8.2.5 Метод Симпсона
При использовании данного метода подинтегральная функция заменяется интерполяционным полиномом второй степени
— параболой, проходящей через три соседних узла. Рассмотрим два шага интегрирования (
), то есть три узла
, через которые проведем параболу, воспользовавшись уравнением Ньютона:

Пусть , тогда

Воспользовавшись полученным соотношением, вычислим интеграл по данному интервалу:

В итоге .
Для равномерной сетки и чётного числа шагов n формула Симпсона принимает вид:



Реализация метода Симпсона средствами Maxima представлена в следующем примере:
(%i1) intsimp(f,n,a,b):=block([h,h2,i,_x,s],h:(b-a)/n, h2:h/2, s:(float(subst(a,x,f))+float(subst(b,x,f)))/2+ 2*float(subst(a+h2,x,f)), _x:a, for i:1 thru n-1 do (_x:_x+h, s:s+2*float(subst(_x+h2,x,f))+float(subst(_x,x,f))), s:s*h/3)$ (%i2) intsimp(x^2,100,-1,1);
