Численное интегрирование
   
Интегрирование методом средних прямоугольников (с заданной точностью)  
Интегрирование методом трапеций (с заданной точностью)  
Интегрирование методом Симпсона (с заданной точностью)  
   
Интегрирование на основе квадратурной формулы Гаусса-Кронрода с 61 узлом
(с заданной точностью)		  


Интегрирование методом средних прямоугольников (с заданной точностью)

Собственное значение определенного интеграла 
    b
S = ∫ F(x) dx
    a
находится методом прямоугольников. Отрезок [a, b] разбивается на n частей x0=a, x1=a+h, ..., xn=b с шагом h = (b-a)/n. На каждом отрезке [xi, xi+1] вычисляется значение интеграла по формуле прямоугольников, таким образом: 
b               n-1
∫ F(x) dx = h * ∑ (F(xi) + h/2) + Rn
a              i = 0
Затем количество точек разбиения удваивается и производится оценка точности вычислений:
Rn = |S2n - Sn|

Если Rn > 3*ε, то количество точек разбиения удваивается. 
function F(x: double): double;
begin
  F := x * x * x * x + 2 * x * x + 4
end;

function IntegralRect(const a, b, Epsilon : Double): Double;
var
  i, n: Integer;
  h, s1, s2: Double;
begin
  n := 1;
  h := b - a;
  s2 := h * F((a+b)/2);
  repeat
    n := 2*n;
    s1 := s2;
    h := h/2;
    s2 := 0;
    i := 1;
    repeat
      s2 := s2+F(a+h/2+h*(i-1));
      i := i + 1;
    until not (i <= n);
    s2 := s2*h;
  until not (Abs(s2-s1) > 3*Epsilon);
  IntegralRect := s2;
end;

begin
  writeln('[Rectangle] result = ', IntegralRect(0, 10, 0.01):12:8);
end.
Получаем результат: 
[Rectangle] result = 20706.66265329
Этот метод является самым медленным методом интегрирования с заданной точностью.


Интегрирование методом трапеций (с заданной точностью)

Собственное значение определенного интеграла 
    b
S = ∫ F(x) dx
    a
находится методом трапеций. Отрезок [a, b] разбивается на n частей x0=a, x1=a+h, ..., xn=b с шагом h = (b-a)/n. На каждом отрезке [xi, xi+1] вычисляется значение интеграла по формуле трапеций, таким образом: 
b                 n-1
∫ F(x) dx = h/2 * ∑ (F(xi) + F(xi+1)) + Rn
a                i = 0
Затем количество точек разбиения удваивается и производится оценка точности вычислений:
Rn = |S2n - Sn| / 3

Если Rn > ε, то количество точек разбиения удваивается. Для уменьшения времени, затрачиваемого на вычисления в алгоритме при получении S2n используется формула: 
                    n-1
S2n = 1/2 Sn + h/2 * ∑ (f(a + ((2*i - 1)*h)/2)
                   i = 0

function F(x: double): double;
begin
  F := x * x * x * x + 2 * x * x + 4
end;

function IntegralTrap(const a, b, Epsilon: double): double;
var
  i, n: Integer;
  h, s1, s2: Double;
begin
  n := 1;
  h := b-a;
  s2 := h*(F(a)+F(b))/2;
  repeat
    s1 := s2;
    s2 := 0;
    i := 1;
    repeat
      s2 := s2+F(a-h/2+h*i);
      i := i+1;
    until not (i <= n);
    s2 := s1/2 + s2*h/2;
    n := 2*n;
    h := h/2;
  until not (Abs(s2-s1) > 3*Epsilon);
  IntegralTrap := s2;
end;

begin
  writeln('[Trapeze] result = ', IntegralTrap(0, 10, 0.01):12:8);
end.
Получаем результат: 
[Trapeze] result = 20706.67469343
Метод трапеций хоть и работает чуть быстрее метода прямоугольников, но все равно является очень медленным.


Интегрирование методом Симпсона (с заданной точностью)

Собственное значение определенного интеграла 
    b
S = ∫ F(x) dx
    a
находится методом Симпсона (парабол). Отрезок [a, b] разбивается на n=2m частей x0=a, x1=a+h, ..., xn=b с шагом h = (b-a)/n.

Вычисляются значения yi = F(xi) функции в точках xi и находится значение интеграла по формуле Симпсона:
S = Sn + Rn
, где

где 

                              m-1          m-1
Sn = h/3 * [y0 + y2n + 2 * ∑ y2i + 4 * ∑ y2i+1]
                              i=1         i=0

              2*m
Rn = - h5/90 * ∑ yIV (ξk), ξk in [xk-1, xk]
              k=1
Затем количество точек разбиения удваивается и производится оценка точности вычислений:

Rn = |S2n-Sn| / 15
Если Rn > ε, то количество точек разбиения удваивается.

Значение суммы 2(y1+y2+...+y2m-1) сохраняется, поэтому для вычисления интеграла при удвоении количества точек разбиения требуется вычислять значения yi лишь в новых точках разбиения. 
function F(x: double): double;
begin
  F := x * x * x * x + 2 * x * x + 4
end;

function IntegralSimps(const a, b, Epsilon: double): double;
var
  i, n: Integer;
  h, s, s1, s2, s3: Double;
  x: Double;
begin
  s2 := 1;
  h := b-a;
  s := F(a)+F(b);
  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 > Epsilon);
  IntegralSimps := s2;
end;

begin
  writeln('[Simpson] result = ', IntegralSimps(0, 10, 0.01):12:8);
end.
Получаем результат: 
[Simpson] result = 20706.66746140
Метод Симпсона является средним по скорости, работает быстрее первых двух методов (прямоугольников и трапеций), но существуют и гораздо более быстрые методы интегрирования.


Интегрирование на основе квадратурной формулы Гаусса-Кронрода с 61 узлом (с заданной точностью)

Недостатком приведенных выше методов (прямоугольников, трапеций и Симпсона) является невысокая скорость работы. Они достаточно просты в использовании, но если попытаться найти с их помощью интеграл с относительной точностью 10-10, то при использовании первых двух методов шанс дождаться конца расчетов крайне мал. Метод Симпсона работает чуть быстрее, но иногда не хватает и той скорости, которую обеспечивает он.

В качестве альтернативного метода может использоваться квадратурная формула Гаусса-Кронрода.

Функция интегрируется на указанном отрезке [A, B] с помощью этой формулы, и мы получаем значение интеграла и оценку погрешности (как правило, взятую с огромным запасом). Если погрешность оказывается слишком велика, то мы разбиваем отрезок пополам, и интегрируем с использованием той же формулы на каждой из половинок отрезка. В случае, если суммарная погрешность по-прежнему велика, то мы делим пополам отрезок с максимальной погрешностью интегрирования. При необходимости процесс повторяется до момента, когда будет достигнута нужная точность.

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

Замечания:
  1. Требуемый размер кучи зависит от передаваемой функции, размера области интегрирования и желаемой точности ответа. Скажем, чтобы проинтегрировать осциллирующую функцию sin(1000*x) на отрезке [0, 1] с абсолютной точностью 10 -6 требуется разбить отрезок интегрирования на 16 суботрезков, а для вычисления интеграла от sin(10000*x) отребуется 128 суботрезков. Более простые функции могут вообще не потребовать разбивать область. Разумным выбором будет начать с размера кучи в несколько сотен элементов.
  2. Следует обратить внимание на то, что алгоритм анализирует не относительную, а абсолютную погрешность. При этом относительная точность, с которой ответ может быть найден, не будет превосходить машинную точность. Скажем, если величина интеграла равна 1020 то на большинстве машин сосчитать его с абсолютной погрешностью 10-3 будет просто невозможно.
  3. Данный алгоритм выводит значение интеграла в любом случае, даже если он не удовлетворяет требованиям к точности. При этом в случае, подобном предыдущему, алгоритм просто переполнит кучу любого размера, рассчитав интеграл с максимально доступной точностью, после чего вернет код ошибки.
(программа немного изменена по сравнению с оригиналом, для обеспечения ее работоспособности при использовании компилятора Турбо Паскаль 7.0): 
function F(x: double): double;
begin
  F := x * x * x * x + 2 * x * x + 4
end;

const
  HeapSize = 100;
  HeapW    = 4;
  N        = 61;

type
  heapType = array[0 .. heapsize - 1, 0 .. heapw] of double;
  arrType  = array[0 .. pred(n)] of double;

(*
  Удаление из кучи самого большого элемента.
  
  Алгоритм перемещает элемент за пределы кучи (в элемент
  с номером HeapSize-1), затем вновь упорядочивает кучу.
*)
procedure MHeapPop(var Heap: heapType;
                   HeapSize, HeapWidth: Integer);
var
  I, p: Integer;
  T: Double;
  maxCP : Integer;
begin
  if HeapSize=1 then Exit;

  I := 0;
  while I <= HeapWidth-1 do begin
    T := Heap[HeapSize-1, I];
    Heap[HeapSize-1, I] := Heap[0, I];
    Heap[0, I] := T;
    Inc(I);
  end;

  P := 0;
  while 2*P+1 < HeapSize-1 do begin
    maxCP := 2*P+1;
    if 2*P+2 < HeapSize-1 then begin

      if Heap[2*P+2,0] > Heap[2*P+1,0] then MaxCP := 2*P+2;

    end;

    if Heap[P,0] < Heap[MaxCP,0] then begin
      I := 0;
      while I <= HeapWidth-1 do begin
        T := Heap[P, I];
        Heap[P, I] := Heap[MaxCP, I];
        Heap[MaxCP, I] := T;
        Inc(I);
      end;
      P := MaxCP;
    end
    else break;
  end;
end;

(*
  Добавление элемента в кучу.
  
  Алгоритм добавляет элемент с номером HeapSize,
  затем вновь упорядочивает кучу.
*)
procedure MHeapPush(var Heap: heapType;
          HeapSize, HeapWidth: Integer);
var
  I, p: Integer;
  T: Double;
  Parent: Integer;
begin
  if HeapSize = 0 then Exit;

  P := HeapSize;
  while P <> 0 do begin

    Parent := (P-1) div 2;
    if Heap[P,0] > Heap[Parent,0] then begin
      I := 0;
      while I <= HeapWidth-1 do begin
        T := Heap[P, I];
        Heap[P, I] := Heap[Parent, I];
        Heap[Parent, I] := T;
        Inc(I);
      end;
      P := Parent;
    end
    else Break;

  end;
end;

function IntegralGauss61(const A, B, Eps: Double;
         var Integral: Double; var HeapUsed: Integer): Boolean;

var
  Heap: heapType;
  SumErr : Double;
  X, WG, WK: arrType;

  NG, i, j: Integer;
  H: Integer;
  V: Double;
  C1, C2: Double;
  IntG, IntK: Double;
  TA, TB: Double;
begin

  NG := 15;
  I:=0;
  while I <= N-1 do begin
    X[I] := 0; WK[I] := 0; WG[I] := 0;
    Inc(I);
  end;

  WG[0] := 0.007968192496166605615465883474674;
  WG[1] := 0.018466468311090959142302131912047;
  WG[2] := 0.028784707883323369349719179611292;
  WG[3] := 0.038799192569627049596801936446348;
  WG[4] := 0.048402672830594052902938140422808;
  WG[5] := 0.057493156217619066481721689402056;
  WG[6] := 0.065974229882180495128128515115962;
  WG[7] := 0.073755974737705206268243850022191;
  WG[8] := 0.080755895229420215354694938460530;
  WG[9] := 0.086899787201082979802387530715126;
  WG[10] := 0.092122522237786128717632707087619;
  WG[11] := 0.096368737174644259639468626351810;
  WG[12] := 0.099593420586795267062780282103569;
  WG[13] := 0.101762389748405504596428952168554;
  WG[14] := 0.102852652893558840341285636705415;

  X[ 0] := 0.999484410050490637571325895705811;
  X[ 1] := 0.996893484074649540271630050918695;
  X[ 2] := 0.991630996870404594858628366109486;
  X[ 3] := 0.983668123279747209970032581605663;
  X[ 4] := 0.973116322501126268374693868423707;
  X[ 5] := 0.960021864968307512216871025581798;
  X[ 6] := 0.944374444748559979415831324037439;
  X[ 7] := 0.926200047429274325879324277080474;
  X[ 8] := 0.905573307699907798546522558925958;
  X[ 9] := 0.882560535792052681543116462530226;
  X[10] := 0.857205233546061098958658510658944;
  X[11] := 0.829565762382768397442898119732502;
  X[12] := 0.799727835821839083013668942322683;
  X[13] := 0.767777432104826194917977340974503;
  X[14] := 0.733790062453226804726171131369528;
  X[15] := 0.697850494793315796932292388026640;
  X[16] := 0.660061064126626961370053668149271;
  X[17] := 0.620526182989242861140477556431189;
  X[18] := 0.579345235826361691756024932172540;
  X[19] := 0.536624148142019899264169793311073;
  X[20] := 0.492480467861778574993693061207709;
  X[21] := 0.447033769538089176780609900322854;
  X[22] := 0.400401254830394392535476211542661;
  X[23] := 0.352704725530878113471037207089374;
  X[24] := 0.304073202273625077372677107199257;
  X[25] := 0.254636926167889846439805129817805;
  X[26] := 0.204525116682309891438957671002025;
  X[27] := 0.153869913608583546963794672743256;
  X[28] := 0.102806937966737030147096751318001;
  X[29] := 0.051471842555317695833025213166723;
  X[30] := 0.000000000000000000000000000000000;

  WK[ 0] := 0.001389013698677007624551591226760;
  WK[ 1] := 0.003890461127099884051267201844516;
  WK[ 2] := 0.006630703915931292173319826369750;
  WK[ 3] := 0.009273279659517763428441146892024;
  WK[ 4] := 0.011823015253496341742232898853251;
  WK[ 5] := 0.014369729507045804812451432443580;
  WK[ 6] := 0.016920889189053272627572289420322;
  WK[ 7] := 0.019414141193942381173408951050128;
  WK[ 8] := 0.021828035821609192297167485738339;
  WK[ 9] := 0.024191162078080601365686370725232;
  WK[10] := 0.026509954882333101610601709335075;
  WK[11] := 0.028754048765041292843978785354334;
  WK[12] := 0.030907257562387762472884252943092;
  WK[13] := 0.032981447057483726031814191016854;
  WK[14] := 0.034979338028060024137499670731468;
  WK[15] := 0.036882364651821229223911065617136;
  WK[16] := 0.038678945624727592950348651532281;
  WK[17] := 0.040374538951535959111995279752468;
  WK[18] := 0.041969810215164246147147541285970;
  WK[19] := 0.043452539701356069316831728117073;
  WK[20] := 0.044814800133162663192355551616723;
  WK[21] := 0.046059238271006988116271735559374;
  WK[22] := 0.047185546569299153945261478181099;
  WK[23] := 0.048185861757087129140779492298305;
  WK[24] := 0.049055434555029778887528165367238;
  WK[25] := 0.049795683427074206357811569379942;
  WK[26] := 0.050405921402782346840893085653585;
  WK[27] := 0.050881795898749606492297473049805;
  WK[28] := 0.051221547849258772170656282604944;
  WK[29] := 0.051426128537459025933862879215781;
  WK[30] := 0.051494729429451567558340433647099;

  I := N-1;
  while I >= N div 2 do begin
    X[I] := -X[N-1-I];
    Dec(I);
  end;
  I := N-1;
  while I >= N div 2 do begin
    WK[I] := WK[N-1-I];
    Dec(I);
  end;
  I := NG-1;
  while I >= 0 do begin
    WG[N-2-2*I] := WG[I];
    WG[1+2*I] := WG[I];
    Dec(I);
  end;

  I := 0;
  while I <= N div 2 do begin
    WG[2*I] := 0;
    Inc(I);
  end;
  C1 := 0.5*(B-A);
  C2 := 0.5*(B+A);
  IntG := 0;
  IntK := 0;
  I:=0;
  while I <= N-1 do begin
    V := F(C1*X[I]+C2);
    IntK := IntK+V*WK[I];
    if I mod 2 = 1 then IntG := IntG + V*WG[I];
    Inc(I);
  end;
  IntK := IntK*(B-A)*0.5;
  IntG := IntG*(B-A)*0.5;
  Heap[0, 0] := Abs(IntG-IntK);
  Heap[0, 1] := IntK;
  Heap[0, 2] := A;
  Heap[0, 3] := B;

  SumErr := Heap[0,0];
  if SumErr < Eps then begin
    IntegralGauss61 := True;
    Integral := IntK;
    HeapUsed := 1;
    Exit;
  end;

  HeapUsed := 1;
  H := 1;
  while H <= HeapSize-1 do begin
    HeapUsed := H+1;
    MHeapPop(Heap, H, HeapW);
    SumErr := SumErr-Heap[H-1,0];
    TA := Heap[H-1,2];
    TB := Heap[H-1,3];
    Heap[H-1,2] := TA;
    Heap[H-1,3] := 0.5*(TA+TB);
    Heap[H,2] := 0.5*(TA+TB);
    Heap[H,3] := TB;
    J:=H-1;
    while J <= H do begin
      C1 := 0.5*(Heap[J,3]-Heap[J,2]);
      C2 := 0.5*(Heap[J,3]+Heap[J,2]);
      IntG := 0;
      IntK := 0;
      I:=0;
      while I <= N-1 do begin
        V := F(C1*X[I]+C2);
        IntK := IntK+V*WK[I];
        if I mod 2 = 1 then IntG := IntG+V*WG[I];
        Inc(I);
      end;

      IntK := IntK*(Heap[J,3]-Heap[J,2])*0.5;
      IntG := IntG*(Heap[J,3]-Heap[J,2])*0.5;
      Heap[J,0] := Abs(IntG-IntK);
      Heap[J,1] := IntK;
      SumErr := SumErr+Heap[J,0];
      Inc(J);
    end;
    MHeapPush(Heap, H-1, HeapW);
    MHeapPush(Heap, H, HeapW);
    if SumErr < Eps then Break;
    Inc(H);
  end;
  IntegralGauss61 := (SumErr < Eps);
  Integral := 0;
  J := 0;
  while J <= HeapUsed-1 do begin
    Integral := Integral+Heap[J,1];
    Inc(J);
  end;
end;

var
  isOk: boolean;
  res: double;
  used: integer;
begin
  isOk := IntegralGauss61(0, 10, 0.0000000001, res, used);
  if isOk then writeln('result = ', res:12:8, ' heap used: ', used)
  else writeln('not enough heap to perform operation.');
end.
При запуске программы получаем: 
result = 20706.66666667 heap used: 1




Free Web Hosting