Подпрограмма:  DE05R (модуль DE05R_p)

Назначение

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

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

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

(1)                         Y ' (X)   =   A(X) * Y(X) + φ(X) ,

                                 Y   =   ( y1,..., yM ) , 
                                 A(X)   =   ( ai j(X) ) ,     i, j  =  1, ..., M ,
                                φ(X)    =   ( φ1(X), ... , φM(X) ) 

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

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

Предполагается, что среди характеристических корней матрицы A(X) имеются большие по модулю корни, а функция  φ (x) является достаточно малой. Решение вычисляется в одной точке XK, которая является концом интервала интегрирования. Для интегрирования системы применяется метод Лоусона.

Метод Лоусона является одношаговым методом и заключается в следующем. Допустим, что искомое решение системы (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)  =  A1(x) Z(x)  +  φ1(x)  =  exp [ - ( x - xn ) A0 ] { A(x) - A0 }
                                   exp [ ( x - xn ) A0 ] Z(x) + exp [ ( - ( x - xn ) A0 ] φ(x)

           xn ≤ x ≤ xn + H  

Данное преобразование выполняется самой подпрограммой. В качестве матрицы A⁰ подпрограмма выбирает матрицу A⁰ = A (x_n + H /2). Если шаг H достаточно мал, то преобразование позволяет уменьшить характеристические корни матрицы A1 (x) по сравнению с характеристическими корнями исходной матрицы A (x). Это приводит к уменьшению константы Липшица системы (2) по сравнению с константой Липшица системы (1). Для решения системы (2) применяются формулы классического метода Рунге - Кутта четвертого порядка точности, причем одновременно с решением (2) производится обратное преобразование от функции Z (x) к функции Y (x).

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

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

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

procedure DE05R(FA :Proc_F3_DE; FI :Proc_F3_DE; M :Integer; XN :Real;
                var YN :Array of Real; XK :Real; HMIN :Real;
                EPS :Real; P :Real; var H :Real; var Y :Array of Real;
                var R :Array of Real; var IERR :Integer);

   Параметры 

   Версии 

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

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

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

  1)
       y1'  =  - ( 2 + x ) y1 /( 1 + x )  +  20 x y2 ,
       y1(0) = 2 ,
       y2'  =  -20 x y1  +  ( 2  + x )  y2 /( 1 + x ) ,
       y2(0) = 18 

       0 ≤ x ≤ 6 

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

Unit TDE05R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FADE05R_p, FIDE05R_p, DE05R_p;

function TDE05R: String;

implementation

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

end.
 
Unit fade05r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;

procedure fade05r(var A :Array of Real; X :Real; M :Integer);

implementation

procedure fade05r(var A :Array of Real; X :Real; M :Integer);
begin
A[0] := -(2.0 + X)/(1.0 + X);
A[2] := 20.0*X;
A[1] := -A[2];
A[3] := A[0];
end;

end.

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

procedure fide05r(var R1 :Array of Real; X :Real; M :Integer);

implementation

procedure fide05r(var R1 :Array of Real; X :Real; M :Integer);
begin
R1[0] := 0.0;
R1[1] := 0.0;
end;

end.

   Результаты: 

             Y(1)                                Y(2)                             H 
      5.911150077338-03    -2.487346326713-03     1.600000000001-01 

      IERR = 0 
                                    
   2)
        y1'  =  - 20 x y1 + ( (1 + 2x)/(1 + 3x) ) y2  +  ( 1 / 10 ) x2 ,
        y1(0) = 22 
        y2'  =  19 x y1  +  ( (2 + x)/(1 + x) ) y2  -  ( 9 / 10 ) x2 ,
        y2(0) = 18 
              0 ≤ x ≤ 3 

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

 
Unit tde05r1_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FADE05R1_p, FIDE05R1_p, DE05R_p;

function tde05r1: String;

implementation

function tde05r1: String;
var
M,IERR :Integer;
XN,HMIN,EPS,P,XK,H :Real;
Y :Array [0..1] of Real;
YN :Array [0..1] of Real;
R :Array [0..36] of Real;
begin
Result := '';  { результат функции }
M := 2;
XN := 0.0;
YN[0] := 22.0;
YN[1] := 18.0;
HMIN := 1.E-10;
EPS := 1.E-10;
P := 100.0;
XK := 3.0;
H := 0.01;
DE05R (FADE05R1,FIDE05R1,M,XN,YN,XK,HMIN,EPS,P,H,Y,R,IERR);
Result := Result + Format(' %20.16f  %20.16f  %20.16f   %2d ',
 [Y[0],Y[1],H,IERR]) + #$0D#$0A;
UtRes('tde05r1',Result);  { вывод результатов в файл tde05r1.res }
exit;
end;

end.

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

procedure fade05r1(var A :Array of Real; X :Real; M :Integer);

implementation

procedure fade05r1(var A :Array of Real; X :Real; M :Integer);
begin
A[0] := -20.0*X;
A[2] := (1.0+2.0*X)/(1.0+3.0*X);
A[1] := 19.0*X;
A[3] := -(2.0+X)/(1.0+X);
end;

end.

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

procedure fide05r1(var R1 :Array of Real; X :Real; M :Integer);

implementation

procedure fide05r1(var R1 :Array of Real; X :Real; M :Integer);
begin
R1[0] := 0.1*X*X;
R1[1] := -9.0*R1[0];
end;

end.

   Результаты: 

             Y(1)                                 Y(2)                             H 
      2.134285516536-02     4.227925779755-01     2.500000000001-03 

      IERR = 0