Методы численного интегрирования в MathCAD
Теорию по численному интегрированию можно почитать, например, здесь, а в этой заметке займёмся реализацией в Маткаде основных методов численного интегрирования, которые чаще всего "проходят" в ВУЗах.
Все рассмотренные ниже методы, в сущности, между собой похожи - если одномерный определённый интеграл есть площадь криволинейной трапеции под графиком:
Одномерный определённый интеграл
, то весь вопрос только в том, какой именно из простых зависимостей (прямая, парабола и т.п.) мы заменим подынтегральную функцию, от которой, в общем случае, интеграл не берётся аналитически (или которая нам неизвестна, но приближена интерполяционным полиномом, или интеграл можно взять, но очень трудоёмко и т.д.)
Ясно, что можно заменить и вот так:
Простейшее "интегрирование" - интеграл как площадь прямоугольника :)
, считая, что площадь жирного прямоугольника приблизительно равна искомой площади под кривой, но это будет очень уж неточно, поэтому отрезок интегрирования по оси x всегда разбивают на небольшие интервалы (проще всего, с постоянным шагом h) и находят значение интеграла как сумму площадей простых фигур, например, прямоугольников, нижняя сторона которых равна h, а высота - значению f(x), взятому в некоторой точке интервала (на рисунке - в серединах):
Простейший "метод прямоугольников"
Ясно, что погрешность уменьшится, но останется.
Теперь от слов к Маткаду. Определим тестовую функцию f(x), пределы интегрирования [a,b] и число интервалов n, на которое разбивается отрезок [a,b]. Величину шага h затем вычислим как (b-a)/n. В учебных целях выведем также "точное" значение искомого интеграла. Следует понимать, что "точное" оно лишь в кавычках, MathCAD-то искал его тоже численным методом.
Численное интегрирование: определение тестовых данных
Реализуем три основных метода прямоугольников. Разница между ними в том, в какой точке каждого отрезка на интервале интегрирования - левой, правой или в середине - берётся значение функции f(x).
Интегрирование методами прямоугольников в MathCAD
В методе трапеций мы для каждого отрезка интегрирования [xi,xi+1] соединяем отрезком прямой линии точки f(xi) и f(xi+1), считая интеграл как сумму площадей трапеций. Это всегда точнее, а сам метод ещё достаточно прост. По-моему, близок к оптимуму при массовых расчётах.
Метод трапеций в MathCAD
Наконец, в методе Симпсона (парабол) функцию f(x) на каждом отрезке интегрирования заменяют параболой, то есть, кривой второго порядка. Расчёт становится сложнее, но точность повышается в разы. Существует немало разновидностей формулы для метода Симпсона, вот 2 неплохих способа расчёта:
Метод Симпсона (парабол) в MathCAD
Ниже показаны оценки погрешностей для всех методов.
Оценки погрешностей всех методов при n=10
Увеличивая число интервалов n, можно оценить и порядок точности всех методов.
Например, для метода первого порядка точности (методы левых и правых прямоугольников) при увеличении числа интервалов разбиения по оси x вдвое (n:=20 вместо n:=10 в начале документа) погрешность решения должна уменьшиться примерно в 2 раза. Для методов второго порядка точности (средних прямоугольников, трапеций) при уменьшении шага по x вдвое погрешность уменьшится примерно в 4 раза (второй по h порядок точности и означает, что погрешность уменьшается пропорционально величине h2). Метод Симпсона имеет четвёртый порядок точности, то есть, при уменьшении шага вдвое (увеличении вдвое числа интервалов n) погрешность решения уменьшится примерно в 24=16 раз.
Следует помнить, что на дискретизации по оси x свет клином не сошёлся, существуют красивые альтернативные методы, скажем, метод Монте-Карло :) При многомерном интегрировании он становится, пожалуй, предпочтительней.
Скачать документ "Численные методы интегрирования" (.xmcd, MathCAD 15) (93 Кб)
Показанные методы можно реализовать и без использования инструментов панели программирования, только с помощью оператора суммы и арифметики. Приведу примеры для методов средних прямоугольников, трапеций и Симпсона, вроде бы, всё работает:
Методы средних прямоугольников, трапеций и Симпсона - реализация без модульного программирования
30.10.2013, 18:11 [79559 просмотров]