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

Листинг 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] учитывается симметрия.