БЧА НИВЦ МГУ. DE95R_P (DE91R, DE93R, DE97R). Задача Коши для жестких систем и систем с большой константой Липшица

Текст подпрограммы и версий
de95r_p.zip de95e_p.zip de91r_p.zip de91e_p.zip de93r_p.zip de93e_p.zip de97r_p.zip de97e_p.zip

Тексты тестовых примеров
tde95r_p.zip tde95e_p.zip tde91r_p.zip tde91e_p.zip tde93r_p.zip tde93e_p.zip tde97r_p.zip tde97e_p.zip

Подпрограмма: DE95R (модуль DE95R_p) (версии: DE91R, DE93R, DE97R)

Назначение

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

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

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

                             Y '   =   F(X,Y)  , 
       Y   =   ( y₁,...,y_M )  ,   F   =   ( f₁( X, y₁,..., y_M ),..., f_M( X, y₁,..., y_M ) )

методом типа Розенброка четвертого порядка.

Решение вычисляется в одной точке XK, которая является концом интервала интегрирования. Для всех компонент решения осуществляется контроль точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P (называемой границей перехода), то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.

Артемьев С.С., Демидов Г.В. A - устойчивый метод типа Розенброка четвертого порядка точности решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений. Некоторые проблемы вычислительной и прикладной математики. Новосибирск: "Наука", 1975.

Современные численные методы решения обыкновенных дифференциальных уравнений. Под ред. Дж.Холла и Дж.Уатта. М.: "Мир", 1979.

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

procedure DE95R(F :Proc_F_DE; FJ :Proc_F_DE; FX :Proc_F_DE;
                M :Integer; XN :Real; var YN :Array of Real;
                XK :Real; HMIN :Real; EPS :Real; P :Real;
                var H :Real; var Y :Array of Real;
                var R :Array of Real; var IERR :Integer);

   Параметры

F -	имя подпрограммы вычисления значений правой части дифференциальных уравнений. Первый оператор подпрограммы должен иметь вид: procedure F (X :Real; var Y :Array of Real; var DY :Array of Real; M :Integer);: Здесь X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY; в случае системы уравнений, т.е. когда M ≠ 1, параметры Y и DY представляют одномерные массивы длиной M (тип параметров X, Y и DY: вещественный);
FJ -	имя подпрограммы вычисления значений элементов матрицы Якоби ∂f / ∂y правой части системы. Первый оператор подпрограммы должен иметь вид: procedure FJ (X :Real; var Y :Array of Real; var Z :Array of Real; M :Integer); Здесь X, Y - значения независимой и зависимой переменных, соответственно. В случае системы уравнений, т.е. когда M ≠ 1, параметр Y представляет собой одномерный массив длины M, а параметр Z - двумерный массив размера M*M. Значения элементов матрицы Якоби ∂f / ∂y правой части системы должны быть помещены в массив Z, при этом частная производная от правой части I - го уравнения по J - ой переменной Y (J) запоминается в элементе Z (I,J) (тип параметров X, Y и Z: вещественный);
FX -	имя подпрограммы вычисления значений частных производных по X правой части системы ∂f / ∂x. Первый оператор подпрограммы должен иметь вид: procedure FX (X :Real; var Y :Array of Real; var Z1 :Array of Real; M :Integer); Здесь: X, Y - значения независимой и зависимой переменных, соответственно. В случае системы уравнений, т.е. когда M ≠ 1, параметры Y и Z1 представляют собой одномерные массивы длиной M. Значения частных производных по X правой части системы ∂f / ∂x должны быть помещены в массив Z1; при этом частная производная от правой части I - го уравнения запоминается в элементе Z1 (I) (тип параметров X, Y и Z1: вещественный);
M -	количество уравнений в системе (тип: целый);
XN, YN -	начальные значения аргумента и решения; в случае системы уравнений (т.е. M ≠ 1) YN представляет одномерный массив длины M (тип: вещественный);
XK -	значение аргумента, при котором требуется вычислить решение задачи Коши (конец интервала интегрирования); XK может быть больше, меньше или равно XN (тип: вещественный);
HMIN -	минимальное значение абсолютной величины шага, которое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
EPS -	допустимая мера погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный);
P -	граница перехода, используемая при оценке меры погрешности решения (тип: вещественный);
H -	вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если XN < XK, отрицательным, если XN > XK, или без всякого учета в виде абсолютной величины;
Y -	искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента XK; для системы уравнений (когда M ≠ 1) задается одномерным массивом длины M. В случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: вещественный);
R -	одномерный рабочий массив вещественного типа длины 3MM + 11*M + 1;
IERR -	целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS; в этом случае интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров H и HMIN.

   Версии

DE97R -	вычисление решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом типа Розенброка четвертого порядка. Отличие подпрограммы DE97R от подпрограммы DE95R состоит в том, что при обращении к подпрограмме DE97R не требуется задавать подпрограмму вычисления матрицы Якоби системы и подпрограмму вычисления частных производных правой части системы по x. Все частные производные от правой части вычисляются в подпрограмме DE97R с помощью разностных отношений. Первый оператор подпрограммы DE97R имеет вид: procedure DE97R(F :Proc_F_DE; M :Integer; XN :Real; var YN :Array of Real; XK :Real; HMIN :Real; EPS :Real; P :Real; var H :Real; var Y :Array of Real; var R :Array of Real; var IERR :Integer); Список формальных параметров подпрограммы DE97R отличается от списка параметров подпрограммы DE95R отсутствием параметров FJ и FX. Параметры подпрограммы DE97R имеют тот же смысл, что и одноименные параметры подпрограммы DE95R.
DE91R -	вычисление решения задачи Коши для автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом типа Розенброка четвертого порядка. Данная подпрограмма предназначена для интегрирования систем уравнений, в правые части которых не входит независимая переменная x, т.е. систем вида y ' = f (y). Первый оператор подпрограммы DE91R имеет вид: procedure DE91R(F :Proc_F_DE; FJ :Proc_F_DE; M :Integer; XN :Real; var YN :Array of Real; XK :Real; HMIN :Real; EPS :Real; P :Real; var H :Real; var Y :Array of Real; var R :Array of Real; var IERR :Integer); Список формальных параметров подпрограммы DE91R отличается от списка параметров DE95R отсутствием параметра FX. Параметры подпрограммы DE91R имеют тот же смысл, что и одноименные параметры подпрограммы DE95R, кроме параметра R. В подпрограмме DE91R параметр R представляет одномерный вещественный рабочий массив длины 3M² + 8M + 1 .
DE93R -	вычисление решения задачи Коши для автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Розенброка четвертого порядка. Данная подпрограмма предназначена для интегрирования систем уравнений, в правые части которых не входит независимая переменная x, т.е. систем вида y ' = f (y). Отличие подпрограммы DE93R от подпрограммы DE95R состоит в том, что при обращении к подпрограмме DE93R не требуется задавать подпрограмму вычисления матрицы Якоби системы ∂f / ∂y и подпрограмму вычисления частных производных ∂f / ∂x. Все частные производные от правой части вычисляются в подпрограмме DE93R с помощью разностных отношений. Первый оператор подпрограммы DE93R имеет вид: procedure DE93R(F :Proc_F_DE; M :Integer; XN :Real; var YN :Array of Real; XK :Real; var HMIN :Real; var EPS :Real; var P :Real; var H :Real; var Y :Array of Real; var R :Array of Real; var IERR :Integer); Список формальных параметров подпрограммы DE93R отличается от списка параметров подпрограммы DE95R отсутствием параметров FJ и FX. Параметры подпрограммы DE93R имеют тот же смысл, что и одноименные параметры подпрограммы DE95R, кроме параметра R. В подпрограмме DE93R параметр R представляет одномерный вещественный рабочий массив длины 3M² + 8M + 1 .
DE95E -	вычисление решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом типа Розенброка четвертого порядка с расширенной (Extended) точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE95R; при этом параметры XN, YN, XK, HMIN, EPS, P, H, Y, R и параметры X, Y, DY, Z, Z1 в подпрограммах F, FJ и FX должны иметь тип Extended.
DE97E -	вычисление решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом типа Розенброка четвертого порядка с расширенной (Extended) точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE97R; при этом параметры XN, YN, XK, HMIN, EPS, P, H, Y, R и параметры X, Y, DY в подпрограмме F должны иметь тип Extended.
DE91E -	вычисление решения задачи Коши для автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом типа Розенброка четвертого порядка с расширенной (Extended) точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE91R; при этом параметры XN, YN, XK, HMIN, EPS, P, H, Y, R и параметры X, Y, DY, Z в подпрограммах F и FJ должны иметь тип Extended.
DE93E -	вычисление решения задачи Коши для автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интергирования методом типа Розенброка четвертого порядка с расширенной (Extended) точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE93R; при этом параметры XN, YN, XK, HMIN, EPS, P, H, Y, R и параметры X, Y, DY в подпрограмме F должны иметь тип Extended.

   Вызываемые подпрограммы

DE94R - DE94E	выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка с обыкновенной и расширенной (Extended) точностью. Вызываются при работе подпрограмм DE95R и DE95E соответственно.
DE96R - DE96E	выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка с обыкновенной и расширенной (Extended) точностью. Вызываются при работе подпрограмм DE97R и DE97E соответственно.
DE90R - DE90E	выполнение одного шага численного интегрирования автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка с обыкновенной и расширенной (Extended) точностью. Вызываются при работе подпрограмм DE91R и DE91E соответственно.
DE92R - DE92E	выполнение одного шага численного интегрирования автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка с обыкновенной и расширенной (Extended) точностью. Вызываются при работе подпрограмм DE93R и DE93E соответственно.
UTDE20 - UTDE21	подпрограммы выдачи диагностических сообщений. Подпрограмма UTDE20 вызывается при работе подпрограмм DE91R, DE93R, DE95R, DE97R; подпрограмма UTDE21 вызывается при работе подпрограмм DE91E, DE93E, DE95E, DE97E.

   Замечания по использованию

В общем случае заданная точность не гарантируется.

При работе подпрограммы и ее версий значения параметров M, XN, YN, XK, HMIN, EPS, P сохраняются. При работе подпрограмм F, FJ и FX значения параметров M, X, Y не должны изменяться.

Если после работы подпрограммы нет необходимости иметь начальное значение решения YN, то параметры YN и Y при обращении к ней можно совместить. При этом следует иметь в виду, что в случае аварийного выхода из подпрограммы, т.е. со значением IERR = 65, значение параметра YN будет испорчено.

В подпрограмме DE95R и ее версиях используется одношаговый метод типа Розенброка четвертого порядка точности. На каждом шаге h интегрирования, выполняемом из текущего узла интегрирования x_n, вычисляются четыре значения правой части системы дифференциальных уравнений. Среди этих значений есть одно, которое должно вычисляться при значении независимой переменной x, равном x_n - h. Это значение аргумента находится слева от точки x_n при h > 0 и справа от точки x_n при h < 0.

В частности, при выполнении первого шага h из начальной точки x_N (напомним, что абсолютная величина первого шаг h берется равной абсолютной величине значения, заданного параметром H при обращении к подпрограмме, а знак первого шага определяется соотношением значений параметров XN и XK) значение x_N - h не будет принадлежать интервалу интегрирования, ограниченному значениями, заданными параметрами XN и XK при обращении к подпрограмме DE95R (или ее версиями). Это следует учитывать при составлении подпрограммы F вычисления правой части системы. Если правая часть системы не определена для x, не принадлежащих интервалу интегрирования, то попытка вычислить правую часть для указанного значения аргумента x может привести к аварийному прерыванию.

Если при решении системы диффиренциальных уравнений не задаются подпрограммы FJ вычисления матрицы Якоби и FX вычисления частных производных по x правой части системы (т.е. производится обращение к подпрограммам DE93R, DE97R, DE93E, DE97E), то все частные производные по y и по x в этих подпрограммах аппроксимируются центральными разностными отношениями. Замена точных значений производных разностными аппроксимациями может привести к росту погрешности приближенного решения системы дифференциальных уравнений по сравнению с тем, что было бы, если бы все частные производные вычислялись точно с помощью подпрограмм FJ и FX. Даже если будет мала погрешность аппроксимации частных производных (например, если правая часть системы является линейной функцией своих аргументов, то погрешность аппроксимации частных производных равна нулю), может оказаться значительной вычислительная погрешность, особенно, когда вычисления выполняются с одинарной точностью. При этом вычислительная погрешность может даже превосходить верхний предел погрешности приближенного решения, заданный при обращении к подпрограмме параметром EPS. В этом случае целесообразно использовать не подпрограммы DE93R, DE97R, а их версии, выполняющие вычисления с удвоенным числом значащих цифр, т.е. подпрограммы DE93E, DE97E.

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

      
   Пример 1.
   Решается задача Коши

      y₁'  =  -100 y₁ 
      y₂'  =  -100 y₁  - 2 y₂  + 20 e^{- 100 x} + 2 e^{- x} cos x 
      y₃'  =  -100 y₁  + 9998 y₂ - 9990 y₃ - 10 y₄  + 20 e^{- 100 x} + 2 e^{- x} cos x 
      y₄'  =  -100 y₁  + 9988 y₂  + 20 y₃ - 10010 y₄  + 20 e^{- 100 x} + 2 e^{- x} cos x  

      0 ≤ x ≤ 10 ,  y₁(0) = 10 ,  y₂(0) = 11 ,  y₃(0) = 111 ,  y₄(0) = 11   

   Точное решение задачи имеет вид:

      y₁  =  10 e^{-100 x} 
      y₂  =  10 e^{- 100 x} + e^{- x} cos x + e^{- x} sin x 
      y₃  =  y₂  + 100 e^{- 10000 x} cos 10x 
      y₄  =  y₃  + 100 e^{- 10000 x} sin 10x 

   Матрица Якоби правой части системы имеет вид:

        |   -100       0            0             0      
        |   -100      -2            0             0      
        |   -100       9998     -9990     -10     
        |   -100       9988      20          -10010   

   Частные производные по  x правой части системы имеют вид:

        |                            0                 
        |   -2000 e^{- 100 x} - 2 e^{- x} (cos x + sin x)
        |   -2000 e^{- 100 x} - 2 e^{- x} (cos x + sin x)
        |   -2000 e^{- 100 x} - 2 e^{- x} (cos x + sin x)

Ниже приводятся фрагмент вызывающей программы для DE95R, подпрограммы F, FJ, FX и результаты счета, полученные после нескольких обращений к подпрограмме DE95R.

Unit TDE95R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE95R_p, FJDE95R_p, FXDE95R_p, DE95R_p;

function TDE95R: String;

implementation

function TDE95R: String;
var
M,IH,IERR :Integer;
XN,HMIN,EPS,P,XK,Y1,Y2,Y3,Y4,H :Real;
YN :Array [0..3] of Real;
Y :Array [0..3] of Real;
R :Array [0..92] of Real;
label
_10;
begin
Result := '';  { результат функции }
M := 4;
XN := 0.0;
YN[0] := 10.0;
YN[1] := 11.0;
YN[2] := 111.0;
YN[3] := 111.0;
HMIN := 1.E-10;
EPS := 1.E-2;
P := 1.E3;
ХК := 10.0;
{   ВЫЧИСЛЕНИЕ ТОЧНОГО РЕШЕНИЯ CИCTEMЫ: }
Y1 := 10.0*Exp(-100.0*XK);
Y2 := Y1+Exp(-XK)*(Cos(XK)+Sin(XK));
Y3 := Y2+100.0*Exp(-1.E4*XK)*Cos(10.0*XK);
Y4 := Y3+100.0*Exp(-1.E4*XK)*Sin(10.0*XK);
Result := Result + Format(' %20.16f  %20.16f  %20.16f  %20.16f ',
 [Y1,Y2,Y3,Y4]) + #$0D#$0A;
IH := 0;
_10:
IH := IH + 1;
H := 0.01;
DE95R (FDE95R,FJDE95R,FXDE95R,M,XN,YN,XK,HMIN,EPS,P,H,Y,R,IERR);
Result := Result + Format('%s',[' ']) + #$0D#$0A; 
Result := Result + Format('%s',[' EPS= ']);
Result := Result + Format('%20.16f ',[EPS]) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A; 
Result := Result + Format(' %20.16f  %20.16f  %20.16f  %20.16f ',
 [Y[0],Y[1],Y[2],Y[3]]) + #$0D#$0A;
Result := Result + Format('%s',[' H= ']);
Result := Result + Format('%20.16f ',[H]) + #$0D#$0A;
EPS := EPS*1.E-2;
if ( IH < 3 ) 
 then goto _10;
UtRes('TDE95R',Result);  { вывод результатов в файл TDE95R.res }
exit;
end;

end.

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

procedure fde95r(X :Real; var Y :Array of Real; var DY :Array of Real;
                M :Integer);

implementation

procedure fde95r(X :Real; var Y :Array of Real; var DY :Array of Real;
                M :Integer);
var
T1 :Real;
begin
DY[0] := -100.0*Y[0];
T1 := 20.0*Exp(-100.0*X)+2.0*Exp(-X)*Cos(X);
DY[1] := DY[0]-2.0*Y[1]+T1;
DY[2] := DY[0]+9998.0*Y[1]-9990.0*Y[2]-10.0*Y[3]+T1;
DY[3] := DY[0]+9988.0*Y[1]+20.0*Y[2]-10010.0*Y[3]+T1;
end;

end.

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

procedure fjde95r(X :Real; var Y :Array of Real; var DF :Array of Real;
                M :Integer);

implementation

procedure fjde95r(X :Real; var Y :Array of Real; var DF :Array of Real;
                M :Integer);
begin
DF[0] := -100.0;
DF[4] := 0.0;
DF[8] := 0.0;
DF[12] := 0.0;
DF[1] := -100.0;
DF[5] := -2.0;
DF[9] := 0.0;
DF[13] := 0.0;
DF[2] := -100.0;
DF[6] := 9998.0;
DF[10] := -9990.0;
DF[14] := -10.0;
DF[3] := -100.0;
DF[7] := 9988.0;
DF[11] := 20.0;
DF[15] := -10010.0;
end;

end.

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

procedure fxde95r(X :Real; var Y :Array of Real; var DX :Array of Real;
                M :Integer);

implementation

procedure fxde95r(X :Real; var Y :Array of Real; var DX :Array of Real;
                M :Integer);
begin
DX[0] := 0.0;
DX[1] := -2.E3*Exp(-1.E2*X)-2.0*Exp(-X)*(Cos(X)+Sin(X));
DX[2] := DX[1];
DX[3] := DX[1];
end;

end.

 Результаты: 

              Y1                                   Y2 
      0.000000000000+00    -6.279230870976-05 
              Y3                                   Y4 
     -6.279230870976-05    -6.279230870976-05 

      после первого обращения к подпрограмме -
      EPS = 1.0-02 
              Y(1)                                 Y(2) 
      7.146873880164-13    -6.764660892966-05 
              Y(3)                                 Y(4) 
     -6.764660892122-05    -6.764660892966-05 
      H = 2.560000000001+00 
   
      после второго обращения к подпрограмме -
      EPS = 1.0-04 
              Y(1)                                 Y(2) 
      6.621227792960-20    -6.330159900469-05 
              Y(3)                                 Y(4) 
     -6.330159900392-05    -6.330159900414-05 
      H = 2.560000000001+00 

      после третьего обращения к подпрограмме -
      EPS = 1.0-06
              Y(1)                                 Y(2) 
      0.000000000000+00    -6.286382905407-05 
              Y(3)                                 Y(4) 
     -6.286382905418-05    -6.286382905407-05 
      H = 1.280000000001+00 
   
      после четвертого обращения к подпрограмме -
      EPS = 1.0-08 
              Y(1)                                 Y(2) 
      0.000000000000+00    -6.279451115976-05 
              Y(3)                                 Y(4) 
     -6.279451115976-05    -6.279451116020-05 
      H = 3.200000000002-01 

   Пример 2.

Решается система из примера 1. Приводятся фрагмент вызывающей программы для DE97R, подпрограмма F и результаты счета, полученные после нескольких обращений к подпрограмме DE97R.

   
Unit tde97r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE97R_p, DE97R_p;

function tde97r: String;

implementation

function tde97r: String;
var
M,IH,IERR :Integer;
XN,HMIN,EPS,P,XK,Y1,Y2,Y3,Y4,H :Real;
YN :Array [0..3] of Real;
Y :Array [0..3] of Real;
R :Array [0..92] of Real;
label
_10;
begin
Result := '';  { результат функции }
M := 4;
XN := 0.0;
YN[0] := 10.0;
YN[1] := 11.0;
YN[2] := 111.0;
YN[3] := 111.0;
HMIN := 1.E-10;
EPS := 1.E-2;
P := 1.E3;
XK := 10.0;
{   BЫЧИCЛEHИE TOЧHOГO PEШEHИЯ CИCTEMЫ: }
Y1 := 10.0*Exp(-100.0*XK);
Y2 := Y1+Exp(-XK)*(Cos(XK)+Sin(XK));
Y3 := Y2+100.0*Exp(-1.E4*XK)*Cos(10.0*XK);
Y4 := Y3+100.0*Exp(-1.E4*XK)*Sin(10.0*XK);
Result := Result + Format(' %20.16f  %20.16f  %20.16f  %20.16f ',
 [Y1,Y2,Y3,Y4]) + #$0D#$0A;
IH := 0;
_10:
IH := IH + 1;
H := 0.01;
DE97R (FDE97R,M,XN,YN,XK,HMIN,EPS,P,H,Y,R,IERR);
Result := Result + Format('%s',[' ']) + #$0D#$0A; 
Result := Result + Format('%s',[' ']) + #$0D#$0A; 
Result := Result + Format('%s',[' EPS= ']);
Result := Result + Format('%20.16f ',[EPS]) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A; 
Result := Result + Format(' %20.16f  %20.16f  %20.16f  %20.16f ',
 [Y[0],Y[1],Y[2],Y[3]]) + #$0D#$0A;
Result := Result + Format('%s',[' H= ']);
Result := Result + Format('%20.16f ',[H]) + #$0D#$0A;
EPS := EPS*1.E-1;
if ( IH < 3 ) 
 then goto _10;
UtRes('tde97r',Result);  { вывод результатов в файл tde97r.res }
exit;
end;

end.

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

procedure fde97r(X :Real; var Y :Array of Real; var DY :Array of Real;
                M :Integer);

implementation

procedure fde97r(X :Real; var Y :Array of Real; var DY :Array of Real;
                M :Integer);
var
T1 :Real;
begin
DY[0] := -100.0*Y[0];
T1 := 20.0*Exp(-100.0*X)+2.0*Exp(-X)*Cos(X);
DY[1] := DY[0]-2.0*Y[1]+T1;
DY[2] := DY[0]+9998.0*Y[1]-9990.0*Y[2]-10.0*Y[3]+T1;
DY[3] := DY[0]+9988.0*Y[1]+20.0*Y[2]-10010.0*Y[3]+T1;
end;

end.

      Результаты: 
              Y1                                   Y2 
      0.000000000000+00    -6.279230870976-05 
              Y3                                   Y4 
     -6.279230870976-05    -6.279230870976-05 

      после первого обращения к подпрограмме -
  
      EPS = 1.0-02 
              Y(1)                                 Y(2) 
      7.146874056701-13    -6.764660658287-05 
              Y(3)                                 Y(4) 
     -6.764662603498-05    -6.764660998348-05 
      H = 2.560000000001+00 

      после второго обращения к подпрограмме -
      EPS = 1.0-04 
              Y(1)                                 Y(2) 
      6.614556599143-02    -6.330157681200-05 
              Y(3)                                 Y(4) 
     -6.330157682632-05    -6.330157680712-05 
      H = 2.560000000001+00 
   
      после третьего обращения к подпрограмме -
      EPS = 1.0-06 
              Y(1)                                 Y(2) 
      2.733180996095-20    -6.286373127617-05 
              Y(3)                                 Y(4) 
     -6.286373127695-05    -6.286373127762-05 
      H = 1.280000000001+00

   Пример 3.
   Решается задача Коши
 
      y₁'  =  - 10000 y₁  + 100 y₂ - 10 y₃  + y₄ 
      y₂'  =  - 1000 y₂  + 10 y₃ -10 y₄ 
      y₃'  =  - y₃  + 10 y₄ 
      y₄'  =  - 0.1 y₄ 

      0 ≤ x ≤ 20 ,  y₁(0) = 1 ,  y₂(0) = 1 ,  y₃(0) = 1 ,  y₄(0) = 1 

   Точное решение задачи имеет вид:

      y₁  = (1 - d - f - g) t₃  + d t₂  + f t₁  + g y₄ 
      y₂  = (1 - a - b) t₂  + a t₁  + b y₄ 
      y₃  = (- 9.1 / 0.9) t₁ + (10 / 0.9) y₄
      y₄  = e^{- 0.1 x} 
 
 где   t₁ = e^{- x}  ,    t₂ = e^{- 1000 x}  ,     t₃ = e^{- 10000 x}   ,

  a = - 91 / (0.9*999) ,    b = 91 / (0.9*999.9)  ,    d =  ( 1 - a - b ) / 90  ,
              
  f = [ 100 a + (91 / 0.9) ] / 9999  ,    g = [ 100 b - (100 / 0.9)  + 1 ] / 9999.9  . 

   Матрица Якоби правой части системы имеет вид:

       -10000     100     -10      1     
         0           -1000    10    -10    
         0             0         -1      10    
         0             0          0      -0.1

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

  
Unit tde91r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE91R_p, FJDE91R_p, DE91R_p;

function tde91r: String;

implementation

function tde91r: String;
var
M,IH,IERR :Integer;
XN,HMIN,EPS,P,XK,Y4,T1,Y3,T2,A,B,C,Y2,T3,DD,FF,GG,QQ,Y1,H :Real;
YN :Array [0..3] of Real;
Y :Array [0..3] of Real;
R :Array [0..80] of Real;
label
_10;
begin
Result := '';  { результат функции }
M := 4;
XN := 0.0;
YN[0] := 1.0;
YN[1] := 1.0;
YN[2] := 1.0;
YN[3] := 1.0;
HMIN := 1.E-10;
EPS := 1.E-2;
P := 100.0;
XK := 20.0;
{  BЫЧИCЛEHИE TOЧHOГO PEШEHИЯ CИCTEMЫ: }
Y4 := Exp(-0.1*XK);
T1 := Exp(-XK);
Y3 := -(9.1/0.9)*T1 + (10.0/0.9)*Y4;
T2 := Exp(-1.E3*XK);
A := -91.0/(0.9*999.0);
B := 91.0/(0.9*999.9);
C := 1.0-A-B;
Y2 := C*T2+A*T1+B*Y4;
T3 := Exp(-1.E4*XK);
DD := C/90.0;
FF := (100.0*A+91.0/0.9)/9999.0;
GG := (100.0*B-100.0/0.9+1.0)/9999.9;
QQ := 1.0-DD-FF-GG;
Y1 := QQ*T3+DD*T2+FF*T1+GG*Y4;
Result := Result + Format(' %20.16f  %20.16f  %20.16f  %20.16f ',
 [Y1,Y2,Y3,Y4]) + #$0D#$0A;
IH := 0;
_10:
IH := IH + 1;
H := 0.01;
DE91R (FDE91R,FJDE91R,M,XN,YN,XK,HMIN,EPS,P,H,Y,R,IERR);
Result := Result + Format('%s',[' ']) + #$0D#$0A; 
Result := Result + Format('%s',[' ']) + #$0D#$0A; 
Result := Result + Format('%s',[' ']) + #$0D#$0A; 
Result := Result + Format('%s',[' EPS= ']);
Result := Result + Format('%20.16f ',[EPS]) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A; 
Result := Result + Format(' %20.16f  %20.16f  %20.16f  %20.16f ',
 [Y[0],Y[1],Y[2],Y[3]]) + #$0D#$0A;
Result := Result + Format('%s',[' H= ']);
Result := Result + Format('%20.16f ',[H]) + #$0D#$0A;
EPS := EPS*1.E-1;
if ( IH < 3 ) 
 then goto _10;
UtRes('tde91r',Result);  { вывод результатов в файл tde91r.res }
exit;
end;

end.

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

procedure fde91r(X :Real; var Y :Array of Real; var DY :Array of Real;
                M :Integer);

implementation

procedure fde91r(X :Real; var Y :Array of Real; var DY :Array of Real;
                M :Integer);
begin
DY[0] := -1.E4*Y[0]+1.E2*Y[1]-10.0*Y[2]+Y[3];
DY[1] := -1.E3*Y[1]+10.0*Y[2]-10.0*Y[3];
DY[2] := -Y[2]+10.0*Y[3];
DY[3] := -0.1*Y[3];
end;

end.

      Результаты: 

             Y1                                    Y2 
     -1.353352661869-03     1.368526917891-02 
             Y3                                    Y4 
      1.503725348452+00     1.353352832364-01 

      после первого обращения к подпрограмме -
      EPS = 1.0-02 
             Y(1)                                  Y(2) 
     -1.352902557240-03     1.368071768447-02 
             Y(3)                                  Y(4) 
      1.503225232165+00     1.352902708961-01 
      H = 1.024000000001+01 

      после второго обращения к подпрограмме -
      EPS = 1.0-04 
             Y(1)                                  Y(2) 
     -1.353311838306-03     1.368485637501-02 
             Y(3)                                  Y(4) 
      1.503679988920+00     1.353311999785-01 
      H = 2.560000000001+00 
   
      после третьего обращения к подпрограмме -
      EPS = 1.0-06 
             Y(1)                                  Y(2) 
     -1.353350127546-03     1.368524355264-02 
             Y(3)                                  Y(4) 
      1.503722532527+00     1.353350296827-01 
      H = 1.280000000001+00 

      после четвертого обращения к подпрограмме -
      EPS = 1.0-08 
             Y(1)                                  Y(2) 
     -1.353352646026-03     1.368526901869-02 
             Y(3)                                  Y(4) 
      1.503725330851+00     1.353352816498-01 
      H = 3.200000000002-01 

   Пример 4.

Решается задача Коши из примера 3. Приводятся фрагмент вызывающей программы для DE93R, подпрограмма F и результаты счета, полученные после нескольких обращений к подпрограмме DE93R.

  
Unit tde93r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE93R_p, DE93R_p;

function tde93r: String;

implementation

function tde93r: String;
var
M,IH,IERR :Integer;
XN,HMIN,EPS,P,XK,Y4,T1,Y3,T2,A,B,C,Y2,T3,DD,FF,GG,QQ,Y1,H :Real;
YN :Array [0..3] of Real;
Y :Array [0..3] of Real;
R :Array [0..80] of Real;
label
_10;
begin
Result := '';  { результат функции }
M := 4;
XN := 0.0;
YN[0] := 1.0;
YN[1] := 1.0;
YN[2] := 1.0;
YN[3] := 1.0;
HMIN := 1.E-10;
EPS := 1.E-2;
P := 100.0;
XK := 20.0;
{  BЫЧИCЛEHИE TOЧHOГO PEШEHИЯ CИCTEMЫ: }
Y4 := Exp(-0.1*XK);
T1 := Exp(-XK);
Y3 := -(9.1/0.9)*T1 + (10.0/0.9)*Y4;
T2 := Exp(-1.E3*XK);
A := -91.0/(0.9*999.0);
B := 91.0/(0.9*999.9);
C := 1.0-A-B;
Y2 := C*T2+A*T1+B*Y4;
T3 := Exp(-1.E4*XK);
DD := C/90.0;
FF := (100.0*A+91.0/0.9)/9999.0;
GG := (100.0*B-100.0/0.9+1.0)/9999.9;
QQ := 1.0-DD-FF-GG;
Y1 := QQ*T3+DD*T2+FF*T1+GG*Y4;
Result := Result + Format(' %20.16f  %20.16f  %20.16f  %20.16f ',
 [Y1,Y2,Y3,Y4]) + #$0D#$0A;
IH := 0;
_10:
IH := IH + 1;
H := 0.01;
DE93R (FDE93R,M,XN,YN,XK,HMIN,EPS,P,H,Y,R,IERR);
Result := Result + Format('%s',[' ']) + #$0D#$0A; 
Result := Result + Format('%s',[' ']) + #$0D#$0A; 
Result := Result + Format('%s',[' ']) + #$0D#$0A; 
Result := Result + Format('%s',[' EPS= ']);
Result := Result + Format('%20.16f ',[EPS]) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A; 
Result := Result + Format(' %20.16f  %20.16f  %20.16f  %20.16f ',
 [Y[0],Y[1],Y[2],Y[3]]) + #$0D#$0A;
Result := Result + Format('%s',[' H= ']);
Result := Result + Format('%20.16f ',[H]) + #$0D#$0A;
EPS := EPS*1.E-1;
if ( IH < 3 ) 
 then goto _10;
UtRes('tde93r',Result);  { вывод результатов в файл tde93r.res }
exit;
end;

end.

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

procedure fde93r(X :Real; var Y :Array of Real; var DY :Array of Real;
                M :Integer);

implementation

procedure fde93r(X :Real; var Y :Array of Real; var DY :Array of Real;
                M :Integer);
begin
DY[0] := -1.E4*Y[0]+1.E2*Y[1]-10.0*Y[2]+Y[3];
DY[1] := -1.E3*Y[1]+10.0*Y[2]-10.0*Y[3];
DY[2] := -Y[2]+10.0*Y[3];
DY[3] := -0.1*Y[3];
end;

end.

   Результаты: 

             Y1                                   Y2 
     -1.353352661869-03     1.368526917891-02 
             Y3                                   Y4 
      1.503725348452+00     1.353352832364-01 

   после первого обращения к подпрограмме -
      EPS = 1.0-02 
             Y(1)                                 Y(2) 
     -1.352902572352-03     1.368071785960-02 
             Y(3)                                 Y(4) 
      1.503225234967+00     1.352902709414-01 
      H = 1.024000000001+01 

   после второго обращения к подпрограмме -
      EPS = 1.0-04 
             Y(1)                                 Y(2) 
     -1.353311843749-03     1.368485643758-02 
             Y(3)                                 Y(4) 
      1.503679988458+00     1.353311999721-01 
      H = 2.560000000001+00 

   после третьего обращения к подпрограмме -
      EPS = 1.0-06 
             Y(1)                                 Y(2) 
     -1.353350112996-03     1.368524337161-02 
             Y(3)                                 Y(4) 
      1.503722531763+00     1.353350295963-01 
      H = 1.280000000001+00 

   после четвертого обращения к подпрограмме -
      EPS = 1.0-08 
             Y(1)                                 Y(2) 
     -1.353351314448-03     1.368525627954-02 
             Y(3)                                 Y(4) 
      1.503725330012+00     1.353352828808-01 
      H = 3.200000000002-01 


   Пример 5.

   Решается задача Коши

      y₁'  =  - a₂ y₁ y₂ 
      y₂'  =  2 a₁ y₅ x₉² - a₂ y₁ y₂ - a₅ x₆ y₂  + a₆ x₇ 
      y₃'  =  a₃ x₆ x₈ - a₇ y₃ 
      y₄'  =  a₆ x₇  + a₇ y₃  + a₄ x₆ 
      y₅'  =  - a₁ y₅ x₉²

где
      x₆  =  1 - a₈ - 2 y₁  + y₂ - y₃ - y₄  + 2y₅ 
      x₇  =  a₈ + y₁ - y₂ - 2 y₅ 
      x₈  =  a₉ - y₃ 
      x₉  =  - 0.8745 - a₈  + y₃  + y₄  + 2 y₅ 

      a₁  = 80,   a₂  = 29,   a₃  = 1,   a₄  = 0.288*10^-3  ,   a₅  = 166 ,  
      a₆  = 0.959 * 10^-4  ,   a₇  = 0.232 * 10^-3  ,   a₈  = 0.0477 ,   a₉  = 0.602 , 
            
      0 ≤ x ≤ 3000 .

      y₁(0)  =  1 ,   y₂(0)  =  0.0477 ,   y₃(0)  =  0 ,   y₄(0)  =  0 ,   y₅(0) = 0.5

Этот пример взят из статьи С.С.Артемьева и Г.В.Демидова, которая приведена в списке литературы, данном в разделе "Математическое описание".

Ниже приводятся фрагмент вызывающей программы для DE91R, подпрограммы F и FJ и результаты счета, полученные после двух обращений к подпрограмме DE91R со значениями параметра EPS, равными 10^-6 и 10^-7.

 
Unit tde91r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE91R_p, FJDE91R_p, DE91R_p;

function tde91r: String;

implementation

function tde91r: String;
var
M,IH,IERR :Integer;
XN,HMIN,EPS,P,XK,Y4,T1,Y3,T2,A,B,C,Y2,T3,DD,FF,GG,QQ,Y1,H :Real;
YN :Array [0..3] of Real;
Y :Array [0..3] of Real;
R :Array [0..80] of Real;
label
_10;
begin
Result := '';  { результат функции }
M := 4;
XN := 0.0;
YN[0] := 1.0;
YN[1] := 1.0;
YN[2] := 1.0;
YN[3] := 1.0;
HMIN := 1.E-10;
EPS := 1.E-2;
P := 100.0;
XK := 20.0;
{  BЫЧИCЛEHИE TOЧHOГO PEШEHИЯ CИCTEMЫ: }
Y4 := Exp(-0.1*XK);
T1 := Exp(-XK);
Y3 := -(9.1/0.9)*T1 + (10.0/0.9)*Y4;
T2 := Exp(-1.E3*XK);
A := -91.0/(0.9*999.0);
B := 91.0/(0.9*999.9);
C := 1.0-A-B;
Y2 := C*T2+A*T1+B*Y4;
T3 := Exp(-1.E4*XK);
DD := C/90.0;
FF := (100.0*A+91.0/0.9)/9999.0;
GG := (100.0*B-100.0/0.9+1.0)/9999.9;
QQ := 1.0-DD-FF-GG;
Y1 := QQ*T3+DD*T2+FF*T1+GG*Y4;
Result := Result + Format(' %20.16f  %20.16f  %20.16f  %20.16f ',
 [Y1,Y2,Y3,Y4]) + #$0D#$0A;
IH := 0;
_10:
IH := IH + 1;
H := 0.01;
DE91R (FDE91R,FJDE91R,M,XN,YN,XK,HMIN,EPS,P,H,Y,R,IERR);
Result := Result + Format('%s',[' ']) + #$0D#$0A; 
Result := Result + Format('%s',[' ']) + #$0D#$0A; 
Result := Result + Format('%s',[' ']) + #$0D#$0A; 
Result := Result + Format('%s',[' EPS= ']);
Result := Result + Format('%20.16f ',[EPS]) + #$0D#$0A;
Result := Result + Format('%s',[' ']) + #$0D#$0A; 
Result := Result + Format(' %20.16f  %20.16f  %20.16f  %20.16f ',
 [Y[0],Y[1],Y[2],Y[3]]) + #$0D#$0A;
Result := Result + Format('%s',[' H= ']);
Result := Result + Format('%20.16f ',[H]) + #$0D#$0A;
EPS := EPS*1.E-1;
if ( IH < 3 ) 
 then goto _10;
UtRes('tde91r',Result);  { вывод результатов в файл tde91r.res }
exit;
end;

end.

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

procedure fde91r(X :Real; var Y :Array of Real; var DY :Array of Real;
                M :Integer);

implementation

procedure fde91r(X :Real; var Y :Array of Real; var DY :Array of Real;
                M :Integer);
begin
DY[0] := -1.E4*Y[0]+1.E2*Y[1]-10.0*Y[2]+Y[3];
DY[1] := -1.E3*Y[1]+10.0*Y[2]-10.0*Y[3];
DY[2] := -Y[2]+10.0*Y[3];
DY[3] := -0.1*Y[3];
end;

end.

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

procedure fjde91r(X :Real; var Y :Array of Real; var DF :Array of Real;
                M :Integer);

implementation

procedure fjde91r(X :Real; var Y :Array of Real; var DF :Array of Real;
                M :Integer);
begin
DF[0] := -1.E4;
DF[4] := 1.E2;
DF[8] := -10.0;
DF[12] := 1.0;
DF[1] := 0.0;
DF[5] := -1.E3;
DF[9] := 10.0;
DF[13] := -10.0;
DF[2] := 0.0;
DF[6] := 0.0;
DF[10] := -1.0;
DF[14] := 10.0;
DF[3] := 0.0;
DF[7] := 0.0;
DF[11] := 0.0;
DF[15] := -0.1;
end;

end.
   
   Результаты: 

   после первого обращения к подпрограмме -
      EPS = 1.0-06 
             Y(1)                                 Y(2)                                Y(3) 
      7.743161591600-02     3.837857646993-05     5.035858808897-01 
             Y(4)                                 Y(5) 
      3.578708481186-01     3.240421442632-02 
      H = 1.638400000001+02 

   после второго обращения к подпрограмме -
      EPS = 1.0-07 
             Y(1)                                 Y(2)                                Y(3) 
      7.743302086658-02     3.837918184429-05     5.035845024640-01 
             Y(4)                                 Y(5) 
      3.578709420071-01     3.240497334838-02 
      H = 8.192000000004+01