Принцип метода Симпсона состоит в замене подынтегральной функции f(х) интерполяционным многочленом Ньютона второй степени. Тогда для каждого элементарного отрезка [хi,хi+1] имеем следующее значение площади подынтегральной кривой:
.
Для всего отрезка интегрирования [a,b] формулой Симпсона:
Данное выражение называется формулой Симсона. Оно относится к формулам повышенной точности и является точной для многочленов второй и третьей степени.
Рисунок 31 – Геометрическая интерпретация численного интегрирования методом Симпсона
Приведём программу, реализующую вычисление определённого интеграла методом Симпсона с заданнойточностью. В качестве подынтегральной будем использовать функцию:
.
Рассмотренные формулы численного интегрирования требуют чёткого указания количества разбиений отрезка интегрирования. Однако классическое использование численного метода предполагает вычисление значения (корня, интеграла и т.д.) с заданной точностью.
Точность любой формулы численного интегрирования зависит от величины отрезка разбиения D.
Будем вычислять значение интеграла при разных значениях D (D1, D2, D3,…), где Di+1 = 2Di. Как только разница между значением интеграла, вычисленного при Di и интеграла, вычисленного при Di+1, станет меньше, чем значение e, будем считать, что интеграл вычислен с заданной точностью e.
Данный метод интегрирования с заданной точностью прост в реализации, однако он требует значительных избыточных вычислений, что приводит к повышению затрат времени на вычисление.
program simp;
function f(x: real): real;
Begin
f:=1/x
end;
Var
a,b,e: real;
i: integer;
xa,xab,xb,dx,s1,s: real;
n: integer;
Begin
writeln('[a,b],e');
readln(a,b,e);
{вычисление интеграла с количеством разбиений равным 1, т. е. одной фигурой с основанием равным [a,b]}
n:=1;
dx:=(b-a)/n;
s:=dx*(f(a)+4*f(a+dx/2)+f(b))/6;
Repeat
n:=n*2; {удвоение количества отрезков разбиения}
s1:=s;
s:=0;
{вычисление длины отрезка – основания прямоугольника (дельта) при новом количестве разбиений}
dx:=(b-a)/n;
{суммирование площадей - нахождение интеграла при заданном количестве разбиений}
for i:=0 to n-1 do
Begin
xa:=a+dx*(i);
xb:=xa+dx;
xab:=xa+dx/2;
s:=s+dx*(f(xa)+4*f(xab)+f(xb))/6;
end;
until abs(s-s1)<=abs(e);
writeln('int=',s);
End.
В данной программе используется подпрограмма функция f, которая вычисляет подынтегральную функцию .
Переменная s1 используется для сохранения значения интеграла, вычисленного при вдвое меньшем количестве разбиений.