Методы численного интегрирования
К вычислениям определенных интегралов сводятся многие практические задачи физики, химии, экологии, механики и других естественных наук. На практике взять интеграл аналитически не всегда удается. В этом случае используются методы численного интегрирования. В данной лабораторной работе рассматриваются методы Ньютона – Котеса, в частности методы прямоугольников, трапеций, Симпсона и метод Гаусса. Кроме того, в лабораторной работе рассматриваются способы аналитического и численного отыскания интегралов средствами MATLAB.
Вычислительные формулы для всех рассматриваемых методов приведены в приложении, теоретический материал следует изучать по материалам лекций и рекомендованной литературе. Здесь же будут рассмотрены практические аспекты реализации методов численного интегрирования в MATLAB.
Вам следует внимательно изучить и выполнить приведенные здесь примеры, может быть для функций из ваших вариантов задач.
1. Аналитическое интегрирование средствами MATLAB.
В лабораторной работе, основной задачей которой является исследование методов численного интегрирования, в некоторых случаях для отыскания погрешности результатов требуется точные значения интегралов. Т.е. эти величины необходимо определить аналитически.
Для вычисления определенных и неопределенных интегралов в MATLAB спользуется функция int .
Синтаксис:
Вызов | Описание |
---|---|
int(S) |
Вычисляется неопределенный интеграл от функции S по ее символьной переменной, определенной в syms .
Вычисляется неопределенный интеграл от функции S по ее символьной переменной v , определенной в syms .
Вычисляется определенный интеграл от a до b функции S по ее символьной переменной, определенной в syms . A и b могут быть переменными символьного или вещественного типа.
Пример использования функции int . В данном примере подинтегральная функция задается явным образом.
График подинтегральной функции:
Обратите внимание на форматирование графика прямо в тексте программы.
В следующем примере подинтегральная функция и переменная интегрирования задаются в символьной форме.
Обратите внимание, что в данном примере для построения графика используется m – файл f.m:
Это связано с тем, что аргументом функции fplot не может быть символьная функция.
График интегрируемой функции:
2. Аналитическое дифференцирование средствами MATLAB.
В лабораторной работе, основной задачей которой является исследование методов численного интегрирования, для отыскания теоретической погрешности результатов требуется знание производной от интегрируемой функции до четвертого порядка включительно.
Для вычисления производных аналитически в MATLAB спользуется функция diff .
Синтаксис:
Вызов | Описание |
---|---|
diff(S) | Вычисляется аналитическое выражение для производной от функции S, заданной символьно. |
diff(S,n) | Вычисляется аналитическое выражение для производной от функции S порядка n. |
Пример использования функции diff при вычислении максимума второй производной на отрезке [a, b]
3. Теоретическая оценка погрешности численного интегрирования
Рассмотрим примеры использования теоретической оценки погрешности интегрирования на примере двух задач.
Пример. Определить теоретическую погрешность численного интегрирования методом трапеций в случае одного элементарного отрезка интегрирования..
Теоретическая погрешность для метода трапеций составляет
В случае элементарного отрезка иртегрирования, если заданы пределы интегрирования и подинтегральная функция, задачу можно решить например следующим образом:
В результате получим
M2 =
16.4000
Погрешность
R =
0.1164
Т.о. теоретическая оценка абсолютной погрешности погрешности составляет 0.1 . В задаче 1 лабораторной работы вам предлагается убедиться, что данная теоретическая оценка действительно справедлива.
4. Численное интегрирование
Вычислительные формулы приведены в приложении . Их анализ показывает, что формулы для одного элементарного отрезка интегрирования не требуют каких либо новых знаний. Необходимо только приготовить m-функцию , в которой следует определить интегрируемую функцию. И конечно для оценки погрешности интегрирования следует знать точное значение интеграла (см. выше).
Случай формул для составного отрезка интегрирования более сложный. Здесь требуется вычислять значения сумм. Например в методе Симпсона
В языке программировани MATLAB, как и в других языках программирования, существует оператор цикла for , возможностей которого достаточно для решения поставленных задач.
Здесь для задания интегрируемой функции необходимо создать m – файл f.m
5. Численное интегрирование средствами MATLAB
В MATLAB реализованы множество современных методов численного интегрирования. Мы рассмотрим простейшие из них.
А. Метод трапеций. Метод трапеции реализован в MATLAB несколькими функциями:
Наиболее интересна последняя из них. Данная функция вычисляет интеграл от функции у(х) по х методом трапеций . Аргумент и функция задаются в виде векторов или х — в виде вектора, а у — в виде матрицы. Если у(х) — матрица, то функция возвращает вектор значений интеграла для каждого столбца матрицы.
Пример. Пусть подынтегральная функция имеет вид
у(х) = х*е x + ln(x) + 1
Необходимо вычислить определенный интеграл в диапазоне от 1 до 10 с шагом 0.5.
Решение:
Б. Метод Симпсона. Метод Симпсона реализован в 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)
6. Правило Рунге оценки погрешности интегрирования
В формулах для оценки погрешности квадратурных формул 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 для метода Симпсона.
Примечание. Если решается задача численного вычисления интеграла с заданной точностью, процесс удвоения числа узлов интегрирования продолжается до тех пор, пока величена не станет меньше заданной погрешности.
БлогNot. Методы численного интегрирования в MathCAD
Методы численного интегрирования в MathCAD
Теорию по численному интегрированию можно почитать, например, здесь, а в этой заметке займёмся реализацией в Маткаде основных методов численного интегрирования, которые чаще всего "проходят" в ВУЗах.
Все рассмотренные ниже методы, в сущности, между собой похожи – если одномерный определённый интеграл есть площадь криволинейной трапеции под графиком:
, то весь вопрос только в том, какой именно из простых зависимостей (прямая, парабола и т.п.) мы заменим подынтегральную функцию, от которой, в общем случае, интеграл не берётся аналитически (или которая нам неизвестна, но приближена интерполяционным полиномом, или интеграл можно взять, но очень трудоёмко и т.д.)
Ясно, что можно заменить и вот так:
, считая, что площадь жирного прямоугольника приблизительно равна искомой площади под кривой, но это будет очень уж неточно, поэтому отрезок интегрирования по оси x всегда разбивают на небольшие интервалы (проще всего, с постоянным шагом h ) и находят значение интеграла как сумму площадей простых фигур, например, прямоугольников, нижняя сторона которых равна h , а высота – значению f(x) , взятому в некоторой точке интервала (на рисунке – в серединах):
Ясно, что погрешность уменьшится, но останется.
Теперь от слов к Маткаду. Определим тестовую функцию f(x) , пределы интегрирования [a,b] и число интервалов n , на которое разбивается отрезок [a,b] . Величину шага h затем вычислим как (b-a)/n . В учебных целях выведем также "точное" значение искомого интеграла. Следует понимать, что "точное" оно лишь в кавычках, MathCAD-то искал его тоже численным методом.
Реализуем три основных метода прямоугольников. Разница между ними в том, в какой точке каждого отрезка на интервале интегрирования – левой, правой или в середине – берётся значение функции f(x) .
В методе трапеций мы для каждого отрезка интегрирования [xi,xi+1] соединяем отрезком прямой линии точки f(xi) и f(xi+1) , считая интеграл как сумму площадей трапеций. Это всегда точнее, а сам метод ещё достаточно прост. По-моему, близок к оптимуму при массовых расчётах.
Наконец, в методе Симпсона (парабол) функцию f(x) на каждом отрезке интегрирования заменяют параболой, то есть, кривой второго порядка. Расчёт становится сложнее, но точность повышается в разы. Существует немало разновидностей формулы для метода Симпсона, вот 2 неплохих способа расчёта:
Ниже показаны оценки погрешностей для всех методов.
Увеличивая число интервалов n , можно оценить и порядок точности всех методов.
Например, для метода первого порядка точности (методы левых и правых прямоугольников) при увеличении числа интервалов разбиения по оси x вдвое ( n:=20 вместо n:=10 в начале документа) погрешность решения должна уменьшиться примерно в 2 раза. Для методов второго порядка точности (средних прямоугольников, трапеций) при уменьшении шага по x вдвое погрешность уменьшится примерно в 4 раза (второй по h порядок точности и означает, что погрешность уменьшается пропорционально величине h 2 ). Метод Симпсона имеет четвёртый порядок точности, то есть, при уменьшении шага вдвое (увеличении вдвое числа интервалов n ) погрешность решения уменьшится примерно в 2 4 =16 раз.
Следует помнить, что на дискретизации по оси x свет клином не сошёлся, существуют красивые альтернативные методы, скажем, метод Монте-Карло 🙂 При многомерном интегрировании он становится, пожалуй, предпочтительней.
Скачать документ "Численные методы интегрирования" (.xmcd, MathCAD 15) (93 Кб)
Показанные методы можно реализовать и без использования инструментов панели программирования, только с помощью оператора суммы и арифметики. Приведу примеры для методов средних прямоугольников, трапеций и Симпсона, вроде бы, всё работает:
30.10.2013, 18:11; рейтинг: 50359
Доброго времени суток! Мы продолжаем говорить о численных методах. И сегодня мы поговорим о реализации численных методов интегрирования в среде Matlab.
Численное интегрирование в Matlab
Геометрический смысл интегрирования — это нахождение площади, которая находится под интегрируемой функцией. На рисунке показана площадь для определённого интеграла, ограниченного a и b.
Численное интегрирование не только в Matlab, но и в других средах, строится именно на нахождении площади. Для начала мы разберем простые методы:
Методы прямоугольников
- метод правых прямоугольников
- метод левых прямоугольников
- метод средних прямоугольников
Суть их в построение под кривой прямоугольников одинаковый ширины и нахождение их суммарной площади.Как видите, они различаются только точкой соприкосновения с кривой. Методы достаточны простые в реализации. Однако, погрешности данных методов весьма высоки. Точнее говоря, методы прямоугольников имеют первый порядок точности. Это означает, что ошибка пропорциональна шагу и накапливается со временем. Соответственно, чем меньше шаг, тем меньшую ошибку мы получим.
Также, следует отметить, что метод средних прямоугольников является более точным и предпочтительно использовать именно этот метод численного интегрирования, если у вас стоит выбор из этих трех методов. Эту точность можно доказать с помощью разложения в ряд Тейлора.
Необходимо посчитать интеграл функции f(x) = xe sin(x) x с шагом разбиения h = 0.02 на интервале от 0 до 1.
Функция feval (родственник функции eval) — интерпретирует и вычисляет текстовую строку, которая может содержать либо арифметическое выражение, либо инструкцию, либо обращение к функции, однако, в отличии от eval, интерпретирует и вычисляет текстовую строку, которая может содержать либо арифметическое выражение, либо инструкцию, либо обращение к функции.
Метод трапеций
Ещё одни популярный и в тоже время простой метод — метод трапеций. Аналогично методу прямоугольников строятся трапеции под кривой и находится их суммарная площадь. Данный метод имеет второй порядок точности (ошибка пропорциональна шагу в квадрате).
В Matlab метод трапеций реализован двумя функциями:
- cumtrapz()
- trapz()
Первую функцию обычно используют при работе с табличными данными или векторами. Откликом функции является n-интегралов, где n — число элементов вектора или элементов в каждом столбце матрицы. Следующие примеры отображают работу этой функции.
Пусть функция y(x) имеет значения, представленные в виде следующего вектора: y = [1,2,3,4,5,6,7,8,9,10] . Необходимо вычислить:
При этом a = 1; b = 1, 2, 3, 4 …,10.
Пишем в Matlab:
Теперь рассмотрим вариант работы с вектором и матрицей:
Функция y(x) задана в виде матрицы y(x) = [1 3 5; 3 5 7; 4 6 8; 4 7 9; 5 7 10] . При этом аргумент представляет собой вектор: x = [1,3,7,9,10].
Вторая функция для интегрирования, работающая по методу трапеций Matlab — trapz(). Наиболее используемая студентами, так как позволяет работать не только с векторами и матрицами, но и с аналитической формой подынтегральной функции. Выглядит это примерно так:
Необходимо вычислить определённый интеграл в диапазоне от 1 до 10 с шагом 0.5 для заданной функции:
Как видите, ничего сложного. А иногда даже удобнее некоторых онлайн сервисов для расчёта интегралов.
Метод Симпсона
Преимущество этого метода в том, что точки, взятые на каждом шаге на кривой, интерполируются полиномом второй степени. Проще говоря, соединяются параболой. Это даёт методу четвёртый порядок точности.
В Matlab интегрирование с помощью метода Симпсона производит функция quad. Сразу разберем пример.
Вычислить определённый интеграл с точностью 10 -4 методом Симпсона.
Точность вычислений задается 4 параметром функции quad. Также, следует отметить, что в задании нижним пределом является 0, а мы использовали число 0.001. Это связано с тем, что при подстановке 0 функция не определена, а точнее, натуральный логарифм не существует.
Ну и реализация этого метода вручную приведена здесь для общего развития. Этим я хочу подчеркнуть, что практически любой метод или алгоритм возможно написать самому, а не пользоваться стандартными методами Matlab.
Символьное интегрирование в Matlab
Часто нам необходимо найти интеграл от какой либо функции, не зная пределов интегрирования. Тогда нам нужно взять интеграл в общем или символьном виде. В Matlab за символьное интегрирование отвечает функция int. Она принимает как минимум 2 параметра: 1 — функция, 2 — имя переменной по которой берется интеграл. int(fun, var). Рассмотрим короткий пример:
Вычислить неопределённый интеграл:
Следует отметить, что функция int также может считать и определенные интегралы, для этого нужно задать пределы интегрирования в 3 и 4 параметры функции соответственно.
Заключение
На этом я хочу закончить сегодняшнюю тему «Интегрирование в Matlab». Не забывайте, что Matlab позволяет программировать сложные алгоритмы, а не только использовать встроенный функционал. Любой численный метод можно реализовать и вызывать как функцию. Если у вас остались вопросы, то задавайте их в комментариях.
В этот раз без исходников, примеры небольшие.