Pascal: Вычисление интеграла методом Симпсона

Для того, что бы понять, как это работает нужно обратиться к геометрическому смыслу интеграла.

Интеграл численно равен площади криволинейной трапеции, ограниченной кривой y(x), прямыми x=a, x=b и отрезком [a; b] оси Ox.

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

с геометрической точки зрения означает, что график функции y(x) заменен другой кривой j(x), состоящей из дуг парабол: каждая сдвоенная дуга кривой y(x) заменяется параболой. За приближенное значение интеграла I принимается площадь криволинейной трапеции, ограниченной кривой j(x), прямыми x=x0, x=x2n и отрезком [x0, x2n]. Таким образом, решение сводится к программированию алгоритма вычисления площади криволинейной трапеции.

Программа, которая вычисляет такой интеграл:


program Simpson;
uses crt;
var
x,a,b,h,s:real;
n:integer;

function Y(t:real):real;
begin
Y:=sqrt(1+t*t+t*t*t*t);
end;
begin
clrscr;
write('Отрезок интегрирования [a,b] ?');
read(a,b);
write('На сколько частей нужно разделить отрезок? n=');
read(n);
h:=(b-a)/n;
s:=0; x:=a+h;
while x<b do
begin
s:=s+4*Y(x);
x:=x+h;
s:=s+2*Y(x);
x:=x+h;
end;
s:=h/3*(s+Y(a)-Y(b));
writeln;
writeln('Интеграл = ',s);
readkey;
end.

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

Код на C# здесь

8 Responses

  1. Berlioz 25.06.2012 / 09:50

    Слишком маленькими отрезки делать тоже не стоит. Я сделал отрезкок 1/1000000, вычисления дали большую погрешность; а при отрезке 1/100000 все было приемлемо.

  2. саша 18.09.2012 / 22:07

    все понял кроме этой части
    begin
    s:=s+4*Y(x);
    x:=x+h;
    s:=s+2*Y(x);
    x:=x+h;
    end;
    s:=h/3*(s+Y(a)-Y(b));
    writeln;
    можешь объяснить?

  3. Юрий 24.09.2013 / 00:26

    Спасибо за код.

    Кстате в цикле можно еще так организовать чередование индексов 2 и 4:
    взять любую переменную равную единице, например c:=1;
    и в цикле считать сумму примерно так:
    s:=s+(3+c)*f(x);
    c:= -c;
    тогда каждый раз в скобках (3+с) будет 2 либо 4.

  4. Дмитрий 26.12.2013 / 14:44

    Спасибо, все прекрасно работает, но не могу понять этой формулы в коде:
    s:=h/3*(s+Y(a)-Y(b))
    Ведь исходя из формулы Симпсона значения начала отрезка и его конца СКЛАДЫВАЮТСЯ, а в нашем случае между ними знак минус ?

Добавить комментарий