Подпрограмма:  DE40R (модуль DE40R_p)

Назначение

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

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

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

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

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

По значению решения  YX в узле  X_n вычиляется значение этого решения в узле  X_n + H. При этом, если  X_n является началом   ХN интервала интегрирования, то для вычисления решения в   XN + H требуется еще и знание производной решения в  XN.

Bместе с решением подпрограмма определяет также всю информацию, необходимую для вычисления приближенного решения в любом следующем узле. Значение  H шага интегрирования задается пользователем при обращении к подпрограмме.

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

Жоголев E.A. Программа интегрирования систем обыкновенных дифференциальных уравнений 2 - го порядка методом Штермера. Сб. "Вычислительные методы и программирование", вып. 1 , Изд - во МГУ, 1962.

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

procedure DE40R(F :Proc_F70; M :Integer; var JSTART :Integer;
                var YX :Array of Real; var X :Real;
                var DY :Array of Real; H :Real;
                var DELTY :Array of Real; var DF :Array of Real;
                var RFN :Array of Real; var RF :Array of Real;
                var YP :Array of Real; var DYP :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 TDE40R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE40R_p, DE40R_p;

function TDE40R: String;

implementation

function TDE40R: String;
var
M,JSTART,IH,IERR :Integer;
X,H,Y1,Y2 :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;
YP :Array [0..1] of Real;
DYP :Array [0..1] of Real;
label
_6,_7,_8,_20;
begin
Result := '';  { результат функции }
M := 2;
Х := 0.75*3.14159265359;
H := 0.01;
Y[0] := Sqrt(2.0)/2.0;
Y[1] := -Y[0];
DY[0] := Y[1];
DY[1] := -Y[0];
JSTART := 0;
IH := 0;
//H := -H;
_6:
DE40R(FDE40R,M,JSTART,Y,X,DY,H,DELTY,DF,RFN,RF,YP,DYP,IERR);
IH := IH+1;
Y1 := Sin(X);
Y2 := Cos(X);
Result := Result + Format(' %20.16f %20.16f %20.16f %20.16f %20.16f %20.16f ',
 [X,Y[0],Y[1],H,Y1,Y2]) + #$0D#$0A;
if ( IH = 1 ) 
 then goto _6;
if ( IH = 5 ) 
 then goto _20;
JSTART := -1;
if ( IH = 4 ) 
 then goto _8;
if ( IH = 3 ) 
 then goto _7;
H := 0.005;
H := -H;
goto _6;
_7:
H := 0.03;
goto _6;
_8:
H := 0.01;
goto _6;
_20:
UtRes('TDE40R',Result);  { вывод результатов в файл TDE40R.res }
exit;
end;

end.

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

procedure fde40r(X :Real; var Y :Array of Real; var DY :Array of Real;
                M :Integer);

implementation

procedure fde40r(X :Real; var Y :Array of Real; var DY :Array of Real;
                M :Integer);
begin
DY[0] := -Y[0];
DY[1] := -Y[1];
end;

end.

Результаты:

после первого обращения к подпрограмме -

                  X                             Y(1)                                Y(2)
       2.366194490 + 00      7.000004762 - 01        - 7.141423761 - 01

                  H                             Y1                                  Y2
       1.000000000 - 02      7.000004762 - 01        - 7.141423761 - 01

после второго обращения к подпрограмме -
 
                  X                             Y(1)                                Y(2)
       2.376194490 + 00      6.928241717 - 01        - 7.211065574 - 01

                  H                             Y1                                  Y2
       1.000000000 - 02      6.928241717 - 01        - 7.211065574 - 01

после третьего обращения к подпрограмме -

                  X                             Y(1)                                Y(2)
       2.361194490 + 00      7.035624232 - 01        - 7.106334615 - 01

                  H                             Y1                                  Y2
     - 5.000000000 - 03      7.035624232 - 01        - 7.106334615 - 01

после четвертого обращения к подпрограмме -
  
                  X                             Y(1)                                 Y(2)
       2.386194490 + 00      6.855785854 - 01        - 7.279986286 - 01

                  H                             Y1                                   Y2
       2.000000000 - 02      6.855785854 - 01        - 7.279986286 - 01

после пятого обращения к подпрограмме -

                  X                             Y(1)                                 Y(2)
       2.376194490 + 00      6.928241717 - 01        - 7.211065574 - 01

                  H                             Y1                                   Y2
       1.000000000 - 02      6.928241717 - 01        - 7.211065574 - 01