Подпрограмма:  DE42R (модуль DE42R_p) (версия DE48R)

Назначение

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

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

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

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

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

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

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

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

Бахвалов H.C., Численные методы, т.1, "Hаука", 1975.

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

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

procedure DE42R(F :Proc_F_DE; M :Integer; var JSTART :Integer;
                HMIN :Real; HMAX :Real; EPS :Real;
                P :Real; var YX :Array of Real; var X :Real;
                var DY :Array of Real; var H :Real; var BUL :Boolean;
                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 TDE42R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE42R_p, DE42R_p;

function TDE42R: String;

implementation

function TDE42R: String;
var
M,IH,JSTART,IERR :Integer;
H,HMIN,X,EPS,P,HMAX,Y1,Y2 :Real;
BUL :Boolean;
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,_53,_20;
begin
Result := '';  { результат функции }
H := -0.01;
HMIN := 1.E-10;
M := 2;
X := 0.75*3.14159265359;
Y[0] := Sqrt(2.0)/2.0;
Y[1] := -Y[0];
EPS := 0.0001;
DY[0] := Y[1];
DY[1] := -Y[0];
IH := 0;
BUL := False;
JSTART := 0;
P := 100.0;
НМАХ := 0.1;
_6:
DE42R(FDE42R,M,JSTART,HMIN,HMAX,EPS,P,Y,X,DY,H,BUL,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.02;
goto _6;
_8:
H := 0.01;
goto _6;
_53:
if ( IERR <> 0 ) 
 then goto _20;
_20:
UtRes('TDE42R',Result);  { вывод результатов в файл TDE42R.res }
exit;
end;

end.

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

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

implementation

procedure fde42r(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.346194490 + 00      7.141423761 - 01      - 7.000004762 - 01

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

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

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

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

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

                   X                           Y(1)                               Y(2)
       2.341194490 + 00      7.176334371 - 01      - 6.964210292 - 01

                   H                           Y1                                  Y2
     - 5.000000000 - 03      7.176334371 - 01      - 6.964210292 - 01

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

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

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

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

                   X                           Y(1)                                Y(2)
       2.356194490 + 00      7.071067812 - 01      - 7.071067812 - 01

                   H                           Y1                                   Y2
       1.000000000 - 02      7.071067812 - 01      - 7.071067812 - 01