Подпрограмма:  DE03R (модуль DE03R_p)

Назначение

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

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

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

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

                                   Y    =   ( y1,..., yM ) , 
                                   A    =   ( ai j ) ,     i, j  =  1, ..., M ,
                           U(X, Y)   =   ( U1(X, Y), ... , UM(X, Y) ) , 

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

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

где A - постоянная числовая матрица.

Предполагается, что среди характеристических корней матрицы A имеются большие по модулю корни, а константа Липшица для функции U (X, Y), т.е. константа L из условия Липшица

| Ui ( x, y1(1),..., yM(1) ) - Ui ( x, y1(2), ... , yM(2) ) |  ≤
                                                                                      M
                                                                             ≤  L  ∑   | yj(1) - yj(2)  |  , 
                                                                                     j=1 

независящая от  i,  x,  y⁽¹⁾,  y⁽²⁾ , невелика. Также предполагается, что нелинейный член U (x, y) является достаточно малым. Решение вычисляется в одной точке ХК, которая является концом интервала интегрирования. Для интегрирования системы применяется метод Лоусона.

Метод Лоусона является одношаговым и заключается в следующем. Допустим, что искомое решение системы (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 ) A ] Z(x)  ,

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

(2)   Z ' (x)  =  U1(x, z)  =  exp [ - ( x - xn ) A ]  U( x, exp [( x - xn ) A ] Z(x) ) 

Данное преобразование выполняется самой подпрограммой. Характеристические корни матрицы Якоби

   ∂U1 / ∂Z   =   exp [ - ( x - xn ) A ]  (∂U / ∂y)  exp [ ( x - xn ) A ] 

являются характеристическими корнями матрицы ∂U / ∂y и в силу малости константы Липшица функции U (x, y) невелики. Поэтому система (2) не жесткая и может быть решена традиционными методами численного интегрирования. Данная подпрограмма решает ее методом Рунге - Кутта  4 - ого порядка, при этом одновременно с решением (2) производится обратное преобразование от функции Z (x) к функции  y (x) .

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

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

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

procedure DE03R(FU :Proc_F_DE; M :Integer; XN :Real;
                var YN :Array of Real; XK :Real; var A :Array of Real;
                HMIN :Real; EPS :Real; P :Real; var H :Real;
                var Y :Array of Real; var R :Array of Real;
                var IERR :Integer);

   Параметры 

   Версии 

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

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

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

                                     
        y1'  =  - 500 y1  +  y2  +  sin ( y1 + y2 ) ( 1 + e - x/2 ) ,
        y1(0) = 50 
        y2'  =  - 1000 y1  +  y2  +  sin ( y1 - y2 ) ( 1 + e - x/2 ) ,
        y2(0) = 50 

        0 ≤ x ≤ 5 

Приводятся подпрограмма вычисления функций

        U1 (x, y)  =  sin ( y1 + y2 ) ( 1 + e - x/2 ) ,
        U2 (x, y)  =  sin ( y1 - y2 ) ( 1 + e - x/2 )

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

Unit TDE03R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FUDE03R_p, DE03R_p;

function TDE03R: String;

implementation

function TDE03R: String;
var
M,IERR :Integer;
XN,HMIN,EPS,P,XK,H :Real;
A :Array [0..3] of Real;
Y :Array [0..1] of Real;
YN :Array [0..1] of Real;
R :Array [0..30] of Real;
begin
Result := '';  { результат функции }
M := 2;
XN := 0.0;
YN[0] := 5.E1;
YN[1] := 5.E1;
HMIN := 1.E-10;
EPS := 1.E-8;
P := 1.E2;
ХК := 5.0;
H := 0.01;
A[0] := -5.E2;
A[2] := 1.0;
A[1] := -1.E3;
A[3] := 1.0;
DE03R (FUDE03R,M,XN,YN,XK,A,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('TDE03R',Result);  { вывод результатов в файл TDE03R.res }
exit;
end;

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

procedure FUDE03R(X :Real; var Y :Array of Real; var U :Array of Real;
                M :Integer);
implementation

procedure FUDE03R(X :Real; var Y :Array of Real; var U :Array of Real;
                M :Integer);
var
T :Real;
begin
T := 1.0 + Exp(-X/2.0);
U[0] := T*Sin(Y[0] + Y[1]);
U[1] := T*Sin(Y[0] - Y[1]);
end;

end.
end.

 Результаты: 

             Y(1)                                   Y(2)                               H 
     -3.508894600484-04      -8.340807464992-02       1.250000000001-03 

      IERR = 0