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 
 в область «внутри» поверхности молекулы или нет.