Текст подпрограммы и версий
qsjac_p.zip  		Тексты тестовых примеров
tqsjac_p.zip 

Подпрограмма:  QSJAC (модуль QSJAC_p)

Назначение

Вычисление узлов и весов квадратурной формулы Гаусса - Якоби.

Математическое описание

Подпрограмма QSJАС вычисляет узлы  xi и веса Wi квадратурной формулы Гаусса - Якоби 

      b
     ∫ (1 - x)α (1 + x)β  f(x) dx    ≈
    a
                                                          N
                                                  ≈     ∑   Wi  f(xi )    по методу Ньютона .
                                                         i =1

Использование 

procedure QSJAC(var X :Array of Real; var W :Array of Real;
                ALF :Real; BTA :Real; N :Integer;
                A :Array of Real; B :Array of Real; C :Array of Real);


   Параметры
X - вещественный одномерный массив размерности 512 вычисленных узлов формулы Гаусса - Якоби;
W - вещественный одномерный массив размерности 512 вычисленных весов формулы Гаусса - Якоби;
         ALF -          BTA   заданные параметры  α  и  β , соответственно, в весовой функции (1 - x)α (1 + x) β (тип: вещественный);
N - заданное число узлов (и весов) интегрирования в формуле Гаусса - Якоби (2 ≤ N ≤ 512, тип: целый);
A, B, С - вещественные одномерные массивы размерности 512, заданных коэффициентов рекуррентного соотношения для полиномов Якоби. 
   Версии:  нет 
   Вызываемые подпрограммы
QSJAC1 - служебная подпрограмма.
QSJAC2 - вещественная подпрограмма - функция, вычисляющая значение логарифма от гамма - функции  ln ( Г (Х) ). 
   Замечания по использованию
 

Подпрограмма QSJAC вычисляет узлы и веса квадратурной формулы Гаусса - Якоби 

        N
        ∑   Wi  f(xi )
       i =1

с абсолютной погрешностью 10 - 10 и для значений 2 ≤ N ≤ 512 .

Перед обращением к подпрограмме QSJAC должны быть вычислены коэффициенты  Anα , β,  Bnα , β,  Cnα , β рекуррентного соотношения для полиномов Якоби 

  Jnα , β (x)   =   ( Anα , β  x  +  Bnα , β )  Jn -1α , β (x)  -  Cnα , β  Jn -2α , β (x)

  по формулам:
                             ( α + β + 2n )  ( α + β + 2n - 1 )
     Anα , β     =    -------------------------------------------     ,
                                   2n ( α + β + n )

                             ( α + β ) ( α - β ) ( α + β + 2n - 1 )
     Bnα , β     =    ------------------------------------------------     ,
                               2n ( α + β + n ) ( α + β + 2n - 2 )

                             ( α + β + 2n ) ( α + 2n - 1 ) ( β + 2n - 1 )
     Cnα , β     =    ------------------------------------------------------     .
                                  n ( α + β + n ) ( α + β + 2n - 2 )

Пример использования 

Unit TQSJAC_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, QSJAC_p;

function TQSJAC: String;

implementation

function TQSJAC: String;
var
NN,N,K :Integer;
ALF,BTA,SAB,DBA,P,FN,FN1,Q1,Q2,Q3,Q4,Q5 :Real;
X :Array [0..511] of Real;
W :Array [0..511] of Real;
A :Array [0..511] of Real;
B :Array [0..511] of Real;
C :Array [0..511] of Real;
label
_1,_2;
begin
Result := '';  { результат функции }
ALF := 1.0;
ВТА := 1.0;
NN := 8;
SАВ := ALF+BTA;
DВА := ALF-BTA;
P := SAB*DBA;
A[0] := 0.0;
B[0] := 0.0;
C[0] := 0.0;
for N:=2 to 512 do
 begin
  FN := N;
  FN1 := N-1;
  Q1 := 2.0*FN+SAB;
  Q2 := Q1-1.0;
  Q3 := SAB+FN;
  Q4 := FN*Q3;
  Q5 := Q4*(Q1-2.0);
  A[N-1] := Q1*Q2/(2.0*Q4);
  B[N-1] := P*Q2/(2.0*Q5);
  C[N-1] := Q1*(ALF+FN1)*(BTA+FN1)/Q5;
_1:
 end;
QSJAC(X,W,ALF,BTA,NN,A,B,C);
for K:=1 to NN do
 begin
_2:
  Result := Result + Format('%20.16f %20.16f ',
 [X[K-1],W[K-1]]) + #$0D#$0A;
 end;
UtRes('TQSJAC',Result);  { вывод результатов в файл TQSJAC.res }
exit;
end;

end.

Результаты:

      x1  =   0.9195339081  ,    w1  =  0.0205900956  , 
      x2  =   0.7387738651  ,    w2  =  0.1021477023  ,
      x3  =   0.4779249498  ,    w3  =  0.2253365549  ,
      x4  =   0.1652789576  ,    w4  =  0.3185923136  ,
      x5  = - 0.1652789576  ,    w5  =  0.3185923136  ,
      x6  = - 0.4779249498  ,    w6  =  0.2253365549  ,
      x7  = - 0.7387738651  ,    w7  =  0.1021477023  ,
      x8  = - 0.9195339081  ,    w8  =  0.0205900956  .