Подпрограмма:  DE47R (модуль DE47R_p)

Назначение

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

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

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

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

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

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

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

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

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

Решение и производная вычисляются в одной точке  XK, которая является концом интервала интегрирования.

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

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

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

procedure DE47R(F :Proc_F2_DE; M :Integer; XN :Real;
                var YN :Array of Real; var DYN :Array of Real;
                XK :Real; HMIN :Real; HMAX :Real; EPS :Real;
                P :Real; var H :Real; var Y :Array of Real;
                var DY :Array of Real; var Z :Array of 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 ZP :Array of Real; var IERR :Integer);

   Параметры 

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

   Версии 

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

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

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

     y1'' = 0.5 ( y2 - y1 + y1' + y2' ) - 0.5   ,   y1 (0) = 1   ,   y1' (0) = 1.5   ,
     y2 = y1'' - 0.5 x   ,   y2 (0) = - 1   ,   y2 (0) = - 0.5   .

Unit TDE47R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE47R_p, DE47R_p;

function TDE47R: String;

implementation

function TDE47R: String;
var
M,_i,IERR :Integer;
XN,XK,HMIN,HMAX,EPS,P,H,Y1,Y2,DY1,DY2 :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;
YP :Array [0..1] of Real;
DYP :Array [0..1] of Real;
Z :Array [0..1] of Real;
ZP :Array [0..1] of Real;
begin
Result := '';  { результат функции }
M := 2;
XN := 0.0;
YN[0] := 1.0;
YN[1] := -1.0;
DYN[0] := 1.5;
DYN[1] := -0.5;
ХК := 5.0;
HMIN := 1.E-10;
НМАХ := 0.1;
EPS := 0.0001;
P := 100.0;
H := 0.01;
DE47R(FDE47R,M,XN,YN,DYN,XK,HMIN,HMAX,EPS,P,H,
     Y,DY,Z,DELTY,DF,RFN,RF,YP,DYP,ZP,IERR);
Result := Result + Format('%s',[' IERR=']);
Result := Result + Format('%3d ',[IERR]) + #$0D#$0A;
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 1)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 1 do
 begin
  Result := Result + Format('%20.16f ',[DY[_i]]);
  if ( ((_i+1) mod 1)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Y1 := Sin(XK)+Cos(XK)+0.5*XK;
Y2 := -Y1+XK;
DY1 := Cos(XK)-Sin(XK)+0.5;
DY2 := -DY1+1.0;
Result := Result + Format('%20.16f %20.16f %20.16f %20.16f %20.16f ',
 [XK,Y1,Y2,DY1,DY2]) + #$0D#$0A;
UtRes('TDE47R',Result);  { вывод результатов в файл TDE47R.res }
exit;
end;

end.

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

procedure FDE47R(X :Real; var Y :Array of Real;
                var DY :Array of Real; var Z :Array of Real;
                M :Integer);
implementation

procedure FDE47R(X :Real; var Y :Array of Real;
                var DY :Array of Real; var Z :Array of Real;
                M :Integer);
begin
Z[0] := 0.5*(Y[1]+DY[0]-Y[0]+DY[1])-0.5;
Z[1] := Y[0]-0.5*X;
end;

end.

Результаты:

      IERR = 0

                      XK                              Y(1)                             Y(2)
       5.000000000000 + 00     1.824737966524 + 00    3.175262014483 + 00

                                                          Y1                                Y2
                                              1.824737910809 + 00      3.175262089189 + 00

                                                         DY(1)                            DY(2)
                                              1.742586640374 + 00    - 7.425866599151 - 01

                                                         DY1                              DY2
                                              1.742586460128 + 00    - 7.425864601282 - 01