Текст подпрограммы и версий
id10r_p.zip , id10e_p.zip 		Тексты тестовых примеров
tid10r_p.zip , tid10e_p.zip

Подпрограмма:  ID10R (модуль ID10R_p)

Назначение

Вычисление первой или второй производной вещественной таблично - заданной функции на равномерной сетке с использованием центральных разностей.

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

ID10R вычисляет либо первые, либо вторые производные таблично - заданной функции на равномерной сетке с шагом  h с использованием центральных разностей по формулам Стирлинга.

И.С.Березин, Н.П.Жидков, Методы вычислений, т.1, M., 1962.

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

procedure ID10R(H :Real; var Y :Array of Real; MU :Integer;
                NL :Integer; NS :Integer; NU :Integer; NO :Integer;
                LO :Integer; NTH :Integer; T :Real;
                var A :Array of Real; var E :Array of Real;
                var LA :Integer; var IERR :Integer);


   Параметры
H - заданный шаг равномерной сетки (тип: вещественный);
Y - вещественный двумерный массив, в котоpом задается таблица значений функции и центральных разностей (см. описание параметра LO); размерность массива Y должна быть при LO = 0 pавна MU на (NO + 1), при LO = 1 - MU на [(NO +  3)/2], при LO = 2 - MU на [(NO + 4)/2];
MU - заданное число стpок в описании размерности двумерного массива Y в программе, вызывающей данную программу (тип: целый);
         NL, NS -             NU   определяют расположение в массиве Y заданных значений функций, а именно, значения функции  y0, y1, y2, ..., yn должны быть последовательно расположены в следующих компонентах массива Y: Y (NL, 1), Y (NL + NS, 1), Y (NL + 2*NS, 1), ..., Y (NU, 1); расположение в массиве Y центральных разностей определяется ниже при описании параметра LO (тип: целый);
NO - заданный максимальный порядок разностей, которые используются в конечно - разностных аппроксимациях производных (тип: целый);
LO - определяет порядок расположения заданных центральных разностей в массиве Y (тип: целый); при этом если:
LO = 0 - то δ 2j - 1yi - 1/2 должна быть расположена в элементе массива Y (NL + I* NS - [NS/2], 2*J), а δ 2jyi - в элементе массива Y (NL + I*NS, 2*J + 1);
LO = 1 - то δ 2j - 1yi - 1/2 должна быть расположена в элементе массива Y (NL + I* NS - [NS/2], J + 1), а δ 2jyi - в элементе массива Y (NL + I*NS, J + 1);
LO = 2 - то δ 2j - 1yi - 1/2 должна быть расположена в элементе массива Y(NL + I* NS - [NS/2], J + 1), а δ 2jyi - в элементе массива Y (NL + I*NS, J + 2);
NHT - полагается равным 1, если надо вычислить первую производную, или 2, если - вторую;
T - заданная точность вычисления производных (тип: вещественный);
A - вещественный вектоp длины NU, в котоpом помещаются вычисленные значения первой или второй производной последовательно в следующем порядке A (NL), A (NL + NS), A (NL + 2*NS), ... ,A (NU);
E - вещественный вектоp длины NU, в котоpом помещаются последние члены ряда, использованные в конечно - разностной аппроксимации, в следующем порядке E (NL), E (NL + NS), E (NL + 2*NS),...,E (NU);
LA - целая переменная, значение которой полагается равным наивысшему порядку разностей, фактически использованных в вычислениях;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом:
IERR= 1 - когда для достижения заданной точности требуются разности более высокого порядка, чем они заданы при обращении к подпрограмме; производные вычисляются, хотя требуемая точность не достигнута;

IERR= 2 - когда разности порядка NO pасходятся; производные вычисляются, хотя требуемая точность не достигнута;
IERR= 3 - когда разности порядка LA (LA < NO) расходятся; производные вычисляются, хотя требуемая точность не достигнута;
IERR=65 - когда значение NTH ≠ 1 или ≠ 2;
IERR=66 - когда неправильно задан порядок разностей; порядок разностей должен быть по крайней меpе, не ниже порядка вычисляемой производной;
IERR=67 - когда не заданы разности нечетного порядка и требуется вычислить первую производную;
IERR=68 - когда не заданы разности четного порядка и требуется вычислить вторую производную. 
   Версии
ID10E - вычисление первой или второй производной вещественной таблично - заданной функции на pавномеpной сетке с использованием центральных разностей с расширенной (Extended) точностью. 
   Вызываемые подпрограммы
UTID10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы ID10R.
UTID11 - подпрограмма выдачи диагностических сообщений при работе подпрограммы ID10E. 
   Замечания по использованию
  1. 

Для ID10E параметры H, Y, T, A и E должны иметь тип Extended.

  2. 

Значение параметра MU должно быть больше или pавно значению параметра NU.

  3. 

Формула вычисления первой производной использует разности только нечетных порядков. Поэтому можно задать LO = 2 и NS = 1.

  4. 

Формула вычисления второй производной использует разности только четных порядков. Поэтому можно задать LO = 1 и NS = 1.

  5. 

Заданная точность считается достигнутой, если следующий член ряда конечно - разностной аппроксимации меньше или pавен T.

  6. 

Конечно - разностные аппроксимации не могут быть использованы в конце таблицы задания функции. Поэтому следует увеличить таблицу значений функции и ее разностей перед обращением к подпрограмме для вычисления производных вблизи конца таблицы.

  7.  Для вычисления соответствующего числа центральных разностей рекомендуется использовать подпрограмму I I10R.

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

Unit TID10R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, II10R_p, ID10R_p;

function TID10R: String;

implementation

function TID10R: String;
var
NTH,I,J,MU,NL,NS,NU,LO,NO,LA,IERR,_i :Integer;
H,XN,X,T :Real;
Y :Array [0..49] of Real;
A :Array [0..9] of Real;
E :Array [0..9] of Real;
label
_3,_1,_2;
begin
Result := '';  { результат функции }
for _i:=0 to 9 do //начальное обнуление массива
 begin
  A[_i] := 0.0;
  E[_i] := 0.0;
 end;
H := 0.25;
XN := 2.0;
X := XN-H;
NТН := 1;
T := 1.E-4;
for I:=1 to 10 do
 begin
  for J:=1 to 5 do
   begin
_3:
    Y[(I-1)+(J-1)*10] := 0.0;
   end;
 end;
for I:=1 to 10 do
 begin
  X := X+H;
  Y[(I-1)+0] := Exp(-X);
_1:
 end;
MU := 10;
NL := 1;
NS := 1;
NU := 10;
LO := 2;
NO := 6;
II10R(MU,NL,NS,NU,NO,LO,Y);
Result := Result + #$0D#$0A;
for J:=1 to 5 do
 begin
  for I:=1 to 1 do
   begin
    Result := Result + Format('%20.16f ',[Y[(I-1)+(J-1)*10]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
ID10R(H,Y,MU,NL,NS,NU,NO,LO,NTH,T,A,E,LA,IERR);
for I:=1 to 10 do
 begin
  Result := Result + Format('%s',['       A']);
  Result := Result + Format('%s',['               E' + #$0D#$0A]);
  Result := Result + Format(' %20.16f %20.16f ',
 [A[I-1],E[I-1]]) + #$0D#$0A;
_2:
 end;
UtRes('TID10R',Result);  { вывод результатов в файл TID10R.res }
exit;
end;

end.

Результат:

         X                   A                            E
       2.         -6.207939493-02      0.000000000+00
       2.25     -1.081979987-01      2.160001819-05
       2.50     -8.205740855-02      1.682211116-05
       2.75     -6.392797405-02     -8.498944868-06
       3.00     -4.978715625-02     -6.618985062-06
       3.25     -3.877427628-02     -5.154870644-06
       3.50     -3.019743673-02     -4.014617251-06
       3.75     -2.343970501-02      7.495577201-05
       4.00     -1.832136756-02      7.671346584-05
       4.25     -8.102809959-03      0.000000000+00