5.5. Метод Симпсона

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

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

Рис. 5.7. Метод Симпсона

Для записи полинома воспользуемся интерполяционной формулой Ньютона (3.6) для трех узлов

где

разделенные разности; расстояние между узлами.

Введем новую переменную тогда и полином (5.27) принимает вид

Теперь вычислим интеграл от полинома

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

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

Интегралы (5.30) подставим в формулы (5.15) и (5.16) и получим уточненное значение интеграла

которое совпадает с формулой Симпсона (5.29).

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

Суммируя разложения около точки для функции в узлах получим, что

тогда интеграл (5.31) принимает вид

Первое слагаемое в правой части формулы (5.32) совпадает с формулой Симпсона, значит, второе слагаемое является главным членом погрешности для интеграла на интервале

Если интеграл вычисляется на интервале путем разбиения его

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

тогда полная погрешность запишется в виде

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

В качестве примера применения метода Симпсона вычислим таблицу значений полного эллиптического интеграла второго рода [36]

Программы 5.3 составлены в соответствии с блок-схемой рис. 5.5.

В основном блоке программы в диалоговом режиме (строка 10) вводятся переменные: количество разбиений интервала интегрирования; граничные значения и шаг изменения аргумента интеграла (5.35). Пределы интегрирования задаются с помощью операций присваивания (строка 20). В теле цикла по аргументу (строки 30-50) размещены операторы обращения к подпрограмме метода Симпсона и вывода на дисплей таблицы результатов.

В подпрограмме метода Симпсона вне цикла вычисляется шаг интегрирования и инициализируется переменная половинным значением подынтегральной функции на левой границе (строки 100-110). Чтобы повторно не вычислять значения подынтегральной функции во внутренних узлах на границах интервалов, на каждом из которых применяется формула Симпсона, соответствующие значения функции удваиваются после цикла. Функции в узлах, совпадающих с пределами интегрирования, суммируются без удвоения. В циле по переменной I (строки 120-150) накапливаются значения узловых

значений подынтегральной функции с коэффициентами- формулы Симпсона

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

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

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

В работе [17] приведена программа на языке Бейсик для численного интегрирования методом Симпсона дискретных функций, заданных в равноотстоящих узлах.

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

дает результат который незначительно отличается от табличного значения эллиптического интеграла

(см. скан)

(см. скан)

(см. скан)