Текст подпрограммы и версий
de09r_p.zip , de09e_p.zip 		Тексты тестовых примеров
tde09r_p.zip , tde09e_p.zip

Подпрограмма:  DE09R (модуль DE09R_p)

Назначение

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

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

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

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

с начальными условиями, заданными в точке XN: 

                         Y(XN)  =  YN ,   YN  =  (  y1 0,..., yM 0  ) ,

методом Инглэнда пятого порядка точности. Решение вычисляется в одной точке XK, которая является концом интервала интегрирования.

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

1. 

Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. - М.: Наука, 1987.

2.  England R. Error estimates for Runge - Kutta type solutions to systems of ordinary differential equations. The Computer Journal. 1969 - v.12, No.2.

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

procedure DE09R(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 R1 :Array of Real;
                var R2 :Array of Real; var R3 :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 представляют одномерные массивы длины М (тип параметров X, Y и DY: вещественный);
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 (тип: вещественный);
R1 - двумерный вещественный рабочий массив размера М*5;
R2, R3 - одномерные вещественные рабочие массивы длины М;
IERR - целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS. В этом случае интегрирование системы прекращается; при желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров HMIN и H. 
   Версии
DE09E - вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Инглэнда с удвоенным числом значащих цифр. При этом параметры XN, YN, XK, HMIN, EPS, P, H, Y, R1, R2, R3 и параметры X, Y и DY в подпрограмме F должны иметь тип Extended. 
   Вызываемые подпрограммы
UTDE20 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE09R.
UTDE21 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE09D. 
   Замечания по использованию
 

В общем случае заданная точность не гарантируется. При работе подпрограммы значения параметров M, XN, YN, XK, HMIN, EPS, P сохраняются.

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

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

               y1' = - y1 - 5y2 ,     y1(0) = 1
               y2' =   y1  + y2 ,      y2(0) = 1 ,    0 ≤ x ≤ 2

  Точное решение системы:
               y1  =  cos2x - 3 sin2x ,
               y2  =  cos2x +  sin2x

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

Unit TDE09R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE09R_p, DE09R_p;

function TDE09R: String;

implementation

function TDE09R: String;
var
M,IERR :Integer;
XN,XK,HMIN,EPS,P,H,Y1,Y2 :Real;
Y :Array [0..1] of Real;
YN :Array [0..1] of Real;
R1 :Array [0..9] of Real;
R2 :Array [0..1] of Real;
R3 :Array [0..1] of Real;
begin
Result := '';  { результат функции }
M := 2;
XN := 0.0;
YN[0] := 1.0;
YN[1] := 1.0;
ХК := 1.0;
HMIN := 1.E-12;
EPS := 1.E-10;
P := 100.0;
H := 0.01;
DE09R (FDE09R,M,XN,YN,XK,HMIN,EPS,P,H,Y,R1,R2,R3,IERR);
Y1 := Cos(2.0) - 3.0*Sin(2.0);
Y2 := Cos(2.0) + Sin(2.0);
Result := Result + Format('          %20.16f           %20.16f ',
 [Y[0],Y[1]]) + #$0D#$0A;
Result := Result + Format('          %20.16f           %20.16f ',
 [Y1,Y2]) + #$0D#$0A;
UtRes('TDE09R',Result);  { вывод результатов в файл TDE09R.res }
exit;
end;

end.

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

procedure fde09r(X :Real; var Y :Array of Real; var Z :Array of Real;
                M :Integer);

implementation

procedure fde09r(X :Real; var Y :Array of Real; var Z :Array of Real;
                M :Integer);
begin
Z[0] := -Y[0] - 5.0*Y[1];
Z[1] := Y[0] + Y[1];
end;

end.

 Результаты: 

      IERR = 0 
      Y(1) = -3.144039116632 + 00       Y(2) = 4.931505897994-01
        Y1 = -3.144039117018                  Y2 = 4.931505902741-01