Подпрограмма:  DE13R (модуль DE13R_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 ) , 

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

О.Б.Арушанян, Стандартная программа решения задачи Коши для систем обыкновенных дифференциальных уравнений методом Рунге - Кутта, вып. 31, под общей редакцией В.В.Воеводина, НИВЦ МГУ, 1968.

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

procedure DE13R(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 YP :Array of Real; var DELTY :Array of Real;
                var YR :Array of Real; var DY :Array of Real;
                var IERR :Integer);

   Параметры 

   Версии 

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

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

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

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

          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 TDE13R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE13R_p, DE13R_p;

function TDE13R: String;

implementation

function TDE13R: String;
var
M,_i,IERR :Integer;
XI,XF,H,HMIN,ERR,P :Real;
YI :Array [0..3] of Real;
YF :Array [0..3] of Real;
YP :Array [0..3] of Real;
DELTY :Array [0..3] of Real;
RАВ :Array [0..3] of Real;
RAB1 :Array [0..3] of Real;
begin
Result := '';  { результат функции }
M := 4;
YI[0] := 0.0;
YI[1] := 0.0;
YI[2] := 0.0;
YI[3] := 0.0;
XI := 0.0;
XF := 8.0;
H := 0.01;
HMIN := 1.E-16;
ERR := 1.E-4;
P := 1.E-8;
DE13R(FDE13R,M,XI,YI,XF,HMIN,ERR,P,H,YF,YP,DELTY,RAB,RAB1,IERR);
Result := Result + Format('%s',[' XK=']);
Result := Result + Format('%20.16f ',[XF]);
Result := Result + Format('%s',[#$0D#$0A + ' Y=']);
Result := Result + #$0D#$0A;
for _i:=0 to 3 do
 begin
  Result := Result + Format('%20.16f ',[YF[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' IERR=']);
Result := Result + Format('%3d ',[IERR]) + #$0D#$0A;
UtRes('TDE13R',Result);  { вывод результатов в файл TDE13R.res }
exit;
end;

end.

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

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

implementation

procedure fde13r(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.

Результаты:

          YF(1)  =  0.00923847083381
          YF(2)  =  0.00482097150315
          YF(3)  =  1.66711319966
          YF(4)  =  0.0188434937933

          IERR  =  0