ЛАБОРАТОРНАЯ РАБОТА 3

ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ

1. Введение

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

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

Рассмотрим простейшие, но широко используемые в практических вычислениях формулы: прямоугольников (с центральной точкой), трапеций, Симпсона. Способ их получения состоит в следующем. Разобьем отрезок интегрирования [a, b] на N частей точками

Положим так что

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

После введения шагов интегрирования искомый интеграл можно представить в виде

(3.1)

где

2.1.1 Формула прямоугольников

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

(3.2)

Суммируя в соответствии с (3.1) приближенные значения по всем элементарным отрезкам, получаем формулу прямоугольников для вычисления приближения к I:

(3.3)

В частном случае, когда формула прямоугольников принимает вид

(3.3а)

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

2.1.2 Формула трапеций

На элементарном отрезке заменим подынтегральную функцию интерполяционным полиномом первой степени:

Выполняя интегрирование на отрезке, приходим к локальной формуле трапеций:

(3.4)

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

Суммируя (3.4) по всем отрезкам, получаем формулу трапеций для вычисления приближения к I:

(3.5)

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

(3.5а)

О точности приближения к см. п. 3.3.

2.1.3 Формула Симпсона

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

Напомним, что мы обозначили: а значение в полуцелой точке

Вычисляя интеграл от полинома на отрезке приходим к локальной формуле Симпсона:

(3.6)

Суммируя (3.6) по всем отрезкам, получаем формулу Симпсона для вычисления приближения к I:

(3.7)

Для постоянного шага интегрирования
формула Симпсона принимает вид

(3.8)

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

(3.8а)

К этой записи приходим, если под локальной формулой понимать результат интегрирования по паре элементарных отрезков:

где — интерполяционный полином второй степени для на построенный по значениям в точках Суммируя локальные приближения по всем парам, получим (3.8а). Разумеется, число пар на [a, b] в этом случае должно быть целым, т. е. N — четным.

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

2.2. Погрешность квадратурных формул

Один из возможных способов оценки точности построенных формул состоит в следующем. Рассмотрим интеграл по элементарному отрезку:

Выберем на этом отрезке какую-либо «опорную» точку x = z и разложим подынтегральную функцию в ряд по формуле Тейлора относительно этой точки:

— остаточный член используемой формулы Тейлора.

Вычисляя интеграл от последней суммы, получаем представление в виде:

(3.9)

где коэффициенты A, B, … зависят от значения производных в точке z: …

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

(3.10)

со своими коэффициентами r, s, q.

Заменяя в (3.10) каждое из значений функции f по формуле Тейлора относительно той же точки z, получим представление приближенного значения в виде, аналогичном (3.9):

(3.11)

Сравнивая представления (3.9) и (3.11), обнаруживаем, что кроме первых слагаемых в (3.9), (3.11) совпадает еще некоторое количество (p – 1) слагаемых, так что … Несовпадающие слагаемые характеризуют ошибку квадратурной формулы. Оценивая величины этих слагаемых, приходим к оценке для локальной (на интервале ) погрешности

где D — числовая константа, а — p-я производная функции

Суммируя локальные погрешности по всем интервалам, получим требуемую оценку погрешности рассматриваемой формулы интегрирования:

(3.12)

где по всему отрезку [a, b] (если шаг интегрирования не постоянен, т. е. то ).

Степень p в (3.12) принято называть порядком точности квадратурной формулы.

Для рассмотренных квадратурных формул полученные таким образом оценки погрешности имеют вид:

— для формулы прямоугольников;

— для формулы трапеций;

— для формулы Симпсона в случае, когда используются узлы с дробным индексом (3.8) и

— для (3.8а).

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

Замечание 2. Полученные оценки погрешности, как следует из их вывода, зависят от гладкости подынтегральной функции. Например, при наличии 4-х (и выше) производных y формула Симпсона обеспечивает 4-й порядок точности. Если же только трижды непрерывно дифференцируема на [a, b], то точность формулы Симпсона на порядок уменьшается.

Если известны оценки для абсолютных величин соответствующих производных, то, используя (3.12), можно a priori (до проведения расчета) определить шаг интегрирования h = const, при котором погрешность вычисленного результата гарантировано не превышает допустимого уровня погрешности e. Для этого, как следует из (3.12), достаточно решить относительно h неравенство

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

2.2.1 Контроль за точностью значения интеграла

Пусть шаг интегрирования h = const, — вычисленное с шагом h приближение к I.

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

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

2.2.2 Экстраполяция Ричардсона

Пусть, как и в предыдущем пункте, h = const, — вычисленное с шагом h приближение к I. Пусть использован метод порядка p, тогда можно оценить значение интеграла по элементарному отрезку:

Измельчив шаг вдвое, получаем

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

а для приближения интеграла I с порядком имеем:

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

2.2.3 Счет с автоматическим выбором шагов интегрирования

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

так что

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

2.3. Приемы вычисления несобственных интегралов

Рассмотрим сходящиеся интегралы следующих двух типов:

причем при (первый тип);

(второй тип).

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

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

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

а) Иногда подходящая замена переменной интегрирования позволяет вообще избавиться от особенности.

В рассматриваемом примере после замены получаем

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

б) Та же цель (избавление от особенности) достигается иногда предварительным интегрированием по частям:

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

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

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

Первый интеграл с требуемой точностью вычисляем аналитически, используя представление подынтегральной функции в окрестности особой точки (x = 0) в виде отрезка ряда по степеням x, который получим после замены соответствующим рядом Тейлора:

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

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

(3.13)

( отводится в качестве допустимого уровня погрешности при вычислении ).

Таким образом, один из параметров (m или d) можно задавать по своему усмотрению, второй – определяется из неравенства (3.13). При этом нужно принять в расчет следующее соображение.

Если то существенно ухудшается оценка погрешности для любой квадратурной формулы, которую мы предполагаем использовать для вычисления так как в качестве коэффициента при (где p — порядок точности выбранной формулы) фигурирует максимальное на отрезке значение p-й производной от подынтегральной функции, которое при в рассматриваемом случае имеет порядок

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

Следовательно, целесообразно задать «не слишком малое» d (например, ), а затем m найти из условия (3.13).

Замечание. Разумеется, если поиск последовательных членов разложения подынтегральной функции затруднителен, то приходится ограничиваться доступными членами. В этом случае из условия типа (3.13) находится параметр d.

Рассмотрим пример вычисления интеграла второго типа:

Можно, как уже отмечалось, свести его к интегралу первого типа. Но мы воспользуемся универсальным приемом выделения особенности. Особенность состоит в том, что верхний предел интегрирования — бесконечность. Представим интеграл в виде суммы двух интегралов: где — интеграл по конечному отрезку [a, A]; — интеграл по Вычисление при заданном A затруднений не вызывает.

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

и требуя, чтобы выполнялось условие

найдем

2.4. Библиографическая справка

О методах численного интегрирования см. [1–3, 5, 8, 9, 10]. Отметим, что большей точности метода при небольшом количестве точек (узлов сетки) позволяют достигать квадратурные формулы Гаусса. Про них можно прочитать в [5], теоретические основы методов Гаусса разбираются в [9, 10].

3. Задание

1. Вычислить с постоянным шагом h = 0,1; 0,02 приближенное значение интеграла

по формулам:

а) прямоугольников;

б) трапеций;

в) Симпсона.

2. Сравнить фактическую погрешность с теоретической (точное значение интеграла I = 36.)

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

2. Выполнить задания п. 1.1–1.3 для

(Точное значение: ).

3. То же задание для интеграла

Объяснить результат, полученный при счете с автоматическим контролем точности (режимы с уменьшающимся шагом и с автоматическим выбором шага), когда начальное значение h = 0,1.

Каким должно быть начальное значение h в этом случае? (Точное значение интеграла

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

с шагом h = 0,1; 0,05; 0,02 и с шагом h = 1/4; 1/8; 1/16, 1/32.

5. Вычислить интеграл

Указание. См. п. 3.4.

4. Контрольные вопросы

1. Получить оценки для погрешности квадратурной формулы трапеций.

2. Выполнить то же задание для формулы прямоугольников (с центральной точкой).

3. То же для формулы Симпсона.

2. Описать алгоритм автоматического выбора шага, основанный на экстраполяции по Ричардсону.