Подпрограмма:  DE01R (модуль DE01R_p) (версия DE85R)

Назначение

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

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

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

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

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

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

                         Y(XN)  =  YN ,   YN  =  ( y1,..., yM )  . 

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

                                     Y '  =   AY  , 

где A - постоянная матрица, вообще говоря, своя для каждого отрезка. Для интегрирования системы применяется метод Лоусона.

Метод Лоусона является одношаговым методом и заключается в следующем. Допустим, что искомое решение системы (1) уже вычислено в некоторой точке  x = x_n  интервала интегрирования, т.е. известно  y_n ≈ y(x_n). Для отыскания решения  y(x_n + 1) = y(x_n + H) в следующем узле  x_n + 1 = x_n + H выполняются такие действия. Исходная система уравнений с помощью замены искомой функции  y (x) на  x_n ≤ x ≤ 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 DE01R(F :Proc_F_DE; FJ :Proc_F_DE; M :Integer; XN :Real;
                var YN :Array of Real; XK :Real; HMIN :Real;
                EPS :Real; P :Real; N :Integer; var H :Real;
                var Y :Array of Real; var R :Array of Real;
                var IERR :Integer);

   Параметры 

   Версии 

      procedure DE85R(F :Proc_F_DE; M :Integer; XN :Real;
                var YN :Array of Real; XK :Real; HMIN :Real;
                EPS :Real; P :Real; N :Integer; var H :Real;
                var Y :Array of Real; var R :Array of Real;
      

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

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

                     Y '  =  AY 

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

                                     
     y1'  =  (- 500 + 2 x) y1  +  (1 - x /10 ) y2  +  (1 + 1/(x +    ) 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 TDE01R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE01R_p, FJDE01R_p, DE01R_p;

function TDE01R: String;

implementation

function TDE01R: String;
var
M,N,IERR :Integer;
XN,HMIN,EPS,P,H,ХК :Real;
YN :Array [0..1] of Real;
Y :Array [0..1] of Real;
R :Array [0..34] of Real;
begin
Result := '';  { результат функции }
M := 2;
XN := 1.0;
YN[0] := 49.0;
YN[1] := 50.0;
HMIN := 1.E-10;
EPS := 1.E-8;
P := 100.0;
H := 0.01;
ХК := 6.0;
N := 2;
DE01R (FDE01R,FJDE01R,M,XN,YN,XK,HMIN,EPS,P,N,H,Y,R,IERR);
Result := Result + Format('     %20.16f      %20.16f      %20.16f      %1d ',
 [Y[0],Y[1],H,IERR]) + #$0D#$0A;
UtRes('TDE01R',Result);  { вывод результатов в файл TDE01R.res }
exit;
end;

end.

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

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

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

end.

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

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

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

end.
Результаты: 

             Y(1)                                    Y(2)                              H 
     -4.175740488392-02       -5.087148371158+01      2.500000000001-03 
      IERR = 0 

Результаты решения этой же задачи, полученные с помощью подпрограммы DE85R

     -4.175740472033-02       -5.087148363417+01      5.000000000003-03 
      IERR = 0

Результаты решения этой же задачи, полученные с помощью подпрограммы DE25R

     -4.175740429849-02       -5.087148203468+01      1.756831014061-03 
      IERR = 0