2.2. Чисельне інтегрування. Метод Ейлера (прямокутників). Метод трапецій. Метод Сімпсона (парабол). Метод Монте-Карло. Метод Філона.
2.3. Нахождение минимума функции многих переменных. Метод градиентного спуска
Чисельне інтегрування
При вирішенні завдань наукового і інженерно-технічного характеру математичними методами часто виникає необхідність проінтегрувати яку-небудь функцію. Є функції, які неможливо інтегрувати аналітично, тобто тільки в деяких випадках по заданій функції можна знайти первинну. Загальним способом інтегрування будь-яких функцій є чисельне інтегрування, методи якого в більшості своїй прості і легко перекладаються на алгоритмічні мови
Чисельне інтегрування (історична назва: квадратура) - обчислення значення визначеного інтеграла (як правило, наближене), засноване на тому, що величина інтеграла чисельно рівна площі криволінійної трапеції, яка обмежена віссю абсцис, графіком інтегруємої функції і відрізками прямих x = а і x = b, де а і b - межі інтегрування (див. рис.).
HYPERLINK "http://ru.wikipedia.org/wiki/%D0%A4%D0%B0%D0%B9%D0%BB:Integral_as_region_under_curve.svg" \o "Определённый интеграл как площадь фигуры" INCLUDEPICTURE "http://upload.wikimedia.org/wikipedia/commons/thumb/f/f2/Integral_as_region_under_curve.svg/200px-Integral_as_region_under_curve.svg.png" \* MERGEFORMATINET
Рис. . Визначений інтеграл як площа фігури.
Чисельні методи інтегрування використовують заміну площі криволінійної трапеції на кінцеву суму площ простіших геометричних фігур, які можуть бути обчислені точно. У цьому сенсі говорять про використання квадратурних формул (по аналогії із завданням про квадратуру круга - побудова квадрата з площею, рівній площі круга з певним радіусом).
У більшості методів використовується наближеною представлення інтеграла у вигляді кінцевої суми (квадратурна формула):
де ci - постійні, звані вагами, а xi - належать інтервалу [а, b] і називаються вузлами.
В основі квадратурних формул лежить ідея апроксимації на відрізку інтегрування графіка підінтегрального виразу функціями простішого вигляду, які легко можуть проінтегрувати аналітично і, таким чином, легко обчислені.
http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integration.html
Формула прямоугольников (метод Ейлера).
.
Интегрирование методом прямоугольников (метод Эйлера).
Пусть функцию (рисунок справа ) необходимо проинтегрировать численным методом на отрезке [a, b]. Разделим отрезок на N равных интервалов (на рисунке N = 4).
Площадь каждой из 4-х криволинейных трапеций можно заменить на площадь прямоугольника.
Ширина всех прямоугольников одинакова и равна INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image029.gif" \* MERGEFORMATINET
В качестве выбора высоты прямоугольников можно предложить выбрать значение функции на левой границе. В этом случае высота первого прямоугольника составит f(a), второго – f(x1), третьего – f(x2), последнего – f(x3). Приближенное значение интеграла получается суммированием площадей прямоугольников
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image031.gif" \* MERGEFORMATINET
Если в качестве выбора высоты прямоугольников взять значение функции на правой границе, то в этом случае высота первого прямоугольника составит f(x1), второго – f(x2), третьего – f(x3), последнего – f(b).
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image033.gif" \* MERGEFORMATINET
Как видно, в этом случае, одна из формул дает приближение к интегралу с избытком, а вторая c недостатком. Можно предложить еще один способ, очевидно лучший, чем обе эти формулы – использовать для аппроксимации значение функции в середине отрезка интегрирования.
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image035.gif" \* MERGEFORMATINET
В общем виде, если отрезок [a, b] разбить на N равных интервалов интегрирования (h) и к каждому интервалу применить формулу прямоугольников, то получим
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image037.gif" \* MERGEFORMATINET
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image039.gif" \* MERGEFORMATINET INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image041.gif" \* MERGEFORMATINET INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image043.gif" \* MERGEFORMATINET
Формула трапеций
Использование для интерполяции полинома первой степени (прямая линия, проведенная через две точки) приводит к формуле трапеций. В качестве узлов интерполирования берутся концы отрезка интегрирования. Таким образом, криволинейная трапеция заменяется на обычную трапецию, площадь которой может быть найдена как произведение полусуммы оснований на высоту.
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image045.gif" \* MERGEFORMATINET
В случае N отрезков интегрирования для всех узлов, за исключением крайних точек отрезка, значение функции войдет в общую сумму дважды (так как соседние трапеции имеют одну общую сторону).
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image047.gif" \* MERGEFORMATINET INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image049.gif" \* MERGEFORMATINET
Интересно, что формула трапеций может быть получена, если взять половину суммы формул прямоугольников по правому и левому краям отрезка
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image051.gif" \* MERGEFORMATINET INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image053.gif" \* MERGEFORMATINET
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image055.gif" \* MERGEFORMATINET
Проиллюстрируем использование формулы трапеций на примере рисунка 1
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image057.gif" \* MERGEFORMATINET
Величину I можно представить как сумму площадей трапеций (в данном случае четырех)
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image059.gif" \* MERGEFORMATINET
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image061.gif" \* MERGEFORMATINET
Проверка устойчивости решения
Как правило, чем меньше длина каждого интервала, т.е. чем больше число этих интервалов, тем меньше различаются приближенное и точное значение интеграла. Это справедливо для большинства функций. В методе трапеций ошибка вычисления интеграла (δ) приблизительно пропорциональна квадрату шага интегрирования h 2
δ ~ h 2
Таким образом, для вычисления интеграла некоторой функции в пределах a, b необходимо разделить отрезок [a, b] на n0 интервалов и найти сумму площадей трапеций. Затем нужно увеличить число интервалов (n1), опять вычислить сумму трапеций и сравнить полученное значение с предыдущим результатом. Это следует повторять до тех пор (ni), пока не будет достигнута заданная точность результата (критерия сходимости).
Для методов прямоугольников и трапеций обычно на каждом шаге итерации число интервалов увеличивается в 2 раза, т.е. ni+1 = 2 ni. Алгоритм процедуры интегрирования можно записать следующим образом:
интеграл (I) рассчитывается по формуле
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image063.gif" \* MERGEFORMATINET , где INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image065.gif" \* MERGEFORMATINET ,
а критерий сходимости
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image067.gif" \* MERGEFORMATINET
Главное преимущество правила трапеций – его простота. Однако если при вычислении интеграла требуется высокая точность, применение этого метода может потребовать слишком большого количества итераций или машинного времени.
Пример:
Пользуясь правилом трапеций вычислить интеграл INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image069.gif" \* MERGEFORMATINET . (Точное решение 1/3)
Для n = 1 INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image071.gif" \* MERGEFORMATINET
Для n = 2 INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image073.gif" \* MERGEFORMATINET
Для n = 4 INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image075.gif" \* MERGEFORMATINET
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image077.gif" \* MERGEFORMATINET
Для n = 64 INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image079.gif" \* MERGEFORMATINET
Правило Симпсона
Использование трех точек для интерполирования подынтегрального выражения позволяет использовать параболическую функцию (полином второй степени). Это приводит к формуле Симпсона приближенного вычисления интеграла.
Рассмотрим произвольный интеграл
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image081.gif" \* MERGEFORMATINET
Воспользуемся заменой переменной таким образом, чтобы границы отрезка интегрирования вместо [a,b] стали [-1,1], для этого введем переменную z:
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image083.gif" \* MERGEFORMATINET , Тогда INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image085.gif" \* MERGEFORMATINET и INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image087.gif" \* MERGEFORMATINET
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image089.gif" \* MERGEFORMATINET
Рассмотрим задачу интерполирования полиномом второй степени (параболой) подынтегральной функции, используя в качестве узлов три равноудаленные узловые точки – z = -1, z = 0, z = +1 (шаг равен 1, длина отрезка интегрирования равна 2). Обозначим соответствующие значения подынтегральной функции в узлах интерполяции
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image091.gif" \* MERGEFORMATINET INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image093.gif" \* MERGEFORMATINET INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image095.gif" \* MERGEFORMATINET
Система уравнений для нахождения коэффициентов полинома
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image097.gif" \* MERGEFORMATINET , проходящего через три точки INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image099.gif" \* MERGEFORMATINET , INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image101.gif" \* MERGEFORMATINET и INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image103.gif" \* MERGEFORMATINET
Примет вид
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image105.gif" \* MERGEFORMATINET или INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image107.gif" \* MERGEFORMATINET
Коэффициенты легко могут быть получены
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image109.gif" \* MERGEFORMATINET
Вычислим теперь значение интеграла от интерполяционного многочлена
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image111.gif" \* MERGEFORMATINET
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image113.gif" \* MERGEFORMATINET
Путем обратной замены переменной вернемся к исходному интегралу. Учтем, что
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image115.gif" \* MERGEFORMATINET соответствует INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image117.gif" \* MERGEFORMATINET
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image119.gif" \* MERGEFORMATINET соответствует INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image121.gif" \* MERGEFORMATINET
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image123.gif" \* MERGEFORMATINET соответствует INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image125.gif" \* MERGEFORMATINET
Получим формулу Симпсона для произвольного интервала интегрирования:
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image127.gif" \* MERGEFORMATINET , и INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image129.gif" \* MERGEFORMATINET
При необходимости, исходный отрезок интегрирования может быть разбит на N сдвоенных отрезков, к каждому из которых применяется формула Симпсона. Шаг интерполирования при этом составит
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image131.gif" \* MERGEFORMATINET
Для первого отрезка интегрирования узлами интерполирования будут являться точки a, a+h, a+2h, для второго – a+2h, a+3h, a+4h, третьего a+4h, a+5h, a+6h и т.д. Приближенное значение интеграла получается суммированием N площадей:
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image133.gif" \* MERGEFORMATINET
В данную сумму входят одинаковые слагаемые (для внутренних узлов с четным значением индекса - 2i). Поэтому можно перегруппировать слагаемые в этой сумме таким образом
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image135.gif" \* MERGEFORMATINET , что эквивалентно
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image137.gif" \* MERGEFORMATINET , так как INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image139.gif" \* MERGEFORMATINET
Погрешность этого приближенного метода уменьшается пропорционально длине шага интегрирования в четвертой степени, т.е. при увеличении числа интервалов вдвое ошибка уменьшается в 16 раз
δ ~ h 4
Пример:
Пользуясь правилом Симпсона вычислить интеграл INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image141.gif" \* MERGEFORMATINET . (Точное решение - 0,2)
Для n = 1
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image143.gif" \* MERGEFORMATINET
Для n = 2
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image145.gif" \* MERGEFORMATINET
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image147.gif" \* MERGEFORMATINET
Правило Симпсона позволяет точно рассчитать интеграл не только от квадратичной функции, но и для полинома третей степени
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image149.gif" \* MERGEFORMATINET
Метод Монте-Карло
В некоторых случаях из-за особенности подынтегральной функции (например, из-за ее большой сложности, неявном способе задания и т.д.), описанные выше методы нельзя или нецелесообразно использовать. В задачах, не требующих высокой точности, широкое распространение получил метод Монте-Карло.
Проиллюстрируем применение этого метода на примере приближенного вычисления следующего интеграла:
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image255.gif" \* MERGEFORMATINET
График подынтегрального выражения приведен на рисунке. Очевидно, что точное значение интеграла равно четверти площади круга единичного радиуса.
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image256.gif" \* MERGEFORMATINET
Иллюстрация метода Монте-Карло
Построим прямоугольную область, которая будет полностью включать в себя искомый интеграл. В данном случае это будет квадрат с единичным ребром, показанный на рисунке. Далее, с помощью датчика случайных чисел генерируются точки
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image258.gif" \* MERGEFORMATINET ,
попадающие в эту область. В данном случае абсциссы и ординаты точек должны быть случайными числами, равномерно распределенными на отрезке [0, 1].
Для каждой точки проверяется, попадает ли она в область под или над графиком функции, то есть проверяется условие:
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image260.gif" \* MERGEFORMATINET
Если условие выполняется, то выбранная точка соответствует «успеху», если нет – то «промаху». Таким образом, процедура может быть описана как игра в «попадание» случайно выбранной точки в область под графиком (отсюда и название метода - Монте-Карло).
Вполне очевидно, что отношение числа «попаданий» (Nусп) к общему числу попыток (N) должно в пределе стремиться к доли площади прямоугольной области (Sпр), которую занимает область под интегрируемой функцией (значение интеграла, I).
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image262.gif" \* MERGEFORMATINET
Отсюда получается формула метода Монте-Карло:
INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image264.gif" \* MERGEFORMATINET
Для реализации метода существенное значение имеет качество используемого датчика случайных чисел. Идеальный датчик должен давать равномерное распределение чисел в заданном диапазоне. Точность расчета интеграла определяется так же числом точек (N), используемых при вычислениях и, очевидно, должна увеличиваться при его росте.
Метод Монте-Карло широко используется в современных методах моделирования динамики молекулярных систем, взаимодействия растворенного вещества с молекулами растворителя, кинетики адсорбции веществ на твердых поверхностях и т.д.
Пример. Вычислить объем молекулы, ограниченный Ван-дер-Ваальсовыми радиусами образующих ее атомов. Очевидно, что из рассмотренных методов единственно возможным является метод Монте-Карло. Подынтегральная функция в этом случае не имеет аналитического представления, зато достаточно легко можно проверить, попадает ли точка с произвольными координатами INCLUDEPICTURE "http://www.physchem.chimfak.rsu.ru/Source/NumMethods/Integr-new1.files/image266.gif" \* MERGEFORMATINET в область «внутри» поверхности молекулы или нет.