// (c) Copydown 2002, all left reserved. http://world.fpm.kubsu.ru.
 
{$APPTYPE CONSOLE}
 
{$F+} {разрешение передачи функций, как параметров}
 
type FunctionType = function(x: real): real;
 
{интегрируемая функция}
function f(x: real): real; begin f := x end;
 
{интегрирование от a до b функции f с точностью e}
function IntegralSimpson(a, b: real; f: FunctionType; e: real): real;
var
   h, x, s, s1, s2, s3, sign: real;
begin
 
if (a = b) then
   begin
     IntegralSimpson := 0; exit
   end;
 
if (a > b) then
   begin
     x := a; a := b; b := x; sign := -1
   end
  else sign:=1;
 
h := b - a; s := f(a) + f(b); s2 := s;
 
repeat
   s3 := s2; h := h/2; s1 := 0; x := a + h;
 
   repeat
     s1 := s1 + 2*f(x); x := x + 2*h;
   until (not(x < b));
 
   s := s + s1; s2 := (s + s1)*h/3; x := abs(s3 - s2)/15
until (not(x > e));
 
IntegralSimpson := s2*sign;
end;
 
begin
{вывод результата интегрирования от 0 до 1 функции f с точностью 0.001}
writeln(IntegralSimpson(0, 1, f, 0.001));
writeln; writeln('Press Enter'); readln;
end.

{**** UBPFD *********** by kladovka.net.ru ****
>> Взятие интеграла методом Симпсона
 
Интеграл методом Симпсона
A, B - интервал интегрирования
N - число точек на интервале
Func - функция, от которой берется интеграл.
 
Возвращаемое значение - значение интеграла
 
Зависимости: System
Автор:       Mystic, mystic2000@newmail.ru, ICQ:125905046, Харьков
Copyright:   Mystic
Дата:        25 апреля 2002 г.
********************************************** }
 
type
TFunction = function(X: Extended; Arg: Pointer): Extended;
 
function Simpson(A, B: Extended; N: Cardinal; Func: TFunction; Arg: Pointer): Extended;
var
h: Extended;
X: Extended;
K: Extended;
I: Integer;
begin
Assert(N>0);
h := 0.5 * (B-A) / N;
Result := Func(A, Arg);
X := A + h;
for I := 1 to 2*N-1 do
begin
   if I mod 2 = 0
     then K := 2
     else K := 4;
   Result := Result + K*Func(X, Arg);
   X := X + h;
end;
Result := Result + Func(B, Arg);
Result := h * Result / 3;
end;

{ **** UBPFD *********** by kladovka.net.ru ****
>> Вычисление определённого интеграла методом Симпсона
 
A, B - границы интегрирования
Eps - заданная относительная точность вычисления
F - подинтегральная функция
 
Зависимости: нет
Автор:       Dimka Maslov, mainbox@endimus.ru
Copyright:   Dimka Maslov
Дата:        26 ноября 2003 г.
********************************************** }
 
type
TDoubleFunc = function (X: Double): Double;
 
function Integral(A, B, Eps: Double; F: TDoubleFunc): Double;
 
function InternalCalc(A, B: Double; F: TDoubleFunc; N: Integer): Double;
var
x, dx: Double;
i: Integer;
begin
dx := (B - A) / N;
Result := 0;
x := A;
for i := 1 to N do begin
  Result := Result + dx * (F(x) + 4*F(x+dx/2) + F(x+dx)) / 6;
  x := x + dx;
end; 
end;
 
var
N: Integer;
Prev: Double;
begin
Result := InternalCalc(A, B, F, 4);
N := 4;
repeat
Prev := Result;
N := N shl 1;
Result:= InternalCalc(A, B, F, N);
until (Result = 0) or (Abs((Result-Prev) / Result) < Eps);
end;

function F(X: Double): Double;
begin
Result := X * X * X;
end;
 
procedure TForm1.Button1Click(Sender: TObject);
begin
Label1.Caption := FloatToStr(Integral(-10, 10, 0.00001, F));
end;

{ **** UBPFD *********** by kladovka.net.ru ****
>> Вычисление определенного интеграла методом левых и правых прямоугольников с заданной точностью
 
"Просто расчет площади под функцией, параметры: a,b - пределы интегрирования, a<=b
eps - допустимая погрешность, практически гарантируется, что расхождение результата с истинным значением интеграла не
превосходит по модулю указанную величину. Только не переборщите :-))
intF - подинтегральная функция. Естественно, желательно задавать функции, интегрируемые в смысле Римана. Объявление
смотри в примере.
Примечание: Несобственные интегралы не считаем :-)
Проверок на переполнение нет, да и вообще нет проверок..."
(Romkin (Москва))
 
Модуль сделан на основе функции вычисления опред. интеграла методом трапеций от Romkin'а (Москва).
 
Зависимости: Нет
Автор:       Алексей Глеб, noodlesf@mail.ru, Чернигов
Copyright:   с подачи Romkin'а (Москва)
Дата:        18 мая 2003 г.
********************************************** }
 
Unit IntPram;
 
Interface
 
Type
TIntFunc=Function(X: Double): Double;
 
Function LeftPramInt(a, b: Double; eps: Double; IntF: TIntFunc): Double;
Function RightPramInt(a, b: Double; eps: Double; IntF: TIntFunc): Double;
 
Implementation
 
Function LeftPramInt(a, b: Double; eps: Double; IntF: TIntFunc): Double;
Var
//S - площадь на предыдущей итерации,
//step - "толщина" прямоугольника
//gran - передвигаемая от a до b граница
//n - число прямоугольников, удваивается на каждой итерации
S, step, gran: Double;
n: integer;
Begin
//Сначала приближение одного прямоугольника
step:=b-a;
Result:=IntF(a)*step;
n:=1;
Repeat
   S:=Result;
   n:=n*2;
   step:=(b-a)/n;
   Gran:=a;
   Result:=0;
   //Ниже - просто вычисляем площади новых прямоугольников
   while gran<b do
   Begin
     Result:=Result+IntF(gran)*step;
     gran:=gran+step;
   End;
Until abs(S-Result)<=eps;
End;
 
Function RightPramInt(a, b: Double; eps: Double; IntF: TIntFunc): Double;
Var
//S - площадь на предыдущей итерации,
//step - "толщина" прямоугольника
//gran - передвигаемая от a до b граница
//n - число прямоугольников, удваивается на каждой итерации
S, step, gran: Double;
n: integer;
Begin
//Сначала приближение одного прямоугольника
step:=b-a;
Result:=IntF(b)*step;
n:=1;
Repeat
   S:=Result;
   n:=n*2;
   step:=(b-a)/n;
   Gran:=b;
   Result:=0;
   //Ниже - просто вычисляем площади новых прямоугольников
   while a<gran do
   Begin
     Result:=Result+IntF(gran)*step;
     gran:=gran-step;
   End;
Until abs(S-Result)<=eps;
End;
 
End.

uses IntPram;
 
function IntSqrt(x: Double): Double;
begin
Result:=Sqrt(x);
end;
 
procedure TForm1.Button1Click(Sender: TObject);
begin
label1.Caption:=FloatToStr(LeftPramInt(0, Pi, 0.00001, S));
label2.Caption:=FloatToStr(RightPramInt(0, Pi, 0.00001, S));
end;

{ **** UBPFD *********** by kladovka.net.ru ****
>> Вычисление определенного интеграла методом трапеций с заданной точностью
 
Просто расчет площади под функцией, параметры: a,b - пределы интегрирования, a<=b
eps - допустимая погрешность, практически гарантируется, что расхождение результата с истинным значением интеграла не
превосходит по модулю указанную величину. Только не переборщите :-))
intF - подинтегральная функция. Естественно, желательно задавать функции, интегрируемые в смысле Римана. Объявление
смотри в примере.
Примечание: Несобственные интегралы не считаем :-)
Проверок на переполнение нет, да и вообще нет проверок...
 
Зависимости: Ну какие могут быть зависимости?
Автор:       Romkin, romkin@pochtamt.ru, Москва
Copyright:   Romkin 2002
Дата:        19 ноября 2002 г.
********************************************** }
 
unit intfunc;
 
interface
 
type
TIntFunc = function(X: Double):Double;
 
function TrapezeInt(a,b:Double; eps: Double; IntF: TIntFunc): Double;
 
implementation
 
function TrapezeInt(a,b:Double; eps: Double; IntF: TIntFunc): Double;
var
//S - площадь на предыдущей итерации,
//x - текущее значение аргумента
//base - высота трапеции
//n - число трапеций, удваивается на каждой итерации
S, x, base: Double;
i, n: Integer;
begin
//Сначала приближение одной трапецией
base := b-a; 
Result := (IntF(a) + IntF(b))/2 * base;
eps := eps / 10; //Вообще говоря, величина делителя зависит от функции
n := 1;
repeat
   S := Result; 
   base := base / 2;
   n := n * 2;
   //Новая площадь вычисляется на основе старой
   Result := Result / 2;
   //Ниже - просто вычисляем площади новых трапеций
   for i := 1 to n div 2 do
   begin
     x := a + base * (i * 2 - 1);
     Result := Result + IntF(x) * base;
   end;
until abs(S-Result) <= eps;
end;
 
end.

uses intFunc;
 
function intSin(x: Double): Double;
begin
Result := sin(x);
end;
 
procedure TForm1.Button1Click(Sender: TObject);
begin
Edit1.Text := FloatToStr(TrapezeInt(0,Pi, 0.00001, intSin));
//результат у меня получился 1.99999990195429 - с запасом
//Точный ответ - 2.0
end;