Г Л А В А   5

 

численное  Дифференцирование

и  интегрирование


 

 

§ 1. Численное дифференцирование

 

Аппроксимация производных. Напомним, что производной функции  называется предел отношения приращения функции  к соответствующему приращению аргумента  при стремлении  к нулю:

 

,    .             (5.1)

 

Обычно для вычисления производных используют готовые формулы (таблицу производных) и к выражению (5.1) не прибегают. Однако в численных расчетах на компьютере использование этих формул не всегда удобно и возможно. В частности, функция  может быть задана в виде таблицы значений (полученных, например, в результате численного расчета). В таких случаях производную приближенно можно найти опираясь на формулу (5.1). Полагая  равным некоторому конечному числу, получают приближенное равенство для вычисления производной

 

.                                         (5.2)

 

Это соотношение называется аппроксимацией (приближением) производной с помощью конечных разностей (значения  и  в формуле (5.2) конечны в отличие от их бесконечно малых значений в (5.1)).

Соотношение (5.2) может быть использовано для приближенного вычисления производной от функции, заданной как аналитическим выражением, так и таблично. В первом случае выбор величины  произволен и определяется характером поведения функции. Для получения хорошей точности величину  выбирают достаточно малой, такой чтобы на интервале  функция  была бы монотонна и менялась не существенно. В случае, когда функция задана таблично, величина  равна разности между соседними узлами таблицы в окрестности которых вычисляется производная. При этом, если количество узлов невелико и узлы расположены на большом расстоянии друг от друга, формула (5.2) может давать существенную погрешность. Вопросы оценки погрешностей, возникающих при численном дифференцировании, будут рассмотрены ниже.

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

Использование интерполяционных полиномов. Пусть в точках  известны значения функции : . По табличным данным аппроксимируем функцию  интерполяционным полином степени n:

 

.

 

Тогда для k-той производной от функции  на отрезке интерполирования  получим приближенную формулу

 

.                                      (5.3)

 

Однако на практике редко прибегают к аппроксимации функции одним интерполяционным полиномом, т.е. к глобальной интерполяции, в частности,  из-за свойственной ей большой погрешности. Как правило, пользуются локальной интерполяцией. При этом в окрестности точки, в которой нужно вычислить производную, функцию интерполируют полиномом невысокой степени (например, ).

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

а) Приблизим в рассматриваемой окрестности функцию  интерполяционным полиномом первой степени , т.е. прямой, проходящей через точки  и :

 

Тогда

.                         (5.4)

 

С другой стороны в рассматриваемой окрестности функцию  можно приблизить и так

 

 

В этом случае

.                        (5.5)

 

Мы получили простейшие приближенные формулы для первой производной от функции, заданной таблично. Формулу (5.4) называют левым разностным отношением, а формулу (5.5) – правым разностным отношением. Смысл этих названий нетрудно понять из рисунка 5.1. Заметим, что эти соотношения можно было написать сразу, опираясь на формулу (5.2), полагая, например,  и , и не привлекая интерполяцию в качестве промежуточного звена. Понятно также,  что для получения приближенных формул для второй и высших производных линейного приближения функции недостаточно.

б) Приблизим в рассматриваемой окрестности функцию  интерполяционным полиномом второй степени , т.е. параболой, проходящей через значения функции в точках :

 

.

 

Эта форма записи интерполяционного полинома несколько отличается о ньютоновской, которая была рассмотрена ранее в главе 3 при описании кусочно-квадратичной интерполяции (формула (3.12)). Такая форма записи является более компактной.

Дифференцируя это выражение один раз, получим новую приближенную формулу для первой производной:

 

.     (5.6)

Здесь, в отличие от (5.4) и (5.5) , приближение зависит от x. В частности, для  имеем

.                                    (5.7)

 

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

Дифференцируя полином   два раза, получаем приближенную формулу для второй производной:

 

.                          (5.8)

 

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

Так, например, в случае интерполяции функции  полиномом четвертой степени  можно получить следующие формулы (центральные разностные соотношения) для первой и второй производной:

 

,

 

.

 

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

 

,

 

где  – ошибка интерполяции. Следовательно, погрешность вычисления производной функции  по формуле (5.3) будет определяться производной от ошибки интерполяции .

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

Прежде всего заметим, что все рассмотренные нами формулы для приближенного вычисления производной в конкретной точке имеют следующую структуру:

 

,                                  (5.9)

 

где  – постоянные коэффициенты; суммирование производится по некоторому диапазону табличных данных.

Анализ погрешности формулы (5.9) сводится к следующему. Заменяем все входящие в правую часть (5.9) значения  разложением  по формулам Тейлора относительно точки, для которой рассматривается приближение (5.9). После проведения несложных арифметических преобразований в качестве главного члена правой части получаем приближенное значение производной. Остальные члены будут характеризовать погрешность.

В качестве примера рассмотрим оценку погрешности формулы  (5.5) для точки . Заменим значения  и  по формулам Тейлора относительно точки x:

 

,

 

 .

 

Подставляя эти выражения в правую часть (5.5) получим

 

 

 .

Таким образом, погрешность (ошибка) формулы (5.5) для произвольной точки  будет определяться следующим выражением

 

 

Для точки  получаем

 

 ,  т.е. ,

 

следовательно, для  формула (5.5) является формулой первого порядка точности[1]. Можно показать, что для любой точки  формула (5.5) также будет иметь первый порядок точности [12]. Исключение составляет точка  – середина интервала , для которой

,  т.е.

 

и, следовательно, в этом случае формула (5.5) является формулой второго порядка точности.

Мы рассмотрели лишь один из источников погрешности численного дифференцирования – погрешность аппроксимации, которая определяется величиной остаточного члена. Как уже отмечалось выше, анализ остаточного члена нетривиален. Отметим, лишь, что погрешность аппроксимации при уменьшении шага h, как правило, уменьшается. Рассмотренный пример, позволяет также сделать вывод о том, что погрешность формул численного дифференцирования зависит от того, в какой точке вычисляется производная.

Другой источник погрешности численного дифференцирования связан с погрешностями вычисления значений функции   в узлах и с погрешностями округлений при проведении расчетов на компьютере. Обусловленные этими причинами погрешности, в отличие от погрешности аппроксимации, возрастают с уменьшением шага h. Действительно, если при вычислении значения функции  абсолютная погрешность равна ,  то при вычислении дробей в (5.4) и (5.5) она составит . Поэтому суммарная погрешность численного дифференцирования может убывать при уменьшении шага лишь до некоторого предельного значения, после чего дальнейшее уменьшение шага не повысит точности результата.

 

 

 

§ 2. Численное интегрирование

 

Классификация методов. Ставится задача вычислить определенный интеграл вида

 

,                                       (5.10)

 

где  – непрерывная функция на интервале .

Во многих случаях, когда подынтегральная функция задана в аналитическом виде, определенный интеграл удается вычислить непосредственно с помощью неопределенного интеграла по формуле Ньютона-Лейбница. Она состоит в том, что определенный интеграл равен приращению первообразной  на отрезке интегрирования:

 

.

 

Однако на практике этой формулой часто нельзя воспользоваться по двум основным причинам:

1) вид функции  не позволяет аналитически выразить первообразную  через элементарные функции;

2) подынтегральная функция  задана таблично.

В этих случаях используют приближенные, методы интегрирования.

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

 

,                                (5.11)

 

где R – погрешность вычисления интеграла.

С целью уменьшения погрешности, связанной с аппроксимацией подынтегральной функции, интервал интегрирования  разбивают на n отрезков и на каждом частичном отрезке заменяют подынтегральную функцию аппроксимирующей функцией . Тогда приближенное значение интеграла определяется суммой  частичных интегралов от функций , взятых в пределах от  до :

 

.                                (5.12)

 

Используемые на практике методы численного интегрирования можно сгруппировать в зависимости от способа аппроксимации подынтегральной функции. Ниже приведем краткую характеристику классов наиболее распространенных методов.

Методы Ньютона-Котеса основаны на полиномиальной аппроксимации (интерполяции) подынтегральной функции. Методы этого класса отличаются друг от друга  степенью используемого полинома, от которого зависит количество узлов, где необходимо вычислять функцию . В методах Ньютона-Котеса отрезок интегрирования разбивается, как правило, на отрезки равной длины, величина которых определяется как  и называется шагом интегрирования. Алгоритмы этих методов просты и легко поддаются программной реализации.

Сплайновые методы базируются на аппроксимации подынтегральной функции сплайнами. Методы различаются по типу выбранных сплайнов. Такие методы имеет смысл использовать в задачах, где алгоритмы сплайновой аппроксимации применяются для обработки данных.

В методах наивысшей алгебраической точности  используют неравноотстоящие узлы, расположенные так, чтобы обеспечить  минимальную погрешность интегрирования при заданном количестве узлов. Методы различаются способом выбора узлов. Наиболее широкое применение получили методы Гаусса, в которых узлы интегрирования выбираются как корни полиномов Лежандра.

В методах Монте-Карло узлы выбираются с помощью генератора случайных чисел, результат в итоге носит вероятностный характер. Методы оказываются особенно эффективны при вычислении кратных интегралов.

Независимо от выбранного метода в процессе численного интегрирования необходимо вычислить приближенное значение интеграла и оценить погрешность R.

Погрешность, возникающая при численном интегрировании, также как и при численном дифференцировании, имеет два основных источника. Первым источником погрешности является приближенная замена подынтегральной функции аппроксимирующей функцией – погрешность аппроксимации. Как будет показано ниже, погрешность аппроксимации уменьшается с увеличением количества отрезков разбиения интервала интегрирования n за счет более точной аппроксимации подынтегральной функции. Другой источник погрешности – неточности в вычислении подынтегральной функции в узловых точках и ошибки округления. Эта погрешность возрастает с ростом n и с некоторого значения  начинает преобладать над погрешностью аппроксимации. Это обстоятельство должно предостеречь от выбора чрезмерно большого числа n.

Интерполяционные методы Ньютона-Котеса. Способ получения формул для вычисления приближенного значения интеграла в методах Ньютона-Котеса состоит в следующем. Разобьем интервал интегрирования  на n элементарных отрезков. Точки разбиения  будем называть узлами интегрирования, а расстояния между узлами  шагами  интегрирования. В частном случае шаг интегрирования может быть постоянным . На каждом частичном отрезке интегрирования  будем аппроксимировать подынтегральную функцию интерполяционным полиномом. В этом случае  вычисление частичных интегралов в формуле (5.12) не составит труда.

Методы прямоугольников. Рассмотрим сначала простейшие методы из класса методов Ньютона-Котеса, в которых подынтегральную функцию  на отрезке интегрирования  заменяют полиномом нулевой степени, т.е. константой: . Подставляя  в формулу (5.12) и выполняя интегрирование, получаем

 

.                 (5.13)

 

Таким образом, приближенное значение интеграла определяется суммой площадей прямоугольников, одна из сторон которых есть длина отрезка интегрирования , а другая – аппроксимирующая константа. Отсюда происходит и название методов.

Ø Замечание. К формуле (5.13) можно также придти, если заменить определенный интеграл интегральной суммой, непосредственно опираясь на определение понятия определенного интеграла.<


Заметим, что замена подынтегральной функции константой неоднозначна, так как константу  можно выбрать равной значению подынтегральной функции в любой точке интервала интегрирования.

Выбирая константу  равную значению подынтегральной функции в левой (рис. 5.2 а) или правой (рис. 5.2 б) границах отрезка  приходим к формулам левых и правых  прямоугольников, соответственно:

 

,                (5.14)

 

.                   (5.15)

 

В формулах (5.14) и (5.15) введены обозначения , , которыми мы будем пользоваться и в дальнейшем.

В случае постоянного шага интегрирования формулы (5.14) и (5.15) приобретают вид

 

,             (5.16)

 

.              (5.17)

 

Наиболее широко на практике используется формула средних прямоугольников, в которой значение константы  (высота прямоугольника) выбирается равной значению подынтегральной функции в средней точке отрезка интегрирования (рис. 5.3):

 

,     .                (5.16)

 

В случае  постоянного шага интегрирования, когда ,  формула средних прямоугольников (5.16) будет иметь следующий вид

 

,     (5.17)

 

где    .

Из трех рассмотренных методов прямоугольников метод средних прямоугольников является наиболее точным.

Ø Замечание. Если подынтегральная функция  задана не аналитическим выражением, а таблично, то формула средних прямоугольников оказывается неприменима (без привлечения дополнительной интерполяции), так как значения функции известны лишь в узловых точках. В этом случае пользуются либо формулами левых или правых прямоугольников, либо используют другие методы.<

Метод трапеций. На частичном отрезке интегрирования  заменим подынтегральную функцию  интерполяционным полиномом первой степени, т.е. прямой, проходящей через точки  и  (рис. 5.4):

 

.

 

Поставляя это выражение в формулу (5.12) и выполняя интегрирование по частичным отрезкам, приходим к формуле трапеций:

 

.

(5.18)

 

Название метода связано с тем, что интеграл по отрезку  заменяется площадью трапеции с основаниями, равными значениям функции  на краях отрезка, и высотой, равной  (рис. 5.4).

В случае постоянного шага интегрирования формула трапеций принимает вид

 

.  (5.19)

 

Для удобства вычислений и оптимизации вычислительных погрешностей формулу (5.19) часто записывают следующим образом:

 

.                         (5.20)

 

Ø Замечание. Несмотря на то, что в методе трапеций используется более высокий порядок интерполяции подынтегральной функции (кусочно-линейная интерполяция), по сравнению с методом прямоугольников, который использует интерполяцию нулевого порядка (кусочно-постоянная), точность метода трапеций оказывается ниже точности метода средних прямоугольников. Более высокая точность метода средних прямоугольников объясняется “удачным” выбором узловых точек, в которых вычисляется высота элементарного прямоугольника.<

Метод Симпсона (метод парабол). Разобьем интервал интегрирования  на четное число n равных отрезков с шагами . На каждом отрезке , содержащем три узла, заменим подынтегральную функцию интерполяционным полиномом второй степени, т.е. параболой

 

 

Интегрируя это выражение на отрезке , получим

 

 

.                                                   (5.21)

 

Приближенное значение интеграла на интервале  получим суммированием частичных интегралов (5.21) по всем отрезкам :

 

 

.     (5.22)

 

Полученное соотношение называется формулой Симпсона или формулой парабол.

Если подынтегральную функцию  интерполировать полиномом второй степени на отрезке  с привлечением дополнительной точки  – середины отрезка, то в этом случае число отрезков разбиения n может быть произвольным (не обязательно четным), а формула Симпсона в этом случае будет иметь вид

 

 

.     (5.23)

 

Напомним, что , .

Ø Замечание 1. Формула (5.22) пригодна для вычисления интегралов от функций, заданных как аналитическим выражением, так и таблично, тогда как формула (5.23) применима только в тех случаях, когда подынтегральная функция задана аналитически.<

Ø Замечание 2. Формулы, используемые для приближенного вычисления интеграла, называются квадратурными формулами. Нетрудно заметить, что все полученные нами формулы имеют следующую структуру:

,

 

где  – значения подынтегральной функции в узловых точках;  весовые коэффициенты, независящие от функции .<

 

Погрешность формул Ньютона-Котеса. Рассмотрим подробно один и возможных способов оценки погрешности формул численного интегрирования на примере метода средних прямоугольников. Для этого запишем выражение для интеграла на отрезке , полученное методом средних прямоугольников

 

,                                (5.24)

откуда

 

.                                (5.25)

 

Для оценки погрешности интегрирования  разложим подынтегральную функцию  в ряд Тейлора около средней точки

 

.       (5.26)

 

В малой окрестности точки x в разложении (5.26) можно ограничится небольшим количеством членов ряда. Поэтому, подставляя в (5.25) вместо функции  ее тейлоровское разложение (5.26) и интегрируя его почленно, можно вычислить интеграл с любой заданной точностью

 

 

                            (5.27)

 

При интегрировании и подстановке пределов получаем, что все интегралы от членов ряда (5.26), содержащих нечетные степени , обращаются в нуль.

Подставляя полученное соотношение в формулу (5.25), получим

 

.

 

При малой величине шага интегрирования h основной вклад в погрешность  будет вносить первое слагаемое, которое называется главным членом погрешности  вычисления интеграла на отрезке

 

.                                    (5.28)

 

Главный член полной погрешности для интеграла на всем интервале  определиться путем суммирования погрешностей на каждом отрезке :

.                             (5.29)

 

На основании (5.29) можно получить оценку для абсолютного значения погрешности:

 

 ,                (5.30)

 

где

.

 

Формула (5.30) представляет собой теоретическую оценку погрешности вычисления интеграла методом средних прямоугольников, эта оценка является априорной, так как не требует знания значения вычисляемого интеграла. Степень шага h, которой пропорциональна величина, называется порядком метода интегрирования. Метод средних прямоугольников имеет второй порядок точности.

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

 

§        методы левых и правый прямоугольников – ;

§        методы средних прямоугольников –  ;

§        метод трапеций – ;

 

§        метод Симпсона (5.22) – ;

§        метод Симпсона (5.23) – .

 

Методы левых и правых прямоугольников являются методами первого порядка. Методы средних прямоугольников и трапеций имеют второй порядок точности, при этом метод трапеций обладает вдвое большей по абсолютной величине погрешностью по сравнению с методом средних прямоугольников. Поэтому, если подынтегральная функция задана аналитически, то предпочтительнее из методов второго порядка применять метод средних прямоугольников вследствие его меньшей погрешности. Метод Симпсона имеет четвертый порядок точности с очень малым численным коэффициентом. Формула Симпсона позволяет получить очень высокую точность, если четвертая производная подынтегральной функции не слишком велика. В противном случае, методы второго порядка точности могут дать большую точность, чем метод Симпсона. Например, для функции  формула трапеций при  для интеграла в пределах  дает точный результат, равный 4, тогда как по формуле Симпсона получим результат, несовпадающий даже по знаку ().

Вычисление интеграла с заданной точностью. Используя приведенные выше оценки можно априори (до проведения расчета) определить шаг интегрирования h, при котором погрешность вычисленного результата гарантированно не превысит допустимый уровень погрешности e. Однако на практике пользоваться априорными оценками погрешности не всегда удобно. Тогда контроль за точностью получаемого результата можно организовать следующим образом. Пусть вычисления проводились с постоянным шагом h, – вычисленное с шагом h приближенное значение интеграла I.  Если затем вычислить приближенное значение  с шагом , то в качестве оценки погрешности последнего вычисленного значения можно рассматривать величину

 

.

 

При необходимости вычислить результат с заданной точностью e  вычисления повторяют с последовательно уменьшающимся (вдвое) шагом до тех пор, пока не выполнится условие

 

.                                 (5.30)

 

Можно применить указанное правило для контроля за локальной погрешностью на каждом частичном интервале. При этом длина очередного интервала , посредством последовательного уменьшения (или увеличения) начальной длины вдвое, устанавливается такой, чтобы выполнялось неравенство

 

.                          (5.31)

 

Тогда в худшем случае ошибка вычисленного значения интеграла на всем интервале интегрирования не будет превосходить сумму локальных погрешностей , т.е. не будет превосходить заданного уровня погрешности.

Способ вычисления интеграла с автоматическим выбором шага имеет то преимущество, что он “приспосабливается” к особенностям подынтегральной функции: в областях резкого изменения функции шаг уменьшается, а там, где функция меняется слабо, – увеличивается. Такого рода алгоритмы называются адаптивными, т.е. приспосабливающимися. Использование адаптивных алгоритмов позволяет сократить затраты машинного времени  без потери точности.

Особые случаи численного интегрирования.

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

 

.

Для вычисления каждого из интегралов в правой части можно использовать любой из рассмотренных выше методов.

б) Несобственные интегралы. Напомним, что к такому типу интегралов относятся интегралы, которые имеют хотя бы  одну бесконечную границу интегрирования или подынтегральную функцию, обращающуюся бесконечность хотя бы в одной точке интервала интегрирования.

Рассмотрим сначала интеграл с бесконечной границей интегрирования, например, интеграл вида

 

.

 

Один из универсальных способов вычисления такого интеграла состоит в следующем. Представим интеграл в виде суммы двух интегралов:

 

,

 

где A – некоторое большое число. Вычисление первого интеграла на конечном отрезке  не вызывает вопросов. Выбор числа A производят таким образом, чтобы в пределах допустимой погрешности вторым интегралом можно было пренебречь.

В случае, когда подынтегральная функция обращается в бесконечность в некоторой точке  (), то интеграл приближенно можно представить в виде суммы двух интегралов

 

,

 

где d – малая величина. На отрезках  и  подынтегральная функция не содержит особенностей и указанные интегралы могут быть вычислены любым из рассмотренных выше методов.

Ø Замечание. При вычислении несобственных интегралов не стоит забывать о том, что иногда подходящая замена переменной интегрирования позволяет вообще избавится от особенности. Так, например, при вычислении интеграла

 

 

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

 

,

 

который не имеет особенностей и вычисляется с требуемой точностью по любой квадратурной формуле.<

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

 

,    .                       (5.31)

Предполагается, что     не сильно меняющаяся на отрезке  функция, и частота колебаний подынтегральных функций в (5.31) определяется величиной w. Будем считать, что . Если воспользоваться комплексным представлением результатов  (i – мнимая единица), то две формулы (5.31) можно объединить в одну:

 

.                                   (5.32)

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

Рассмотрим один из возможных подходов к построению квадратурных формул, учитывающих специфику подынтегральной функции.

Как обычно отрезок интегрирования разбивается на n интервалов. При этом узлы интегрирования определяются как , , , . В частном случае шаг интегрирования может быть постоянным . Искомый интеграл представим в виде суммы интегралов по каждому частичному отрезку:

 

,                  (5.33)

где

.                                 (5.34)

Для вычисления интегралов  будем использовать аппроксимацию не всей подынтегральной функции, а только неосциллирующей ее части, т.е. .

Заменим  на отрезке  полиномом нулевой степени, т.е. константой. Причем значение константы выберем равным значению функции  в середине отрезка : .

Тогда

.

Интеграл от функции  может быть вычислен аналитически:

 

 ,   (5.35)

 

где .

Следовательно, сумма (5.33) может быть записана в виде

 

.         (5.36)

 

Для случая постоянного шага интегрирования

 

,   .       (5.37)

 

Заметим, что при  формулы (5.36) и (5.37) естественным образом переходят в формулы средних прямоугольников (5.16) и (5.17), полученные нами ранее. По сути, отличие данных результатов от результатов, полученных ранее состоит в единственном множителе . Оказывается, учет этого множителя позволяет существенно снизить требования к величине шага интегрирования.

Если неосциллирующую часть подынтегральной функции аппроксимировать полином первой или второй степени, то можно получить аналоги формул трапеций и Симпсона для вычисления интегралов от быстроосциллирующих функций. Рассмотренный нами подход к вычислению интегралов вида (5.32) впервые был предложен в работах Филона[2].

Вычисление кратных интегралов. Для вычисления кратных интегралов вида

 

 

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

 

.

 

Если областью интегрирования G является прямоугольник: , , то для вычисления двойного интеграла можно использовать любую из рассмотренных выше формул. Например, применение формулы средних прямоугольников дает следующий результат

 

.          (5.32)

 

Если область интегрирования G непрямоугольная, например,

 

,

 

то её можно привести к прямоугольному виду с помощью замены

 

,  .

 

С повышением кратности интегралов резко возрастает объем вычислений и рассмотренный подход становится неэффективным. Например, если мы разбиваем интервал изменения каждой переменной всего на десять частей, то для вычисления тройного интеграла нам потребуется вычислить сумму тысячи слагаемых, а при вычислении 10-кратного интеграла, потребуется сумма, количество слагаемых в которой определяется числом . Вычисление такой суммы затруднительно даже на самых быстродействующих современных компьютерах. В этом случае применяют другие методы численного интегрирования, среди которых особое место занимает метод статистических испытаний.

Методы Монте-Карло. Методы статистических испытаний, называемые также методами Монте-Карло[3], применяются к решению разнообразных задач вычислительной математики и, в том числе, для вычисления интегралов. Рассмотрим вначале один из простых вариантов метода Монте-Карло, который можно интерпретировать как статистический метод прямоугольников (рис. 5.5). На отрезке интегрирования  выберем N случайных точек , являющихся значениями случайной величины x с равномерным распределением на данном отрезке. Для каждой точки вычислим площадь прямоугольника, одна сторона которого равна , а вторая равна значению функции в данной точке : . Вследствие случайности узла , значение площадей  также будет носить случайный характер. В качестве приближенного значения интеграла можно принять результат усреднения площадей :

 

.         (5.33)

 

Погрешность вычисления интеграла будет уменьшаться с ростом числа испытаний N по закону  [8].

Полученная формула формально совпадает с формулой правых прямоугольников (5.17), но отличие состоит в том, что в формуле (5.17) узлы интегрирования расположены регулярно, а в данном случае расположение узлов носит случайный характер.

Формула (5.33) непосредственно обобщается на кратные интегралы

 

,

 

где  – объем области интегрирования. Например, для двукратного интеграла с прямоугольной областью интегрирования имеем

 

.              (5.34)

 

Теоретическое обоснование рассмотренного варианта метода Монте-Карло для вычисления интегралов состоит в следующем. Пусть x – случайная величина с равномерным распределением на отрезке . Это означает, что ее плотность вероятности задается соотношением

 

 

Тогда любая функция  также будет случайной величиной, и ее математическое ожидание равно

 

.

 

Или, читая это равенство наоборот, получим

 

.

 

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

 

,

 

что в итоге приводит к соотношению

 

,

 

которое совпадает с формулой (5.33).

В другом варианте метода Монте-Карло исходный интеграл приводится к виду

 

,

 

где  на интервале . Генерируются пары значений  и  двух случайный величин x и y, которые можно рассматривать как координаты точек в единичном квадрате (рис. 5.6). При равномерном распределении точек в квадрате за приближенное значение интеграла принимается отношение количества точек S, попавших под кривую , к общему числу испытаний (точек) N

 

.

 

Этот алгоритм также обобщается на кратные интегралы.

Для использования методов Монте-Карло при вычислении интегралов, как и в других их приложениях, необходимо вырабатывать последовательности случайных чисел с заданным законом распределения вероятностей. Наиболее распространенный способ выработки случайных чисел на компьютере состоит в использование специального алгоритма получения таких чисел. Поскольку эти числа генерируются по заданному алгоритму, то получаемую последовательность (одну и ту же) можно воспроизвести сколько угодно раз. В этом смысле получаемые числа не совсем случайны (псевдослучайны). Тем не менее, в рамках данной последовательности, получающиеся числа обладают всеми необходимыми статистическими свойствами, характерными для значений истинно случайной величины.

 

 

Вопросы для самоконтроля

1.     В чем особенность численного дифференцирования с использованием интерполяционных полиномов?

2.     Каковы основные источники погрешности численного дифференцирования?

3.     В чем состоит идея методов Ньютона-Котеса для приближенного вычисления определенных интегралов?

4.     Как влияет на точность численного интегрирования величина шага h?

5.     Можно ли добиться неограниченного уменьшения погрешности интегрирования путем последовательного уменьшения шага?

6.     Почему метод средних прямоугольников неприменим для численного интегрирования таблично заданных функций?

7.     Каковы преимущества формулы Симпсона  по сравнению с формулой трапеций и следствием чего являются эти преимущества?

8.     В каком случае формула Симпсона дает точное значение интеграла?

9.     Какой подход используется на практике для вычисления интеграла с заданной точностью?

10. Какова идея вычисления определенного интеграла методом Монте-Карло?

 

 

 

 



[1] Порядком точности называется степень h, которой пропорциональна погрешность.

[2] L.N.G. Filon, Proc. Roy. Soc., Edinb., 1928-1929, p. 49.

[3] Название метода произошло от названия города Монте-Карло, знаменитого своими казино, в которых играют в рулетку, являющуюся одним из простейших генераторов случайных чисел.