Подпрограмма:  DE43R (модуль DE43R_p) (версия DE49R)

Назначение

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

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

Решается задача Коши для системы  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' )   ,   -

многошаговым предсказывающе - исправляющим методом пятого порядка точности.

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

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

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

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

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

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

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

procedure DE43R(F :Proc_F_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 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'' = - 6y1 - 7y2   ,
     y2'' = - 3y1 - 2y2 + 2x   ;

     y1 (0) = 0   ,           y2 (0) = 0   , 
     y1' (0) = - 1 / 9   ,   y2' (0) = - 3   . 

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

Unit TDE43R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE43R_p, DE43R_p;

function TDE43R: String;

implementation

function TDE43R: String;
var
M,K,IERR :Integer;
H,HMIN,XN,XK,EPS,Y3,Y1,Y2,P,НМАХ :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;
label
_200,_20;
begin
Result := '';  { результат функции }
H := 0.01;
HMIN := 1.E-10;
M := 2;
XN := 0.0;
ХК := 5.0;
YN[0] := 0.0;
YN[1] := 0.0;
DYN[0] := -1.0/9.0;
DYN[1] := -3.0;
EPS := 0.0001;
Y3 := (Exp(5.0)-Exp(-5.0))/3.0;
Y1 := Y3-7.0*Sin(15.0)/9.0+70.0/9.0;
Y2 := -Y3-Sin(15.0)/3.0-20.0/3.0;
P := 100.0;
НМАХ := 0.1;
for K:=1 to 2 do
 begin
  DE43R(FDE43R,M,XN,YN,DYN,XK,HMIN,HMAX,EPS,P,H,Y,DY,DELTY,
       DF,RFN,RF,YP,DYP,IERR);
  Result := Result + Format(' %20.16f %20.16f %20.16f %20.16f %20.16f %20.16f ',
 [XK,Y[0],Y[1],H,Y1,Y2]) + #$0D#$0A;
  if ( IERR <> 0 ) 
   then goto _20;
  Y1 := -Y1;
  Y2 := -Y2;
  ХК := -5.0;
  H := 0.01;
_200:
 end;
_20:
UtRes('TDE43R',Result);  { вывод результатов в файл TDE43R.res }
exit;
end;

end.

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

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

implementation

procedure fde43r(T :Real; var Y :Array of Real; var Z :Array of Real;
                M :Integer);
begin
Z[0] := -6.e0*Y[0]-7.e0*Y[1];
Z[1] := -3.e0*Y[0]-2.e0*Y[1]+2.e0*T;
end;

end.

Результаты:

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

                   XK                           Y(1)                              Y(2)
       5.000000000 + 00      5.674103844 + 01      - 5.655217994 + 01

                   H                              Y1                                Y2
       1.000000000 - 01      5.674080540 + 01      - 5.635223633 + 01

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

                   XK                            Y(1)                             Y(2)
     - 5.000000000 + 00    - 5.674103838 + 01       5.635217988 + 01

                   H                               Y1                                Y2
     - 1.000000000 - 01    - 5.674080540 + 01       5.635223633 + 01