Текст подпрограммы и версий
de22r_p.zip , de22e_p.zip 		Тексты тестовых примеров
tde22r_p.zip , tde22e_p.zip

Подпрограмма:  DE22R (модуль DE22R_p)

Назначение

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

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

Решается задача Коши для системы M обыкновенных дифференциальных уравнений первого порядка 

          Y ' = F (X, Y) ,

          Y = ( y1, ..., yM ) ,
          F = ( f1 (X, y1,..., yM),..., fM (X, y1,..., yM) )

 с начальными условиями, заданными в точке XN:
          Y(XN) = YN ,     YN = ( y10,..., yM0 ) ,

методом рациональной экстраполяции Грэгга - Булирша - Штера. Решение вычисляется в одной точке XK, которая является концом интервала интегрирования. Bсе компоненты решения могут проверяться на точность по относительной погрешности, абсолютной погрешности и стандартной погрешности. Тип погрешности специфицируется пользователем при обращении к подпрограмме.

Пусть Uin и Vin являются двумя последовательными приближениями к значению точного решения Yin в узле  xn. Тогда Vin принимается за искомое значение решения в узле  xn, если выполнен один из критериев сходимости:

1.  | Uin - Vin | ≤ EPS * | Zin | , где Zin = max{ Yi1,Yi2,...,Yin - 1} для  i = 1,...,M (стандартная погрешность)
2.  | Uin - Vin | ≤ EPS * | Vin | , для  i = 1,...,M (относительная погрешность)
3.  | Uin - Vin | ≤ EPS  , для  i = 1,...,M (абсолютная погрешность).

Bulirsch R., Stoer J., Numerical treatment of ordinary differential equations by extrapolation methods. Numerische Mathematik, 8(1) 1966, 1-13.

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

procedure DE22R(F :Proc_F_DE; M :Integer; XN :Real;
                var YN :Array of Real; XK :Real; HMIN :Real;
                EPS :Real; IORDER :Integer; IU :Integer;
                var H :Real; var Y :Array of Real;
                var YPM :Array of Real; var DELTY :Array of Real;
                var RAB :Array of Real; var IERR :Integer);


   Параметры
F - имя подпрограммы вычисления значений правой части дифференциального уравнения. Первый оператоp подпрограммы должен иметь вид:
procedure F (X :Real; var Y :Array of Real; var DY :Array of Real; M :Integer);
Здесь: X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY. B случае системы уравнений, т.е. когда M ≠ 1, параметры Y и DY представляют массивы длины M (тип параметров X, Y и DY: вещественный);
M - количество уравнений в системе (тип: целый);
XN, YN - начальные значения аргумента и решения. B случае системы уравнений (т.е. M ≠ 1) YN представляет массив длины M (тип: вещественный);
XK - значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования). XK может быть больше, меньше или pавно XN (тип: вещественный);
HMIN - минимальное значение абсолютной величины шага, который разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
EPS - допустимая погрешность, с которой требуется вычислить все компоненты решения; тип погрешности специфицируется с помощью параметpа IU (тип: вещественный);
IORDER - целая переменная, указывающая максимальный допустимый порядок метода рацональной экстраполяции; IORDER должен быть меньше 7;
IU - целый указатель типа погрешности численного решения:
IU = 1 - для стандартной погрешности;
IU = 2 - для относительной погрешности;
IU = 3 - для абсолютной погрешности;
H - вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если XK > XN, отрицательным, если XK < XN, или без такого учета в виде абсолютной величины;
Y - искомое решение задачи Коши, вычисленное подпрограммой для значения аргумента XK. Для системы уравнений (когда M ≠ 1) Y задается массивом длины M. B случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: вещественный);
         YPM -
       DELTY   		одномерные вещественные рабочие массивы длины M;
RAB - одномерный вещественный рабочий массив длины 29*M;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом:
IERR= 1 - когда неправильно задан параметр IORDER, т.е.IORDER < 1 или IORDER > 6 максимальный допустимый порядок метода); в этом случае интегрирование данной системы уравнений ведется со значением IORDER=6;
IERR=65 - когда какая - нибудь компонента pешения не может быть вычислена с требуемой точностью EPS; в этом случае интегрирование прекращается. При желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров H, HMIN и IORDER. 
   Версии
DE22E - вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом рациональной экстраполяции Грэгга - Булирша - Штера с расширенной (Extended) точностью. При этом параметры XN, YN, XK, HMIN, EPS, H, Y, YPM, DELTY, RAB и параметры X, Y и DY в подпрограмме F должны иметь тип Extended. 
   Вызываемые подпрограммы
DE20R - выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений первого порядка методом рациональной экстраполяции Грэгга - Булирша - Штера; вызывается при работе подпрограммы DE22R.
DE20E - выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений первого порядка методом рацональной экстраполяции Грэгга - Булирша - Штера с расширенной (Extended) точностью; вызывается при работе подпрограммы DE22E.
UTDE10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE22R.
UTDE11 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE22E. 
   Замечания по использованию
 

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

При работе подпрограммы значения параметров M, XN, YN, XK, HMIN, EPS, IORDER, IU сохраняются. Если после работы подпрограммы нет необходимости иметь начальное значение решения YN, то параметры YN и Y при обращении к ней можно совместить.

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

Использование подпрограммы иллюстрируется на примере: 

           y1'  =  0.2 ( y4 - y1 )
           y2'  =  y1  +  2 ( y2 - y2 y3 )
           y3'  =  y4  -  ( y3 - y2 y3 )
           y4'  =  10 y1 -  ( 61 - 0.13 x ) y4  +  0.13 x ,    0 ≤ x ≤ 8

          y1(0)  =  y2(0)  =  y3(0)  =  y4(0)  =  0 .

Приводятся подпрограмма вычисления значений правой части и фрагмент вызывающей программы, а также результаты счета. 

Unit TDE22R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE22R_p, DE22R_p;

function TDE22R: String;

implementation

function TDE22R: String;
var
M,IORDER,IU,_i,IERR :Integer;
XN,XK,H,HMIN,EPS :Real;
RАВ :Array [0..115] of Real;
Y :Array [0..3] of Real;
DELTY :Array [0..3] of Real;
YРМ :Array [0..3] of Real;
begin
Result := '';  { результат функции }
M := 4;
Y[0] := 0.0;
Y[1] := 0.0;
Y[2] := 0.0;
Y[3] := 0.0;
XN := 0.0;
ХК := 8.0;
H := 0.01;
HMIN := 1E-6;
EPS := 1.E-5;

IORDER := 4;
IU := 2;
DE22R(FDE22R,M,XN,Y,XK,HMIN,EPS,IORDER,IU,
     H,Y,YPM,DELTY,RAB,IERR);
Result := Result + Format(' %20.16f',[XN]) + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 3 do
 begin
  Result := Result + Format('%20.16f ',[Y[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TDE22R',Result);  { вывод результатов в файл TDE22R.res }
exit;
end;

end.

Unit FDE22R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;

procedure FDE22R(T :Real; var Y :Array of Real; var Z :Array of Real;
                M :Integer);
implementation

procedure FDE22R(T :Real; var Y :Array of Real; var Z :Array of Real;
                M :Integer);
var
R23,CT :Real;
begin
R23 := Y[2-1]*Y[3-1];
CT := 0.13*T;
Z[1-1] := 0.2*(Y[4-1]-Y[1-1]);
Z[2-1] := Y[1-1]+2.0*(Y[2-1]-R23);
Z[3-1] := Y[4-1]-(Y[3-1]-R23);
Z[4-1] := 10.0*Y[1-1]-(61.0-CT)*Y[4-1]+CT;
end;

end.

Результаты:

          Y(1)  =  0.9238470846 * 10-2
          Y(2)  =  0.4820971549 * 10-2
          Y(3)  =  1.667112950
          Y(4)  =  0.188434937 * 10-1

          IERR  =  0