Подпрограмма:  DE00R (модуль DE00R_p) (версия DE84R)

Назначение

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

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

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

(1)                             Y '  =   F(X, Y) ,

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

Предполагается, что среди характеристических корней матрицы Якоби ∂F / ∂y функции F имеются большие по модулю корни. По заданному значению решения YX в узле  x_n  вычисляется значение этого решения в узле  x_n + H . Предполагается также, что система (1) на [ x_n, x_N + H ] может быть достаточно хорошо аппроксимирована линейной однородной системой с постоянными коэффициентами, т.е. системой вида

                                     Y '  =   AY  , 

где A - некоторая постоянная матрица. Вычисление производится по методу Лоусона. Метод Лоусона заключается в следующем. Исходная система уравнений с помощью замены искомой функции  y (x) на [ x_n, x_N + H ] по формуле

                             y(x)  =  exp [ ( x - xn ) A0 ] Z(x)  ,

где A⁰ - некоторая постоянная матрица, преобразуется в систему уравнений относительно новой неизвестной функции Z (X) :

(2) Z ' (x)  =  F1(x, z)  =  exp [ - ( x - xn ) A0 ] { F(x, exp [ ( x - xn ) A0 ] Z(x)) -
                                      - A0 exp [ ( x - xn ) A0 ] Z(x) } 

Матрица Якоби ∂ F1 / ∂Z системы (2) и матрица Якоби ∂F / ∂y системы (1) связаны между собой соотношением

         ∂ F1 / ∂Z  =  exp [ - (x - xn) A0] { ∂F / ∂y - A0 } exp [ ( x - xn ) A0 ] 

Указанное преобразование выполняется самой подпрограммой. В качестве матрицы A⁰ подпрограмма выбирает матрицу ∂F / ∂y , вычисленную либо в точке ( x_n, y_n ), либо в некоторой предыдущей точке ( x_n - k, y_n - k ). Если шаг H достаточно мал, то данное преобразование позволяет уменьшить характеристические корни матрицы ∂ F1 / ∂Z по сравнению с характеристическими корнями матрицы ∂F / ∂y . А это приводит к уменьшению константы Липшица системы (2) по сравнению с константой Липшица системы (1).

Для решения системы (2) применяются классический метод Рунге - Кутта четвертого порядка точности, причем одновременно с решением (2) производится обратное преобразование от функции Z (x) к функции Y (x).

Все компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Абсолютная погрешность приближенного решения оценивается по правилу Рунге.

J.Douglas Lowson, Generalized Runge - Kutta processes for stable systems with large lipshitz constants, SIAM Journal on Numerical Analisys, 1967, Vol.4, No 3.

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

procedure DE00R(F :Proc_F_DE; FJ :Proc_F_DE; M :Integer;
                var JSTART :Integer; HMIN :Real; EPS :Real; P :Real;
                var YX :Array of Real; var X :Real; var H :Real;
                var BUL :Boolean; var XP :Real; var YP :Array of Real;
                var A :Array of Real; var E1 :Array of Real;
                var E2 :Array of Real; var E3 :Array of Real;
                var E4 :Array of Real; var R :Array of Real;
                var R0 :Array of Real; var R1 :Array of Real;
                var R2 :Array of Real; var R3 :Array of Real;
                var YD :Array of Real; var IERR :Integer);

   Параметры 

   Версии 

procedure DE84R(F :Proc_F_DE; M :Integer; var JSTART :Integer;
                HMIN :Real; EPS :Real; P :Real; var YX :Array of Real;
                var X :Real; var H :Real; var BUL :Boolean;
                var XP :Real; var YP :Array of Real;
                var A :Array of Real; var E1 :Array of Real;
                var E2 :Array of Real; var E3 :Array of Real;
                var E4 :Array of Real; var R :Array of Real;
                var R0 :Array of Real; var R1 :Array of Real;
                var R2 :Array of Real; var R3 :Array of Real;
                var YD :Array of Real; var IERR :Integer);
implementation

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

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

               Y '  =  AY ,
  где   A - некоторая постоянная матрица. 

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

     y1'  =  (- 500 + 2 x) y1  +  (1 - x / 10 )y2  +  (1 + 1 / ( x+1 ) ) sin y1 ,
     y1(   = 49 
     y2'  =  (- 1000 - 10 x) y1  +  (1  + x / 10 ) y2  +  (2  + 1 / ( x+2 ) ) cos y2 ,
     y2(   = 50 

              1 ≤ X ≤ 6 

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

Unit TDE00R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE00R_p, FJDE00R_p, DE00R_p;

function TDE00R: String;

implementation

function TDE00R: String;
var
M,JSTART,IH,IERR :Integer;
X,HMIN,EPS,P,H,ХР :Real;
BUL :Boolean;
YX :Array [0..1] of Real;
YP :Array [0..1] of Real;
A :Array [0..3] of Real;
E1 :Array [0..3] of Real;
E2 :Array [0..3] of Real;
E3 :Array [0..3] of Real;
E4 :Array [0..3] of Real;
R :Array [0..1] of Real;
R0 :Array [0..1] of Real;
R1 :Array [0..1] of Real;
R2 :Array [0..1] of Real;
R3 :Array [0..1] of Real;
YD :Array [0..1] of Real;
label
_100,_101,_102,_103,_104,_105,_106;
begin
Result := '';  { результат функции }
M := 2;
X := 1.0;
YX[0] := 49.0;
YX[1] := 50.0;
HMIN := 1.E-10;
EPS := 1.E-8;
P := 1.E2;
JSTART := 0;
H := 0.01;
IH := 0;
_100:
IH := IH+1;
DE00R (FDE00R,FJDE00R,M,JSTART,HMIN,EPS,P,YX,X,H,BUL,XP,YP,A,E1,
     E2,E3,E4,R,R0,R1,R2,R3,YD,IERR);
Result := Result + Format('     %20.16f      %20.16f      %20.16f ',
 [X,YX[0],YX[1]]) + #$0D#$0A;
Result := Result + Format('%s',['-']) + #$0D#$0A; 
Result := Result + Format('     %20.16f      %20.16f      %20.16f ',
 [H,R1[0],R1[1]]) + #$0D#$0A;
Result := Result + Format('%s',['**']) + #$0D#$0A; 
case IH of
 1: goto _101;
 2: goto _102;
 3: goto _103;
 4: goto _104;
 5: goto _105;
 6: goto _106;
end;
_101:
H := 1.0;
JSTART := 3;
goto _100;
_102:
EPS := 1.E-2;
goto _100;
_103:
JSTART := -1;
EPS := 1.E-8;
goto _100;
_104:
EPS := 1.E-2;
H := 1.E-8;
JSTART := 3;
goto _100;
_105:
H := -1.E-3;
EPS := 1.E-8;
JSTART := 0;
goto _100;
_106:
UtRes('TDE00R',Result);  { вывод результатов в файл TDE00R.res }
exit;
end;

end.

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

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

procedure FDE00R(X :Real; var Y :Array of Real; var DY :Array of Real;
                M :Integer);
var
E2,E,E3,E1 :Real;
begin
DY[1-1] := (2.0*X-5.E2)*Y[1-1]+(1.0-X*1.E-1)*Y[2-1]+
     (1.0+1./(1.0+X))*Sin(Y[1-1]);
DY[2-1] := (-1.E3-1.E1*X)*Y[1-1]+(1.0+X*1.E-1)*Y[2-1]+
     (2.0+1./(2.0+X))*Cos(Y[2-1]);
end;

end.

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

procedure FJDE00R(X :Real; var Y :Array of Real;
                var DF :Array of Real; M :Integer);
implementation

procedure FJDE00R(X :Real; var Y :Array of Real;
                var DF :Array of Real; M :Integer);
var
E2,E,E3,E1 :Real;
begin
DF[(1-1)+(1-1)*2] := 2.0*X-5.E2+(1.0+1./(1.0+X))*Cos(Y[1-1]);
DF[(1-1)+(2-1)*2] := 1.0-X*1.E-1;
DF[(2-1)+(1-1)*2] := -1.E3-1.E1*X;
DF[(2-1)+(2-1)*2] := 1.0+X*1.E-1-(2.0+1./(2.0+X))*Sin(Y[2-1]);
end;

end.

   Результаты: 

   После первого обращения к подпрограмме
             X                                    YX(1)                             YX(2) 
      1.000039062499+00       4.895766071531+00       4.808757932455+00 
             H                                     R1(1)                             R1(2) 
      3.906250000002-05        9.597251282538-10       -7.233271996176-09 

   После второго обращения к подпрограмме
             X                                    YX(1)                             YX(2) 
      1.000069580077+00       4.733408100787+01       4.661901490984+01 
             H                                     R1(1)                             R1(2) 
      3.051757812500-05        1.319373647371-10        6.946114202335-02 

   После третьего обращения к подпрограмме
             X                                    YX(1)                             YX(2) 
      1.000100097655+00       4.662140584405+01       4.517257692956+01 
             H                                     R1(1)                             R1(2) 
      6.103515625000-05       -6.208817164108-11        1.830362084864-09 

   После четвертого обращения к подпрограмме
             X                                    YX(1)                             YX(2) 
      1.000100097655+00       4.662140584411+01       4.517257692973+01 
             H                                     R1(1)                             R1(2) 
      3.051757812500-05       -5.820766091346-11        1.851003617047-09 

   После пятого обращения к подпрограмме
             X                                    YX(1)                             YX(2) 
      1.000100107654+00       4.662117408356+01       4.517210655950+01 
             H                                     R1(1)                             R1(2) 
      2.000000000004-08        0.000000000000 + 00      0.000000000000 + 00 

   После шестого обращения к подпрограмме
             X                                    YX(1)                             YX(2) 
      1.000068857653+00       4.735108317091+01       4.665352269082+01 
             H                                     R1(1)                             R1(2) 
     -3.125000000004-05        3.880510727560-11       -2.250696221985-09