Подпрограмма:  DE34R (модуль DE34R_p) (версия DE18R)

Назначение

Построение начальных значений при интегрировании системы обыкновенных дифференциальных уравнений второго порядка методом Штермера с контролем точности.

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

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

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

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

          Y (XN) = YN,      YN = ( y10, ..., yM0 )   ,
          Y '(XN) = DYN,  DYN = ( y10', ..., yM0' )   ,

отыскиваются необходимые при численном интегрировании этой системы многошаговым методом Штермера значение решения  Y в первом (после начальной точки  XN) узле интегрирования  ХN + Н, его конечная разность первого порядка и конечные разности правой части системы до четвертого порядка включительно в том же узле  XN + H.

Четвертый порядок разностей соответствует шестому порядку точности используемого метода Штермера.

Указанная процедура вычисления необходимых для счета начальных значений называется разгоном.

B качестве критерия точности при разгоне принята оценка первого отброшенного члена экстраполяционной формулы Адамса - Штермера, при этом все компоненты решения проверяются на точность по меpе погрешности, который заключается в следующем: если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы  P, называемой границей перехода, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.

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

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

procedure DE34R(F :Proc_F_DE; M :Integer; IORDER :Integer;
                XN :Real; var YN :Array of Real;
                var DYN :Array of Real; HMIN :Real; EPS :Real;
                P :Real; var H :Real; var X :Real;
                var YX :Array of Real; var DY :Array of Real;
                var DF :Array of Real; var RFN :Array of Real;
                var RF :Array of Real; var R :Array of Real;
                var IERR :Integer);

   Параметры 

   Версии 

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

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

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

       y1'' = - y1   , 
       y2'' = - y2   ,
       y1 (3/4 * π) = √2 / 2,        y2 (3/4 * π) = - √2 / 2   ,
       y1' (3/4 * π) = - √2 / 2,     y2' (3/4 * π) = - √2 / 2   .

Unit TDE34R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE34R_p, DE34R_p;

function TDE34R: String;

implementation

function TDE34R: String;
var
IORDER,M,_i,J,I,IERR :Integer;
H,P,EPS,HMIN,XN,X :Real;
YN :Array [0..1] of Real;
DYN :Array [0..1] of Real;
Y :Array [0..1] of Real;
DY :Array [0..1] of Real;
DELTY :Array [0..1] of Real;
DF :Array [0..9] of Real;
RFN :Array [0..1] of Real;
RF :Array [0..1] of Real;
begin
Result := '';  { результат функции }
H := 0.01;
IORDER := 5;
P := 100.0;
EPS := 0.0001;
HMIN := 1.E-10;
M := 2;
XN := 0.75*3.14159265359;
YN[0] := Sqrt(2.0)/2.0;
YN[1] := -YN[0];
DYN[0] := YN[1];
DYN[1] := -YN[0];
DE34R(FDE34R,M,IORDER,XN,YN,DYN,HMIN,EPS,P,H,X,Y,
     DY,DF,RFN,RF,DELTY,IERR);
RЕSult := Result + Format(' %20.16f ',[X]) + #$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 + #$0D#$0A;
FОR _i:=0 to 1 do
 begin
  Result := Result + Format('%20.16f ',[DY[_i]]);
  if ( ((_i+1) mod 2)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for J:=1 to 5 do
 begin
  for I:=1 to 2 do
   BЕgin
    Result := Result + Format('%20.16f ',[DF[(I-1)+(J-1)*2]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
UtRes('TDE34R',Result);  { вывод результатов в файл TDE34R.res }
exit;
end;

end.

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

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

implementation

procedure fde34r(T :Real; var Y :Array of Real; var Z :Array of Real;
                M :Integer);
begin
Z[0] := -Y[0];
Z[1] := -Y[1];
end;

end.

Результаты:

       IERR  =  0

       X  =  2.366194490 + 00

       Y(1)   =   7.000004762 - 01        Y(2)   =  - 7.141423761 - 01

       DY(1)  =  - 7.106305006 - 03        DY(2)  =  - 7.035594917 - 03

       DF(1, 1)  =  - 5.83333730151 - 06       DF(1, 2)  =     5.92192083634 - 08
       DF(1, 3)  =    5.89251453187 - 10        DF(1, 4)  =  - 5.86490578325 - 12
       DF(1, 5)  =  - 5.76760861293 - 14

       DF(2, 1)  =    5.95118646749 - 06       DF(2, 2)  =     5.86299576250 - 08
       DF(2, 3)  =  - 5.89250072347 - 10       DF(2, 4)  =  - 5.92374760355 - 12
       DF(2, 5)  =    6.01324545713 - 14