Подпрограмма:  DE21R (модуль DE21R_p) (версия DE24R)

Назначение

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

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

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

          Y' = F (X, Y) ,

          Y = ( y1, ..., yM ) ,
          F = ( f1 (X, y1,..., yM),..., fM (X, y1,..., yM) )     методом Гира.

Метод Гира для нежесткой системы является многошаговым предсказывающе - исправляющим методом Адамса, записанным в форме Нордсика, при этом предсказание и исправление имеют один и тот же порядок.

B случае, когда система уравнений является жесткой, интегрирование осуществляется специальным методом основанным на методе типа Адамса и использующим якобиан ( ∂F/∂Y ) системы, который вычисляется подпрограммой по формулам численного дифференцирования.

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

                M   
             (  ∑ ( eI / YPM(I) )2 )1/2  ≤  EPS  ,
                I=1 

где e_I - вычисляемая подпрограммой оценка абсолютной погрешности приближенного значения I - й компоненты решения в узле X_n + H, YPM (I) - задается при обращении к подпрограмме и, по существу, позволяет использовать разный тип погрешности для каждой компоненты. Если, например, YPM(I) = 1, то для I - ой компоненты YX(I) решения будет использоваться абсолютная погрешность; если YPM (I) = | YX (I) | ≠ 0 , где YX(I) - значение I - й компоненты YX(I) в узле X_n, то для этой компоненты берется отношение абсолютной погрешности приближенного ее значения в узле X_n + H к абсолютной величине приближенного значения этой компоненты в предыдущем узле Х_n. Значение H может быть меньше или pавно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.

Gear C.W., The automatic integration of ordinary differential equations. Communications of the ACM, 14, 3 (March 1971), 176-179.

Gear C.W., Numerical initial value problems in ordinary differential equations, Prentice - Hall, Englewood Cliffs, N.J., 1971.

Gear C.W., The automatic integration of stiff ordinary differential equations. Information Processing 68, A.J.H.

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

procedure DE21R(F :Proc_F_DE; M :Integer; ISTIFJ :Integer;
                IORDER :Integer; var JSTART :Integer; HMIN :Real;
                HMAX :Real; EPS :Real; var YX :Array of Real;
                var X :Real; var H :Real; var BUL :Boolean;
                var YPM :Array of Real; var DELTY :Array of Real;
                var RAB :Array of Real; var YP :Array of Real;
                var IERR :Integer);

   Параметры 

   Версии 

      procedure DE24R(F :Proc_F_DE; FJ :Proc_F_DE; M :Integer; IORDER :Integer;
                var JSTART :Integer; HMIN :Real; HMAX :Real;
                EPS :Real; var YX :Array of Real; var X :Real;
                var H :Real; var BUL :Boolean; var YPM :Array of Real;
                var DELTY :Array of Real; var RAB :Array of Real;
                var YP :Array of Real; var IERR :Integer);
      

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

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

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

         y1'  =   y3
         y2'  =   y4
         y3'  =  - ( 1 - 0.002 x ) y1 / ( y12  +  y22 )3/2
         y4'  =  - ( 1 - 0.002 x ) y2 / ( y12  +  y22 )3/2  ,

        y1(0)  =  0.09411764706 ,     y2(0)  =  0 ,     y3(0)  =  0 ,     y4(0)  =  4.5 

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

Unit TDE21R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE21R_p, DE21R_p;

function TDE21R: String;

implementation

function TDE21R: String;
var
M,IORDER,JSTART,ISTIFJ,I,_i,IERR :Integer;
X,H,HMAX,HMIN,EPS :Real;
BUL :Boolean;
RАВ :Array [0..67] of Real;
Y :Array [0..3] of Real;
DELTY :Array [0..3] of Real;
YРМ :Array [0..3] of Real;
YP :Array [0..31] of Real;
label
_10;
begin
Result := '';  { результат функции }
BUL := False;
M := 4;
X := 0.0;
H := 0.01;
Y[0] := 0.09411764706;
Y[1] := 0.0;
Y[2] := 0.0;
Y[3] := 4.5;
IORDER := 4;
JSTART := 0;
НМАХ := 0.1;
HMIN := 1.E-18;
EPS := 1.E-5;
ISTIFJ := 0;
for I:=1 to M do
 begin
_10:
  YPM[I-1] := 1.0;
 end;
DE21R(FDE21R,M,ISTIFJ,IORDER,JSTART,HMIN,HMAX,EPS,Y,X,H,
     BUL,YPM,DELTY,RAB,YP,IERR);
Result := Result + Format(' %20.16f  %20.16f ',[X,H]) + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 3 do
 begin
  Result := Result + Format('%20.16f ',[Y[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['  DELTY(1)']);
Result := Result + Format('%s',['          DELTY(2)']);
Result := Result + Format('%s',['          DELTY(3)']);
Result := Result + Format('%s',['          DELTY(4)' + #$0D#$0A]);
Result := Result + #$0D#$0A;
for _i:=0 to 3 do
 begin
  Result := Result + Format('%20.16f ',[DELTY[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TDE21R',Result);  { вывод результатов в файл TDE21R.res }
exit;
end;

end.

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

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

implementation

procedure fde21r(X :Real; var Y :Array of Real; var DY :Array of Real;
                M :Integer);
var
R2,R :Real;
begin
R2 := Y[0]*Y[0]+Y[1]*Y[1];
R2 := R2*Sqrt(R2);
R := -(1.0-0.002*X)/R2;
DY[0] := Y[2];
DY[1] := Y[3];
DY[2] := R*Y[0];
DY[3] := R*Y[1];
end;

end.

Результаты:

          X  =   4.057708603-05
          H  =   4.057708603-05

          Y(1)  =   9.411746119-02
          Y(2)  =   1.825965265-04
          Y(3)  =  -4.580764462-03
          Y(4)  =   4.499991113

          IERR  =   0