Подпрограмма:  DE23R (модуль DE23R_p) (версия DE25R)

Назначение

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

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

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

          Y' = F (X, Y) ,

          Y = ( y1, ..., yM ) ,
          F = ( f1 (X, y1,..., yM),..., fM (X, y1,..., yM) )
 с начальными условиями, заданными в точке XN:
          Y(XN) = YN ,     YN = ( y10,..., yM0 ) , 

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

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

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

                M   
             (  ∑ δI2 )1/2  ≤  EPS  ,
                I=1 

где δ_I - одна из погрешностей следующих типов: абсолютная, относительная или стандартная. При этом под относительной погрешностью приближенного значения I - й компоненты решения в узле  x_n подразумевается отношение абсолютной погрешности  e_I этого значения в узле  x_n к абсолютной величине значения I - й компоненты в предыдущем узле  x_n - 1, т.е.  e_I / | y_I^n - 1|, а под стандартной погрешностью - отношение  e_I / YPM (I), где

           YPM(I)  =  max ( | yI0 |, | yI1 |,..., | yIn-1 | ) 

Tип погрешности специфицируется пользователем при обращении к подпрограмме.

Gear C.W. The automatic integration of ordinary differential equations. Communicatuons 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 DE23R(F :Proc_F_DE; M :Integer; XN :Real;
                var YN :Array of Real; XK :Real; HMIN :Real;
                HMAX :Real; EPS :Real; ISTIFJ :Integer;
                IORDER :Integer; IU :Integer; var H :Real;
                var Y :Array of Real; var YPM :Array of Real;
                var DELTY :Array of Real; var RAB :Array of Real;
                var YP :Array of Real; var IERR :Integer);

   Параметры 

   Версии 

      procedure DE25R(F :Proc_F_DE; FJ :Proc_F_DE; M :Integer; XN :Real;
                var YN :Array of Real; XK :Real; HMIN :Real;
                HMAX :Real; EPS :Real; IORDER :Integer;
                IU :Integer; var H :Real; var Y :Array of Real;
                var YPM :Array of Real; var DELTY :Array of Real;
                var RAB :Array of Real; var YP :Array of Real;
                var IERR :Integer);
      

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

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

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

 
         y1'  =  -20 y1  +  y2
         y2'  =   19 y1 - 2 y2 ,     0 ≤ x ≤ 1

         y1(0)  =  2 ,      y2(0)  =  18
 
 Собственные значения якобиевой матрицы системы
   λ1  =  -21 ,    λ2  =  -1 .
 Поэтому эту систему можно рассматривать как жесткую.
    
Unit TDE23R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE23R_p, DE23R_p;

function TDE23R: String;

implementation

function TDE23R: String;
var
IU,IORDER,M,ISTIFJ,IERR :Integer;
EPS,HMIN,XN,XK,HMAX,H,X :Real;
RАВ :Array [0..37] of Real;
Y :Array [0..1] of Real;
YN :Array [0..1] of Real;
DELTY :Array [0..1] of Real;
YРМ :Array [0..1] of Real;
YP :Array [0..15] of Real;
begin
Result := '';  { результат функции }
IU := 3;
EPS := 0.00001;
IORDER := 6;
HMIN := 1.E-10;
M := 2;
XN := 0.0;
YN[0] := 2.0;
YN[1] := 18.0;
ХК := 1.0;
ISTIFJ := 1;
НМАХ := 0.1;
H := 0.01;
DE23R(FDE23R,M,XN,YN,XK,HMIN,HMAX,EPS,ISTIFJ,IORDER,
     IU,H,Y,YPM,DELTY,RAB,YP,IERR);
Result := Result + Format(' %20.16f %20.16f %20.16f ',
 [X,Y[0],Y[1]]) + #$0D#$0A;
UtRes('TDE23R',Result);  { вывод результатов в файл TDE23R.res }
exit;
end;

end.

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

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

implementation

procedure fde23r(X :Real; var Y :Array of Real; var DY :Array of Real;
                M :Integer);
begin
DY[0] := -20.0*Y[0]+Y[1];
DY[1] := 19.0*Y[0]-2.0*Y[1];
end;

end.

Результаты:

          IERR  =  0

          Y(1)  =  3.678796482-01
          Y(2)  =  6.989709436