Численные методы вычисленияинтегралов. Метод Ньютона-Котеса. Метод Гаусса
1. Численные методы вычисленияинтегралов. Постановка задачи
Решая физические задачи, часто приходитсявычислять значения определённых интегралов от функций />. Во многих случаях, в видутого, что подлежащий вычислению интеграл не выражается через элементарныефункции, прибегают к приближённым численным методам.
Прежде всего, рассмотрим случай, когда /> - конечный интервал.
В таком случае, как известно, функция /> является ограниченной,т.е. />. В этом случае наиболеечасто применяемый численный метод интегрирования состоит в том, что интеграл от/> заменяется некоторойлинейной комбинацией значений /> в /> точках />:
/>(1)
Формула (1) называется квадратурнойформулой, а коэффициенты /> -квадратурными коэффициентами или весами, абсциссы /> -узлами квадратурной формулы.
Методы численного интегрированияклассифицируются в зависимости от того, заданы ли значения аргумента черезравные промежутки или нет. Так методы Ньютона-Котеса требуют, чтобы значения /> были заданы с постояннымшагом, а методы Гаусса не налагают такого ограничения. Перейдём к рассмотрениюэтих методов.
2. Методы Ньютона-Котеса
Пусть /> различныеточки отрезка />, служащие узламиинтерполяции для некоторой интерполирующей функцию /> функции/>. Тогда имеем:
/>(2)
где /> -остаточный член. Предположим, что
/>(3)
причём /> подобранытак, чтобы все интегралы
/>(4)
можно вычислить точно. Тогда мы получаемквадратурную формулу
/> (5)
2.1 Формула трапеций
/>
/> Частным случаем методов Ньютона-Котесаявляется квадратурная формула трапеции. Подынтегральную функцию будеминтерполировать по формуле Лагранжа, в том случае, когда на каждом отрезкеделения принимается линейная интерполяция, а результаты суммируются (рис 1):
/>
Рис. 1.
а) графический вывод:
Определённый интеграл />, как известно, задаётплощадь /> криволинейной трапеции />, поэтому, вписав ломаную вдугу кривой />, мы получаем, что площадькриволинейной трапеции можно приближённо вычислить как сумму площадей трапеций:
/>(6)
Между тем, очевидно, что
/>(7)
Так как, в методахНьютона-Котеса, />, учитывая (6)получаем:
/> (8)
или, соединяя подобныечлены, имеем:
/> (9)
Формула (9) – называется формулойтрапеций.
б) Аналитический вывод:
Выведем формулу трапециианалитическим способом. Для этого используем интерполяционный многочленЛагранжа для отрезка />, построиммногочлен первой степени, который на концах отрезка принимает заданные значения/>. Ясно, что в таком случаеинтерполирующая функция /> имеетвид:
/>(10)
т.к. в методеНьютона-Котеса />, учитывая (3) и(4), из (10) получаем:
/>
/> (11)
Аналогично, />, т.е.
/> (12)
Таким образом, получаемформулу:
/> (13)
тогда, используя свойствоаддитивности оператора интегрирования, имеем:
/> (14)
где />. Получили формулу (14)трапеций, которая естественно, совпадает с (9).
2.2 Формула Симпсона
Рассмотрим методНьютона-Котеса (т.е. />), в случаеинтерполяции подинтегральной функции квадратичными функциями на каждоминтервале деления. В данном случае мы имеем дело с параболическиминтерполированием, поэтому на каждом интервале />,необходимо знание значения функции /> в трёхточках (т.к. /> имеет 3неизвестных параметра – коэффициенты />). Вкачестве третьей точки на каждом отрезке /> -выбирается середина этого отрезка, т.е. точка />.
Вывод формулы Симпсонабудем производить аналитически. Как и в предыдущем случае применяеминтерполяционный многочлен Лагранжа, для интерполирования функции />, на отрезке />, при чём считаем, что намизвестны значения />. Тогда,очевидно, что многочлен Лагранжа имеет вид квадратичной функции:
/> (15)
Интегрируя (15) наотрезке /> будем иметь формулу:
/>(16)
используя свойствоаддитивности интеграла, получаем:
/>(17)
где />является четным числом(/> — число делений отрезка />, т.е. число равных отрезковразбиения).
Формула (17)-называется формулойСимпсона.
Приняв обозначения />, получаем привычный вид квадратурныхформул:
а) Формула трапеций:
/>(18)
б) Формула парабол(Симпсона) (при />)
/>(19)
2.3 Метод Ромберга
Пусть промежутокинтегрирования разбит на /> равныхчастей и для этого разбиения по формуле трапеции получено значение />. Значение /> - совпадает со значениемвычисляемого интеграла, если интегрируемая функция /> линейна,т.е. является многочленом первой степени. По формуле:
/>(20)
называемой формулойРомберга, построим /> — схему:
/>(21)
Оказывается, что дляинтегрируемых по Риману функций, все столбцы и строки /> — схемы сходятся кисходному значению интеграла.
Пример: Выписать явные формулы для фрагмента /> — схемы:
/>
/>
/>/>
/>
/>
Решение:
Пусть />Тогда
/>
/>
/>
/>
/>
/>
3. Квадратурныеформулы Гаусса
Во всех приведенных досих пор формулах численного интегрирования Ньютона-Котеса и во всех формулах,получаемых методом Ромберга, используются равноотстоящие узлы. В случаеквадратурных формул Гаусса это уже не так. Иначе говоря, смысл квадратурныхформул Гаусса состоит в том, чтобы при наименьшем возможном числе узлов точноинтегрировать многочлены наивысшей возможной степени. Можно показать, что при />гауссовых узлах пополученной формуле можно точно интегрировать многочлены степени />.
/>(22)
Для количества узлов исоответствующих значений />и /> — составлены таблицы,которые позволяют вычислять интегралы по формуле (22).
Для понимания сути этихтаблиц рассмотрим пример.
Пример:
Пусть нам нужно составитьквадратурную формулу с двумя узлами />, покоторой точно интегрируются многочлены до />степеньвключительно.
Решение: Искомая формула имеет вид:
/>,(23)
где /> - остаток, которыйобращается в нуль, для
/>, при />.
Тогда, подставляя в (23)имеем:
/>(24)
Отсюда, приравниваякоэффициенты при /> />, справа и слева, получаемсистему уравнений:
/>(25)
Ее решение имеет вид:
/>(26)
Следовательно, искомаяквадратурная формула такова:
/>.(27)
Ясно, что если нам нужновычислить интеграл со многими узловыми точками, действуем следующим образом:
а) промежутокинтегрирования /> делим на /> — равных промежутков и накаждом маленьком промежутке /> применяемформулу Гаусса с неравноотстоящими узлами (27);
б) полученные результатыскладываем.
В случае, когда />, оказывается, что узловымиточками при делении отрезка на /> — частейявляются корни соответствующих многочленов Лежандра.
Для вычисления кратныхинтегралов, их сводят обычно к повторным интегралам, а далее применяют те жесамые кубатурные формулы для каждого значения узловых точек, что и в одномерномслучае. Однако, надо иметь в виду, что кратные интегралы значительно сложнеевычислять с заданной точностью.
Точность произведённыхвычислений зависит от точности аппроксимации подынтегральной функциимногочленами.
4. Оценка интегралов
При численном интегрировании наряду сприближёнными формулами представляет также интерес нахождение нижних и верхнихграниц интегралов. Рассмотрим два метода оценки интегралов:
а) оценка интеграла вслучае, когда подинтегральная функция />,удовлетворяет условию:
/> для /> (28)
б) общий случай.
Рассмотрим интеграл:
/> (29)
где />, />. Не умоляя общность, будемсчитать, что />, />, тогда (Рис. 1) ясно, что
/>/>/>/>/>/>/>/>/> />
К Е
N
М
0 /> /> /> />
Рис. 1
0 /> />/>/>
Площадь криволинейнойтрапеции /> заключена между площадями aMNb и aKEb, т.е.
/>(30)
Очевидно, что
/>(31)
/>(32)
Таким образом, для оценкиинтеграла в случае />, имеем:
/>(33)
если же />, неравенство (33)заменяется на обратное.
б) Другой принцип грубой, но зато общейоценки значения интеграла, основан на «монотонности» интеграла. При этомспособе подынтегральную функцию приближают снизу и сверху интегрируемыми взамкнутом виде функциями /> и />, т.е.
/>, />(34)
Тогда
/>(35)
5. Вычислениеинтегралов методом Монте-Карло
Пусть нам нужно вычислитьинтеграл:
/>(36)
В случае, когда методы Ньютона-Котеса иГаусса работают плохо, приходится обращаться к вероятностным методам случайногопоиска. К таким методам относится метод Монте-Карло.
Для вычисления интеграла(36) методом Монте-Карло, заменим переменную интегрирования /> таким образом, чтобыпределы интегрирования /> отобразилисьсоответственно в />. Для этого нужновоспользоваться преобразованием:
/>(37)
тогда интеграл (36) принимает вид:
/>(38)
Для вычисления жеинтеграла на /> имеем формулу:
/>(39)
где /> - случайные числа,равномерно распределённые на />. Такимобразом, по методу Монте-Карло, интеграл (36) считается по формуле:
/>(40)
где /> - равномернораспределённые случайные числа из промежутка />.
Аналогично, для кратныхинтегралов. Получаем:
/>(41)
где /> - случайные точки,равномерно распределённые на квадрате /> (Здесьзнак «/>» означает декартовоепроизведение).
В случае, когда областьинтегрирования является сложным множеством /> (рис.6), пользуемся прямоугольником />,который описывается вокруг множества />. Иинтеграл по множеству /> заменяеминтегралом по прямоугольнику />,который уже умеем вычислять по формуле (41). Замена интеграла по множеству /> производится соотношением:
/> (42)
где
/> (43)
таким образом:
/> (44)
который легкорассчитывается по формуле (41).
Аналогично вычисляются итрёхкратные интегралы. Этот подход легко обобщается для n-кратных интегралов.
Литература
1. Р.В. Хемминг.Численные методы, Наука, М.,1998
2. Коллатц.,Ю.Альбрехт. Задачи по прикладной математике. Мир, М.,1998.
3. Т.Шуп. Решениеинженерных задач на ЭВМ. Мир, М., 1992.
4. К.Бреббия, Ж.Теллес, Л. Врубел.Методы граничных элементов. Мир, М.,1987.
5. И.С.Берехин.,Н.П.Жидков. Методы вычислений, ч.1., М.,1982.