Быстропеременные функции. Метод Филона

Рассмотренные выше методы интегрирования хороши, если подынтегральная функция и её производные непрерывны и ограничены на отрезке интегрирования.

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

(6.23)

где — высокочастотное колебание (т.е. ω велико), а f (x) — амплитуда. Подынтегральные функции в (6.23) являются быстропеременными (быстро осциллирующими), а их m -ая производная есть величина порядка O (ω ^m).

Функция имеет на отрезке [ a, b ] примерно ω (b – a)/π корней. Так как число корней велико, то для приближенного вычисления интеграла с помощью квадратурных формул придется использовать интерполяционные многочлены высокого порядка или очень большое число узлов интерполирования, что приводит к громоздким вычислениям и, как следствие, большим ошибкам округления.

Будем считать весовой функцией. В качестве узлов интерполирования примем

x_i = (b + a)/2 + (b – a) d_i /2, i = 1, …, n, (6.24)

где числа d_i принадлежат отрезку [–1; 1]. Заменим функцию в интеграле (6.23) на интерполяционный многочлен с узлами x_i:

. (6.25)

Здесь

. (6.26)

Мы получили квадратурную формулу

. (6.27)

Рассмотрим частные случаи. При n = 2, d ₁ = –1, d ₂ = 2 из (6.26), (6.27) получим [1]:

(6.28)

(6.29)

. (6.30)

Проведем анализ влияния погрешности машинных округлений на значения D ₁(p), D ₂(p), вычисляемых по формулам (6.28), (6.29). Пусть t — число двоичных разрядов записи числа. Тогда можно считать, что погрешность вычисления значений sin p, cos p есть величина порядка O (2^{– t}), а погрешность вычисления D ₁(p), D ₂(p) имеет порядок O (2^{– t} / p). Отсюда следует, что при малых значениях p погрешности D ₁(p), D ₂(p) могут быть большими. С другой стороны, используя разложения в ряд Тейлора

,

можно вычислить пределы

Чтобы обойти эту ситуацию при вычислении D ₁(p), D ₂(p) в стандартных программах применяется следующее правило[1]:

(6.31)

где p ₂ — некоторое малое положительное число, определяемое подбором.

Разделим отрезок [ a, b ] на N частей точками

x_i = a + ih, i = 0, 1, …, N; h = (b – a)/ N.

Тогда можем записать обобщенную формулу (6.30) в виде:

(6.32)

Сформулируем алгоритм вычисления интеграла (6.23) по формулам (6.31) — (6.32).