Приближенное вычисление интеграла

          sum:=left_integ+middle_integ+right_integ;

          integral_abs:=integral_abs-    Abs(eastimate)+Abs(left_integ)+Abs(middle_integ)+Abs(right_integ);

            if (nest_level>1) and ((nest_level=max_level) or

    (Abs(sum- estimate)<=eps+eta*integral_abs)) then simpson3point:=sum

    else

    Begin

            If nest_level>highest_level then

            Inc(highest_level);

             Eps3:=0.577*eps;

             Eta3:=0.577*eta;

             Left_integ:=simpson3point(x0,dx3,left_integ,integral_abs,eps3,eta3, left,F1,F2);

    middle_integ:=simpson3point(x0+dx3,dx3,middle_integ, integral_abs,eps3,eta3,F2,middle,F4);

right_integ:=simpson3point(x0+2.0*dx3,dx3,right_integ,integral_abs,eps3,eta3,F4,F5,right);

    simpson3poin:=left_integ+middle_integ+right_integ;

    end;

    Dec(nest_level);

    End; {simpson3point}

    Begin {adaptive_simpson}

             nest_level:=1;

             highest_level:=1;

             no_evaluations:=3;

             adaptive_simpson:=simpson3point(x0,x1-x0,0.0,0.0,eps,eta,F(x0),F(x0+0.5*(x1-x0)),F(x1));

    end;{adaptive_simpson} 

    1.3. Метод Ромберга  и его реализация  на языке Pascal. 

    Интегрирование  следующим методом – методом Ромберга – основано на правиле трапеций, использующем кусочно-линейное приближение для интегрируемой функции. Основной факт относительно погрешности в методе трапеций следующий.

    Теорема. Пусть F(x) – гладкая функция на интервале [a,b], и этот интервал делится на т равных частей, каждая длиной h = . Пусть I(h) обозначает соответствующее приближение метода трапеций:

    I(h) =

,

    где f_i=F(a+jh) – значение интегрируемой функции в точке a+jh.

    Тогда

    

,

    Где a_k – некоторая константа.

    Основное здесь то, что погрешность в методе трапеций может быть выражена рядом по четным степеням шага интегрирования h. Построим таблицу значений T_ik:

                                                        

                                         

                                 

                                             

                                       

    В нулевой строке T_0k = I((b – a)/2k), так что T₀₀,T₀₁,… являются последовательными приближениями метода трапеций для интеграла, каждое с удвоенным по сравнению с предыдущим числом интервалов. Согласно приведенной выше теореме,

    

,

    где h = ((b – a)/2k.

    Отсюда  следует, что 

    

,

    Поэтому положим

    

.

    В общем случае строим j-ю строку таблицы Ромберга по формуле

    

,

    а оценка погрешности имеет вид

    

,

    где h = (b – a)/2k.

    Для работы понадобится не целая таблица, а только последняя вычисленная строка. Число точек выборки на каждом шаге удваивается. Обратите внимание на то, что функцию следует вычислять только в новых точках, которые являются средними точками предыдущих подынтервалов:

    F₀ + 2F₁ + 2F₂ + …+ 2F_2n-1 + F_2n =

    = (F₀ + 2F₂ + 2F₄ + …+ 2F_2n-2 + F_2n) + 2(F₁ + F₃ +…+F_2n-1). 

    Таким образом, чтобы модифицировать предыдущее приближение, необходимо вычислить сумму значений функции в новых средних точках. Это делается в цикле со счетчиком k. Метод Ромберга реализован в функции romberg (листинг1.5).

    Листинг 1.5. Функция romberg модуля integral

    Function romberg (F:real_fun;x0,x1,eps,eta:real; min,max:word):real;

      const

              abs_max=30;

    var

         p,dx,error,F_of_x0, F_of_x1, F_of_xk,

         roundoff_error,integral_abs,tolerance,

         previous_estimate,current_estimate,

         mid_sum, temp_sum, mid_sum_abs:real;

         table:array[0..abs_max] of real;

         j,n:word;

         k,r:longint;

        done:Boolean;

        denom:array[1..abs_max] of real;

    begin

           p:=1.0;

           for k:=1 to abs_max do

    begin

           p:=4.0*p;

          denom[k]:=1.0/(p-1.0);

    end;

    dx:=x1-x0;

    F_of_x0:=F(x0);

    F_of_x1:=F(x1);

    current_estimate:=0.0;

    previos_estimate:=0.0;

    done:=False;

    table[0]:=0.5*dx*(F_of_x0+F_of_x1);

    integral_abs:=0.5*Abs(dx)*(Abs(F_of_x0)+Abs(F_of_x1));

    n:=1;

    r:=1;

    repeat

            dx:=0.5*dx;

            mid_sum:=0.0;

            mid_sum_abs:=0.0;

            roundoff_error:=0.0;

           for k:=1 to r do

    begin

            F_of_xk:=F(x0+(2*k-1)*dx);

            mid_sum_abs:=mid_sum_abs+Abs(F_of_xk);

            F_of_xk:=F_of_xk+roundof_eroor;

           temp_sum:=mid_sum+F_of_xk;

           roundof_error:=(mid_sum-temp_sum)+F_of_xk;

           mid_sum:=temp_sum;

        if KeyPressed then

                          Halt;

    end;

             table[n]:=0.5*table[n-1]+dx*mid_sum;

             integral_abs:=0.5*integral_abs+Abs(dx)*mid_sum_abs;

             for j:=n-1 downto 0 do

             table[j]:=table[j+1]+denom[n-j]*(table[j+1]-table[j]);

              if n>=min then

    begin

             tolerance:=eta*integral_abs+eps;

             error:=Abs(table[0]-current_estimate)+Abs(current_estimate-previos_estimate);

             done:=(error<tolerance);

    end;

    Inc(n);

    done:=done or(n>max);

    previous_estimate:=current_estimate;

    current_estimate:=table[0];

    r:=r+r;

    until done;

              romberg:=current_estimate;

    end; 

    1.4. Метод Гаусса и его реализация на языке Pascal. 

    Теперь  перейдем к гауссовским квадратурам – семейству правил интегрирования, основанных на неравномерном разбиении основного интервала интегрирования. Вообще, метод гаусса с n точками точен для полиномов степени 2n – 1. В функции gauss3 (листинг1.6.) основной трехточечный алгоритм Гаусса применяется к каждой  из n равных частей интервала. Для интервала [-1,1] узлами квадратурной формулы являются нули полинома Лежандра третьей степени P₃ = (5x³ – 3x)/2, а коэффициенты выбираются специальным образом.

    Листинг 1.6. Функция gauss3 модуля integral

    Function gaus3(F:real_fun; x0,x1:real; n:word):real;

    var

         t,sum,x,z,dx:real;

         i,k:word;

         gzero,gweight:array[1..3] of real;

    procedure initialize_constants;

    var

         s,t:real;

         j:word;

    begin

             gzero[1]:=-sqrt(0.6);

             gzero[2]:=0.0;

             gzero[3]:=sqrt(0.6);

             gweight[1]:=5.0/9.0;

             gweight[2]:=8.0/9.0

             gweight[3]:=5.0/9.0;

             for j:=1 to 3 do

    begin

            gzero[j]:=0.5*(1.0+gzero[j]);

            gweight[j]:=0.5*gweight[j]);

    end;

    end; {initialize_constants}

    begin {gauss3}

             initialize_constants;

             dx:=(x1-x0)/n;

             x:=x0;

             sum:=0.0;

             for i:=0 to n-1 do

    begin

               t:=0.0;

               for k:=1 to 3 do

    begin

               z:=x+dx*gzero[k];

               t:=t+gweight[k]*F(z);

    end;

               sum:=sum+dx*t;

               x:=x+dx;

    end;

              gauss3:=sum;

    end;{gauss3} 

    Дадим некоторый обзор некоторых свойств  полиномов Лежандра. Рекурсивное определение полиномов Лежандра приводилось ранее в этой работе. Они образуют ортогональное (но не ортонормированное)семейство на промежутке [-1,1], то есть

    

,  m≠n,

    

.

    Величина  второго интеграла определяет нормировку для этих полиномов. Имеет место также следующее представление полиномов Лежандра:

    

 .

    Другая  явная формула:

    

.

    Приведем  несколько первых полиномов Лежандра:

    P₀(x) = 1,

    P₁(x) = x,

    

,

    

,

    

,

    

.

    Очевидно, что в общем случае полиномы Лежандра нечетной степени являются нечетными функциями, а четной степени – четными функциями.

    Нам требуется найти нули полинома P_n(x). Важно здесь то, что эти нули являются простыми и принадлежат открытому интервалу (-1,1). Таким образом,

    -1<x₁<x₂<…<x_n<1, P_n(x_j) = 0.

    Соответствующая формула гаусовского интегрировании (с остатком) имеет следующий вид:

    

.

    В этой формуле

    

,

    Где -1 < <1.

    Веса  задаются несколькими эквивалентными формулами

    

.

    Процедура compute_gauss_coeffs (листинг1.7). предназначена для вычисления нулей и весов квадратурной формулы Гаусса. Подпрограмма legendre_poly вычисляет значения P_n(x), P_n-1(x) и P_n’(x). Последнее получается дифференцированием основной рекуррентной формулы для P_n(x):

    

.

    Нули  находятся предварительным делением интервала и применением метода секущих, после чего следует ньютоновские итерации, в которых используются значения производных. Затем применяется первая формула для весов, в которой вновь используются значения производных. Здесь zero – массив нулей полинома Лежандра n-й степени, а weight – массив соответствующих весов.