Подпрограмма:  DE16R (модуль DE16R_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 качестве абсолютной погрешности решения используется оценка главного члена асимптотического разложения погрешности метода на одном шаге, получаемая при вычитании двух выражений, представляющих приближенные значения решения пятого и четвертого порядка точности. При этом на каждом шаге интегрирования для определения решения и его погрешности используются всего шесть вычислений правой части системы.

Дж.Форсайт, М.Малькольм, К.Моулер, Машинные методы математических вычислений, "Мир", M., 1980.

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

procedure DE16R(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 RAB :Array of Real;
                var IERR :Integer);

   Параметры 

   Версии 

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

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

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

          y1'  =  -y2 - 0.1 x - 0.9
          y2'  =   y1 - 0.1 x - 1.1     0 ≤ x ≤ π
    
          y1 (0)  =   1 ,     y2 (0)  =  - 2
   
 Точное решение системы:

         y1  =    0.1 x  +  sin x  +  1.
         y2  =  - 0.1 x - cos x -    

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

Unit TDE16R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE16R_p, DE16R_p;

function TDE16R: String;

implementation

function TDE16R: String;
var
M,_i,IERR :Integer;
XN,XK,HMIN,EPS,P,H,Y1,Y2 :Real;
YN :Array [0..1] of Real;
Y :Array [0..1] of Real;
RАВ :Array [0..12] of Real;
begin
Result := '';  { результат функции }
M := 2;
XN := 0.0;
YN[0] := 1.0;
YN[1] := -2.0;
ХК := 3.141592653589793238;
HMIN := 1.E-14;
EPS := 1.E-5;
P := 100.0;
H := 0.01;
DE16R(FDE16R,M,XN,YN,XK,HMIN,EPS,P,H,Y,RAB,IERR);
Result := Result + Format('%s',[' IERR=']);
Result := Result + Format('%4d ',[IERR]) + #$0D#$0A;
Y1 := 0.1*XK+Sin(XK)+1.0;
Y2 := -0.1*XK-Cos(XK)-1.0;
Result := Result + Format('%20.16f ',[XK]) + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 1 do
 begin
  Result := Result + Format('%20.16f ',[Y[_i]]);
  if ( ((_i+1) mod 2)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%20.16f  %20.16f  %20.16f ',
 [H,Y1,Y2]) + #$0D#$0A;
UtRes('TDE16R',Result);  { вывод результатов в файл TDE16R.res }
exit;
end;

end.

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

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

implementation

procedure fde16r(X :Real; var Y :Array of Real; var Z :Array of Real;
                M :Integer);
var
R :Real;
begin
R := 0.1*X;
Z[0] := -Y[1]-R-0.9;
Z[1] := Y[0]-R-1.1;
end;

end.

Результаты:

          IERR  =   0

          Y(1)  =   1.314159265929 + 00
          Y(2)  =  -3.141592563520 - 01