Текст подпрограммы и версий
de02r_p.zip , de02e_p.zip
Тексты тестовых примеров
tde02r_p.zip , tde02e_p.zip

Подпрограмма:  DE02R (модуль DE02R_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) ) , 

где 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(1),  y(2) , невелика. Также предполагается, что нелинейный член U (x, y) является достаточно малым. По заданному значению решения YX в узле  xn  вычисляется значение этого решения в узле  xn + H . Вычисление производится по методу Лоусона.

Метод Лоусона заключается в следующем. Исходная система уравнений с помощью замены искомой функции  y (x)  на [ xn , xn + 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) .

Значение H может быть меньше или равно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.

Все компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине меньше некоторой наперед заданной константы 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 DE02R(FU :Proc_F_DE; M :Integer; var JSTART :Integer;
                var A :Array of Real; 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 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);

Параметры

FU - подпрограмма вычисления функции U (x, y) в правой части системы. Первый оператор подпрограммы должен иметь вид:
procedure FU (X :Real; var Y :Array of Real; var U :Array of Real;M :Integer);
Здесь Y и U - одномерные массивы длины M. В массив U помещается значение функции U (x, y), вычисленное при значении аргументов X и Y (тип параметров X, Y, U: вещественный);
M - количество уравнений в системе (тип: целый);
JSTART - целый указатель режима использования подпрограммы, имеющий следующие значения:
0 - первое обращение к подпрограмме должно быть выполнено с нулевым значением JSTART;
+1 - выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных и шага интегрирования, заданных параметрами X, YX и H соответственно; при этом само значение величины шага интегрирования H, на который предполагается продолжить решение системы, должно совпадать с тем значением шага, который был фактически выполнен подпрограммой при предыдущем обращении к ней. Если же требуется выполнить шаг интегрирования, отличный от реально сделанного подпрограммой во время предыдущего обращения к ней, то для этого необходимо обратиться к подпрограмме со значением JSTART = 0;
+2 - то же, что и для JSTART = + 1, но только с той разницей, что величина шага интегрирования, на который предполагается продолжить решение, в два раза больше того значения шага, который был фактически выполнен при предыдущем обращении к подпрограмме;
-1 - повторить последний шаг интегрирования с новыми значениями параметров H и/или HMIN;
  при выходе из подпрограммы параметр JSTART полагается равным 2, если выполнив заданный в H шаг, подпрограмма рекомендовала использовать для интегрирования вдвое большее его значение, и 1 в противном случае, т.е. в том случае, когда рекомендованное значение шага равно только что выполненному шагу; рекомендованное значение шага на выходе из подпрограммы запоминается в переменной H;
A - вещественный двумерный массив размера M*M, содержащий элементы матрицы A;
HMIN - минимальное значение абсолютной величины шага, которое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
EPS - допустимая мера погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный);
P - граница перехода, используемая при оценке меры погрешности решения (тип: вещественный);
YX, X - заданные вещественные значения решения и соответствующее ему значение аргумента; в результате работы подпрограммы в X получается новое значение аргумента, в YX - соответствующее значение решения; в случае системы уравнений т.е. когда M ≠ 1, YX задается одномерным массивом длиной M;
H - вещественная переменная, содержащая значение шага интегрирования; если для этого значения шага точность приближенного решения достигается, то именно он и реализуется подпрограммой, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность EPS; на выходе из подпрограммы H содержит рекомендуемое подпрограммой значение следующего шага интегрирования, определяемое ею с целью достижения более экономного способа интегрирования;
BUL - логическая переменная, значение которой при обращении к подпрограмме полагается равным TRUE, если заданный в H шаг выводит в конец интервала интегрирования, и FALSE в противном случае; в результате работы подпрограммы BUL равно FALSE, если вместо исходного шага интегрирования был реализован меньший шаг; в противном случае, т.е. когда был выполнен именно заданный при обращении в H шаг, значение параметра BUL не меняется;
XP, YP - вещественная рабочая переменная и одномерный рабочий массив длины M соответственно; значения параметров XP, YP на выходе из подпрограммы равны тем значениям, которые имели параметры X, YX при входе в нее (т.е. предыдущий узел и решение в нем);
       E1, E2 -
       E3, E4  
вещественные двумерные рабочие массивы размера M*M;
     R, R0, R1 -
     R2, R3, YD  
вещественные одномерные рабочие массивы длины M;
IERR - целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров H, HMIN и значением JSTART = - 1 .

Версии

DE02E - выполнение одного шага численного интегрирования квазилинейной системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона с удвоенным числом значащих цифр. При этом параметры HMIN, EPS, P, A, YX, X, H, XP, YP, E1, E2, E3, E4, R, R0, R1, R2, R3, YD и параметры X, Y, DY в подпрограмме FU должны иметь тип Extended.

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

AME2R - подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы DE02R.
AME2E - подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы DE02E.
UTDE20 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE02R;
UTDE21 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE02E.
  Кроме того, при работе подпрограмм DE02R и DE02E вызываются рабочие подпрограммы DE02RS, DE04RQ и DE02ES, DE04EQ соответственно.

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

 

Данная подпрограмма предназначена для интегрирования квазилинейных систем, имеющих малый нелинейный член U (x, y).

В общем случае заданная точность не гарантируется.

При работе подпрограммы значения параметров M, HMIN, EPS, P, A сохраняются.

Между последовательными обращениями к подпрограмме со значениями параметра JSTART = 1, 2 пользователь не должен изменять содержимое массивов A, E1, E2.

При работе подпрограммы FU значения параметров X, Y и M не должны изменяться.

При обращении к подпрограмме со значением JSTART = - 1 в качестве исходных значений аргумента и решения принимаются значения параметров XP и YP соответственно, т.е. те значения, которые эти параметры получили после самого последнего обращения к подпрограмме с неотрицательным значением JSTART.

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

Подпрограммы DE02R и DE02E предназначены также для интегрирования жестких дифференциальных уравнений (1).

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

        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 
     
        x ≥ 0  

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

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

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

Unit TDE02R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FUDE02R_p, DE02R_p;

function TDE02R: String;

implementation

function TDE02R: String;
var
M,JSTART,IH,IERR :Integer;
X,HMIN,EPS,P,H,ХР :Real;
BUL :Boolean;
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;
YP :Array [0..1] of Real;
YX :Array [0..1] of Real;
label
_100,_101,_102,_103,_104,_105,_106;
begin
Result := '';  { результат функции }
M := 2;
X := 0.0;
YX[0] := 5.E1;
YX[1] := 5.E1;
HMIN := 1.E-10;
EPS := 1.E-8;
P := 1.E2;
JSTART := 0;
H := 0.01;
A[0] := -5.E2;
A[2] := 1.0;
A[1] := -1.E3;
A[3] := 1.0;
IH := 0;
_100:
IH := IH + 1;
DE02R (FUDE02R,M,JSTART,A,HMIN,EPS,P,YX,X,H,BUL,XP,YP,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;
goto _100;
_102:
EPS := 0.1;
goto _100;
_103:
JSTART := -1;
EPS := 1.E-10;
goto _100;
_104:
EPS := 1.E-4;
H := -1.E-5;
JSTART := -1;
goto _100;
_105:
H := -1.E-4;
EPS := 1.E-8;
goto _100;
_106:
UtRes('TDE02R',Result);  { вывод результатов в файл TDE02R.res }
exit;
end;

end.

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

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

procedure FUDE02R(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.

 Результаты: 
 
 После первого обращения к подпрограмме
              X                                    YX(1)                                  YX(2) 
      1.953125000001 - 05       4.951502369251 + 01       4.902915763238 + 01 
              H                                    R1(1)                                   R1(2) 
      1.953125000001 - 05       3.523503740624 - 09        -1.435788969199 - 10 

 После второго обращения к подпрограмме
              X                                    YX(1)                                  YX(2) 
      5.004882812498 - 05       4.876668035128 + 01       4.753104564862 + 01 
              H                                    R1(1)                                   R1(2) 
      3.051757812500 - 05       5.044663945828 - 10         1.331015179553 - 09 

 После третьего обращения к подпрограмме
              X                                    YX(1)                                  YX(2) 
      8.056640624998 - 05       4.802966483223 + 01       4.605556652445 + 01 
              H                                    R1(1)                                   R1(2) 
      6.103515625000 - 05      -2.398155629634 - 08        -9.972912569839 - 10 

 После четвертого обращения к подпрограмме
              X                                    YX(1)                                  YX(2) 
      5.767822265623 - 05       4.858137986262 + 01       4.716006872657 + 01 
              H                                    R1(1)                                   R1(2) 
      7.629394531250 - 06      -4.268561800317 - 11        -3.880510727560 - 12 

 После пятого обращения к подпрограмме
              X                                    YX(1)                                  YX(2) 
      4.004882812492 - 05       4.901063239173 + 01       4.801943553990 + 01 
              H                                    R1(1)                                   R1(2) 
     -2.000000000005 - 05       8.537123600634 - 11         7.761021455121 - 12 

 После шестого обращения к подпрограмме
              X                                    YX(1)                                  YX(2) 
     -5.995117187529 - 05       4.950228966336 + 01       4.900345166435 + 01 
              H                                    R1(1)                                   R1(2) 
     -1.000000000002 - 04      -7.190586378176 - 09        -8.769954244295 - 10