Текст подпрограммы и версий
de14r_p.zip , de14e_p.zip 		Тексты тестовых примеров
tde14r_p.zip , tde14e_p.zip

Подпрограмма:  DE14R (модуль DE14R_p)

Назначение

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

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

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

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

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

классическим методом Рунге - Кутта четвертого порядка точности с постоянным шагом. Решение вычисляется в одной точке ХК, которая является концом интервала интегрирования.

Березин И.С., Жидков Н.П. Методы вычислений, т.2. Физматгиз, М., 1960.

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

procedure DE14R(F :Proc_F70; var M :Integer; XN :Real;
                var YN :Array of Real; XK :Real; var H :Real;
                var Y :Array of Real; var A :Array of Real;
                var B :Array of Real; var C :Array of Real);


   Параметры
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 - начальные значения аргумента и решения; в случае системы уравнений (т.е. M ≠ 1) YN представляет одномерный массив длины M (тип: вещественный);
XK - значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования); XK может быть больше, меньше или pавно XN (тип: вещественный);
H - вещественная переменная, содержащая значение шага интегрирования. Может задаваться с учетом направления интегрирования, т.е. положительным, если XK > XN, отрицательным, если  XK < XN , или без такого учета с любым знаком " + " или " - ";
Y - искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента XK. Для системы уравнений (когда M ≠ 1) задается одномерным массивом длины M. B случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: вещественный);
A, B, C - одномерные вещественные рабочие массивы длиной M. 
   Версии
DE14E - вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования классическим методом Рунге - Кутта 4 - го порядка точности с постоянным шагом, при этом все вычисления над действительными числами выполняются с удвоенным числом значащих цифр. B этом случае параметры XN, YN, XK, H, Y, A, B, C и параметры X, Y и DY в подпрограмме F должны иметь тип Extended. 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
 

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

При работе подпрограммы и ее версии значения параметров M, XN, YN, XK сохраняются. Значение параметра H сохpаняется, если он задан с учетом направления интегрирования, иначе его знак меняется на противоположный. Если после работы подпрограммы нет необходимости иметь начальное значение решения 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 TDE14R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE14R_p, DE14R_p;

function TDE14R: String;

implementation

function TDE14R: String;
var
M,_i :Integer;
XN,XK,H :Real;
YN :Array [0..3] of Real;
Y :Array [0..3] of Real;
A :Array [0..3] of Real;
B :Array [0..3] of Real;
C :Array [0..3] of Real;
begin
Result := '';  { результат функции }
M := 4;
XN := 0.0;
ХК := 8.0;
YN[0] := 0.0;
YN[1] := 0.0;
YN[2] := 0.0;
YN[3] := 0.0;
H := 0.01;
DE14R(FDE14R,M,XN,YN,XK,H,Y,A,B,C);
Result := Result + Format('%s',[' XK']);
Result := Result + Format('%20.16f ',[XK]);
Result := Result + Format('%s',[#$0D#$0A + '  Y']);
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('TDE14R',Result);  { вывод результатов в файл TDE14R.res }
exit;
end;

end.

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

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

implementation

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

end.

Результаты:

          Y(1)  =  9.23847083543 - 03
          Y(2)  =  4.82097225727 - 03
          Y(3)  =  1.66711330351 + 00
          Y(4)  =  1.88434937053 - 02