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

В MATLAB реализованы множество современных методов численного интегрирования. Мы рассмотрим простейшие из них.

А. Метод трапеций. Метод трапеции реализован в MATLAB несколькими функциями:

cumtrapz(у)
cumtrapz(x,y)
trapz(y)
trapz(x,y)

Наиболее интересна последняя из них. Данная функция вычисляет интеграл от функции у(х) по х методом трапеций. Аргумент и функция задаются в виде векторов или х — в виде вектора, а у — в виде матрицы. Если у(х) — матрица, то функция возвращает вектор значений интеграла для каждого столбца матрицы.

Пример.Пусть подынтегральная функция имеет вид

у(х) = х*е^x + ln(x) + 1

Необходимо вычислить определенный интеграл в диапазоне от 1 до 10 с шагом 0.5.
Решение:

x=1:0.5:10;
y=x.*exp(x)+log(x)+1;
trapz(x,y)

Решение:

ans =
2.032841320958599e+005

Б. Метод Симпсона. Метод Симпсона реализован в MATLAB также несколькими функциями. Мы рассмотрим простейшую из них:

quad('fun', a, b)
quad('fun', a, b, tol)

Обозначения:

'fun'— подынтегральная функция, взятая в одинарные кавычки;
а, b — пределы интегрирования;
tol — относительная погрешность, задаваемая пользователем; по умолчанию tol=10^-3.

Пример. Вычислить интеграл от функции x+e^x на отрезке [1, 2] с точностью 10^-5.

quad('x+exp(x)', 1, 2, 1e-5)

Результат:

ans = 6.170774280645000e+000

Правило Рунге оценки погрешности интегрирования

В формулах для оценки погрешности квадратурных формул R используются значения производных подинтегральной функции, что требует дополнительного анализа и вычислений. В связи с этим получило распространение практическое правило Рунге оценки погрешности.

Пусть

I – точное значение интеграла,
I(n) – значение интеграла вычисленное при n узлах интегрирования h = (b-a)/n,
I(2n) – значение интеграла вычисленное при 2*n узлах интегрирования, h = (b-a)/2n.

Необходимо определить, с какой точностью вычислен итеграл I(2n), т.е. найти абсолютную погрешность

Для непосредственно определения данной погрешности необходимо найти максимум модуля соответствующей производжной от интегрируемой функции на отрезке [a, b]. Часто это достаточно трудоемкий или вообще невозможный процесс. Напрмер если интегрируемая функция задана таблично. В таких случаях оценку погрешности величины I(2n) можно провести следующим образом:

Здесь m = 3 для методов средних прямоугольников и трапеций, m = 15 для метода Симпсона.

Примечание. Если решается задача численного вычисления интеграла с заданной точностью, процесс удвоения числа узлов интегрирования продолжается до тех пор, пока величена не станет меньше заданной погрешности.