Содержание
Введение                                                                                                  2

1. Различные методы вычисления определенных интегралов           3 

1.1. Метод Симпсона для интегрирования функций F(x) по

 заданному промежутку и его реализация на языке Pascal                  4

1.2. Метод Симпсона для интегрирования функции от двух

переменных F(x,y) по прямоугольной двумерной области и его

реализация на языке Pascal                                                                     5

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

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

Заключение                                                                                              16

Литература                                                                                               17


Введение
Система программирования Турбо Паскаль представляет собой единство двух в известной степени самостоятельных начал: компилятора с языка программирования Паскаль (язык назван в честь выдающегося французского математика и философа Блеза Паскаля (1623-1662)) и некоторой инструментальной программной оболочки, способствующей повышению эффективности создания программ.

Паскаль – гибкий и развитый в отношении типов данных язык. Привлекательны его рекурсивные возможности, а также поддержка технологии объектно-ориентированного программирования.

Изучение программирования на языке Паскаль может дать хороший старт в огромный и увлекательный мир программирования. Обучение языку программирования проходит намного более эффективно с изучением примеров.

В данной работе рассмотрен пример использования языка программирования высокого уровня Pascal для вычисления определенных интегралов.




1.     
Различные методы вычисления определенных интегралов.
Приближенное вычисление интеграла,

I = ,

Основано на его замене конечной суммой:

I_n = _k
F(x_k),

где w_k – числовые коэффициенты, а x_k – точки отрезка [x₀,x₁]. Приближенное равенство

I ≈ I_n называется квадратурной формулой, точки x_k – узлами квадратурной формулы, а числа w_k – коэффициентами квадратурной формулы. Разные методы приближенного интегрирования отличаются выбором узлов коэффициентов. От этого выбора зависит погрешность квадратурной формулы.

R_n =‌ ‌.

В модуле integral реализовано несколько методов численного интегрирования как для простых (одномерных), так и для кратных (многомерных) интегралов.

В функции simpson реализован стандартный метод Симпсона для интегрирования функции F(x) по заданному промежутку, когда число разбиений интервала выбирается заранее. Функция double_simpson является прямым обобщением метода Симпсона на случай интегрирования функции от двух переменных F(x,y) по прямоугольной двумерной области.

Функция adaptive_simpson служит для вычисления простых интегралов, она корректирует число и размер разбиений интервала, чтобы ошибка вычисления интеграла попала в заранее заданный интервал. Этот метод называется адаптивным интегрированием. Все современные программы интегрирования так или иначе адаптивны.

В функции romberg  запрограммирован еще один метод адаптивного интегрирования – метод Ромберга, в настоящее время, вероятно, один из наиболее популярных. Имеются также функция gauss – одномерная версия метода интегрирования Гаусса. Интерфейсная секция модуля integral приведена в листинге 1.1.

Листинг 1.1. Интерфейсная секция модуля integral.

Unit integral;

Interface

Const 

         Max_dim =10;

         Max_deg=96;

Type

Real_fun=function(x:real):real;

Real_fun2=function(x,y:real):real;

Real_vec=array[1..max_dim+1] of real;

Index=array[1..max_dim+1] of word;

Vec_fun=function(j:word; x:real_vec):real;

 Var

no_evaluations, highest_level:word;

function simpson(F:real_fun; x0,x1:real; div_no:word):real;

function double_simpson(F:real_fun2; x0,x1,y0,y1:real; x_div,y_div:word):real;

function adaptive_simpson(F:real_fun;x0,x1,eps,eta:real):real;

function romberg(f:real_fun; x0,x1,eps,eta:real; min,max:word):real;

function gauss3(F:real_fun;x0,x1:real; n:word):real;

procedure compute_gauss_coeffs(deg:word);

function gauss(Freal_fun:x0,x1:real; deg:word):real;
1.1. Метод Симпсона для интегрирования функций
F(
x) по заданному промежутку и его реализация на языке
Pascal.
Перейдем к секции реализации. Она начинается описанием функции simpson. Стоит сказать несколько слов о выборе узлов и коэффициентов квадратурной формулы Симпсона. Идея трехточечного метода Симпсона заключается в следующем.

Пусть x_m – это средняя точка интервала [x₀, x₁] и пусть Q(x) – единственный полином второй степени, который интерполирует (приближает) подынтегральную функцию F(x) по точкам x₀, x_m и x₁. Искомый интеграл аппроксимируется интегралом от функции  Q(x):

I≈.

Это оценка точна, если F(x) является полиномом степени 3.

В функции simpson интервал интегрирования делится на div_no равных частей, а трехточечная формула Симпсона применяется к каждому такому интервалу. Параметрами функции simpson (листинг1.2) являются, по порядку, подынтегральная функция, нижняя и верхняя границы интервала интегрирования и количество подынтервалов.
Листинг
1.2. Функция simpson модуля integral

Implementation

Uses crt;

Var

      zero, weight:array [1..max_deg] of real;

Function simpson(F:real_fun; x0,x1:real; div_no:word):real;

Var

      x, dx, sum:real;

      j:word;

begin

       dx:=(x1-x0)/(2.0*div_no);

        sum:=F(x0)+F(x1);

        x:=x0;

        for j:=1 to 2*div_no-1 do

begin

        x:=x+dx;

        if Odd (j) then

        sum:=sum+4.0*F(x)

else

        sum:=sum+2.0*F(x);

    end;

        simpson:=dx*sum/3.0;

end;

 

1.2. Метод Симпсона для интегрирования функции от двух переменных
F(
x,
y) по прямоугольной двумерной области и его реализация на языке
Pascal.
Функция double_simpson (листинг 1.3.) является, по существу, прямым обобщением одномерного метода Симпсона на случай вычисления двойного интеграла  по прямоугольной области.

Листинг
1.3. Функции double_simpson и simple_simpson модуля integral

Function double_simpson(F:real_fun2; x0,x1,y0,y1:real; x_div,y_div:word):real;

var

     dx,dy,x,sum:real;

      i:word;

function simple_simpson(x:real):real;

var

     y,sum:real;

     j,v:word;

begin

     sum:=F(x,y0)+F(x,y1);
      y:=y0;

      for j:=1 to 2*y_div-1 do

      begin

              y:=y+dy;

              if Odd(j) then

              sum:=sum+4.0*F(x,y)

       else

               sum:=sum+2.0*F(x,y);

end;

       simple_simpson:=sum;

end;{simple_simpson}

begin{doudle_simpson}

        dx:=(x1-x0)/(2.0*x_div);

        dy:=(y1-y0)/(2.0*y_div);

         x:=x0;

        sum:=simple_simpson(x0)+simple_simpon(x1);

        for i:=1 to 2*x_div-1 do

        begin

                x:=x+dx;

                if Odd(i) then

                         sum:=sum+4.0*simple_simpson(x)

         else

                         sum:=sum+2.0*simple_simpson(x);

end;

     double_simpson:=dx+dy*sum/9.0;

end;{double_simpson}
Недостатком рассмотренных функций интегрирования является то, что они не дают возможности явно задать точность вычисления интеграла. Точность связана с количеством точек разбиения, но ее значение в этих функциях не определяется с адаптивным выбором шага разбиения. Такой функцией является adaptive_simpson. Параметры eps и eta задают соответственно абсолютную и относительную погрешности. Их роль поясняется следующим неравенством:

.

Функция adaptive_simpson (листинг1.4) использует рекурсивную процедуру simpson3point, которая вычисляет значение интеграла по интервалу [x₀, x₀+δ
x], где x₀ – не обязательно исходная левая граничная точка.

Если трехточечный метод Симпсона не дает достаточную точность на данном интервале, этот интервал делится на три равные части, и метод вновь применяется к каждой из полученных частей. В результате получим 7 точек разбиения, но вычислять функцию F(x) придется только в четырех из них, поскольку значения в других трех точках уже известны.

При адаптивном разбиении имеется одна тонкость. При переходе к подынтервалам, составляющим одну треть от исходного, чтобы получить новые абсолютную и относительную погрешности, надо поделить eps и eta на .

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

Function adaptive_simpson(F:real_fun:x0,x1,eps,eta:real):real;

   const

            max_level=35;

var

        k,nest_level:word;

        integral_abs:real;

function simpson3poin(x0,delta_x, estimate, integral_abs,

                  eps,eta,left,middle,right:real):real;

var

     dx3,sum,eps3,eta3,factor,left_integ,

      middle_integ, right_integ,F1,F2,F4,F5:real;

begin

       Inc(nest_level);

       dx3:=delta_x/3.0;

       F1:=F(x0+0.5*dx3);

       F2:=F(x0+dx3);

       F:=F(x0+2.0*dx3);

       F5:=F(x0+2.5*dx3);

       Inc(no_evaluations,4);

       factor:=dx3/6.0;

       left_integ:=factor*(left+4.0*F1+F2);

       middle_integ:=factor*(F2+4.0*middle+F4);

      right_integ:=factor*(F4+4.0*F5+right);

      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 – массив соответствующих весов.

Метод вычислений нулей полинома заключается в том, чтобы поделить интервал [0,1] на маленькие подынтервалы и проверить каждый из них на изменение знака полинома. Если изменение знака имеет место, то однократное применение метода секущих позволяет достаточно хорошо определить положение нуля. Для уточнения этого значения применяется метод Ньютона. Для обработки интервала [-1,0] учитывается симметрия.

Листинг 1.7. Процедура compute_gauss_coeffs модуля integral

procedure compute_gauss_coeffs(deg:word);

     const

            eps=6.0e-20;

var

     i,index:word;

     P0k, P0k_1,D0k, P1k,P1k_1, D1k,

     x0,x1,y,z,dx,x,u:real;

procedure legendre_poly(n:word; x:real; var Pk,Pk_1,Dk:real);

var

      Pk_2,Dk_1,Dk_2:real;

      i,j,k:word;

begin

        If n=0 then

Begin

         Pk:=1.0;

         Dk:=0.0;

end

     else

begin

        Pk_1:=1.0;

        Pk:=x;

        Dk_1:=0.0;

        Dk:=1.0;

        i:==3;

        j:=1;

        for k:=2 to n do

begin

              Pk_2:=Pk_1;

              Pk_1:=Pk;

              Dk_2:=Dk_1;

              Dk_1:=Dk;

              Pk:=(i*x*Pk_1-j*Pk_2)/k;

              Dk:=(i*(Pk_1+x*Dk_1)-j*Dk_2)/k;

Inc(I,2);

Inc(j);

                     end;

            end

end;{legendre_poly}

begin{computr_gauss_coeffs}

         index:=(deg+1) div 2;

         dx:=1.0/(10.0*deg);

         x0:=0.0;

         x1:=x0+dx;

            if Odd(deg) then

begin

         zero[index]:=0.0;

         legendre_poly(deg,x0,P0k, P0k_1,D0k);

         weight[index]:=2.0/(P0k_1*D0k*deg);

end;

           for i:=0 to 10*deg-1 do

begin

          x0:=x1;

          x1:=x1+dx;

          legendre_poly(deg,x0,P0k_1,D0k);

          legendre_poly(deg,x1,P1k,P1k_1,D1k);

          if P0k*P1k<=0.0 then

begin

          x:=x0-P0k*dx/(P1k-P0k);

          legendre_poly(deg,x,P0k,P0k_1,D0k);

          u:=P0k/D0k;

          y:=x-u;

          while Abs(x-y)>=eps do

begin

         if keyPressed then begin

         writeln(‘>=eps loop:’,x:10:10,’ ‘, y:10:10,’ ‘,Abs(x-y):10);

    readln;

end;

       x:=y;

       legendre_poly(deg,x,P0k,P0k_1,D0k);

       u:=P0k/D0k;

       y:=x-u;

end;

      inc(index);

      legendre_poly(deg,y,P0k,P0k_1,D0k);

      zero[index]:=y;

      weight[index]:=2.0/(P0k_1*D0k*deg);

      if index=deg then

                 Break;

end;

end;

      For i:=1 to deg div 2 do

begin

         Zero[i]:=zero[deg-i+1];

         Weight[i]:=weight[deg-i+1];

end;

end;{compute_gauss_coeffs}
В функции gauss (листинг 1.8.) запрограммирован один гауссовский шаг на заданном интервале. Конечно, серьезная прикладная программа будет делить интервал на меньшие подынтервалы, применять эту процедуру к каждому из них адаптивным способом.

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

Function gauss(F:real_fun;x0,x1:real;deg:word):real;

var

    Index:word;

    a,b,sum:real;

begin

    a:=0.5*(x1-x0);

    b:=0.5*(x1+x0);

     sum:=0.0;

     for index:=1 to deg do begin

      sum:=sum+F(a*zero[index]+b)*weight[index];

              if KeyPressed then

                            Halt;

end;

      gauss:=a*sum;

end;

end.
Заключение
В данной работе были рассмотрены различные методы интегрирования определенных интегралов и их реализация на языке программирования высокого уровня Pascal. Таким образом было показано, что данный язык программирования возможно использовать для решения различных задач из области высшей математики и численных методов. В данной работе затронута лишь одна проблема – проблема вычисления интегралов, но Pascal позволяет решать и такие проблемы как: решение дифференциальных уравнений, вычисление с полиномами, решение нелинейных уравнений, вычисления связанные с линейной алгеброй.
Литература
1.      Немнюгин С.А.  Turbo Pascal, 2 - издание – С-П.: Питер, 2003-544 с.

2.      Большой энциклопедический словарь: под редакцией Ю.В. Прохорова – М.: Большая Российская энциклопедия, 2000. – 845 с.