Подпрограмма:  DE36R (модуль DE36R_p)

Назначение

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

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

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

        Y ' (X) = A(X) * Y(X) + φ(X)   ,
        Y = ( y1, ..., yM )   ,
        A(X) = ( a i j(X) ) ,        i, j = 1, ..., M   ,
        φ(X) = ( φ1(X), ..., φM(X) )   .

По заданному значению решения  YX в узле  X_n вычисляется значение этого решения в узле  X_n + H. Вычисление производится по неявной трехчленной формуле Рунге - Кутта шестого порядка точности. Значение  H может быть меньше или pавно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.

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

Современные численные методы решения обыкновенных дифференциальных уравнений. Ред. Дж.Холл и Дж.Уатт. "Мир", M., 1979.

Butcher J.C. Implicit Runge - Kutta processes. Math. Comp., 18, 50 - 64, 1964.

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

procedure DE36R(FA :Proc_F3_DE; FI :Proc_F3_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 IR :Array of Integer; var R1 :Array of Real;
                var R2 :Array of Real; var R3 :Array of Real;
                var R4 :Array of Real; var R5 :Array of Real;
                var IERR :Integer);

   Параметры 

   Версии 

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

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

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

      y1' = - 20y1 + y2   ,
      y2' = 19y1 - 2y2   ,

   Точное решение системы:

      y1 = e - x + e - 21x   ,
      y2 = 19e - x - e - 21x   .

Unit TDE36R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FB1DE36R_p, FB2DE36R_p, DE36R_p;

function TDE36R: String;

implementation

function TDE36R: String;
var
IH,M,JSTART,_i,IERR :Integer;
X,HMIN,EPS,P,H,Y1,Y2,ХР :Real;
BUL :Boolean;
YX :Array [0..1] of Real;
YP :Array [0..1] of Real;
IR :Array [0..5] of Integer;
R1 :Array [0..41] of Real;
R2 :Array [0..5] of Real;
R3 :Array [0..3] of Real;
R4 :Array [0..1] of Real;
R5 :Array [0..1] of Real;
label
_6,_104,_105,_106,_20;
begin
Result := '';  { результат функции }
IH := 0;
M := 2;
X := 0.0;
YX[0] := 2.0;
YX[1] := 18.0;
HMIN := 1.E-6;
EPS := 1.E-4;
P := 100.0;
JSTART := 0;
H := 0.01;
BUL := False;
_6:
DE36R(FB1DE36R,FB2DE36R,M,JSTART,HMIN,EPS,P,YX,X,H,BUL,XP,YP,
     IR,R1,R2,R3,R4,R5,IERR);
Result := Result + Format('%s',[' IERR=']);
Result := Result + Format('%4d ',[IERR]) + #$0D#$0A;
IH := IH+1;
Y1 := Exp(-X)+Exp(-21.0*X);
Y2 := 19.0*Exp(-X)-Exp(-21.0*X);
Result := Result + Format('%20.16f ',[X]) + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 1 do
 begin
  Result := Result + Format('%20.16f ',[YX[_i]]);
  if ( ((_i+1) mod 2)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%20.16f  %20.16f  %20.16f ',
 [H,Y1,Y2]) + #$0D#$0A;
case IH of
 1: goto _6;
 2: goto _6;
 3: goto _104;
 4: goto _105;
 5: goto _106;
 6: goto _6;
 7: goto _20;
end;
_104:
JSTART := -1;
H := 0.005;
goto _6;
_105:
JSTART := -1;
H := 0.02;
goto _6;
_106:
JSTART := -1;
H := 0.03;
goto _6;
_20:
UtRes('TDE36R',Result);  { вывод результатов в файл TDE36R.res }
exit;
end;

end.

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

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

implementation

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

end.

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

procedure fb2de36r(var FI :Array of Real; X :Real; M :Integer);

implementation

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

end.

Результаты:

после первого обращения к подпрограмме -

                      X                               YX(1)                              YX(2)
       9.999999999991 - 03    1.800634079711 + 00    1.800036259522 + 01
                      H                               Y1                                   Y2
       1.999999999998 - 02    1.800634079724 + 00    1.800036259528 + 01

после второго обращения к подпрограмме -

                      X                               YX(1)                              YX(2)
       2.999999999997 - 02    1.503037334354 + 00    1.790587333651 + 01
                      H                               Y1                                   Y2
       3.999999999996 - 02    1.503037334556 + 00    1.790587333636 + 01

после третьего обращения к подпрограмме -

                      X                               YX(1)                              YX(2)
       6.999999999994 - 02    1.162319294410 + 00    1.748555710353 + 01
                      H                               Y1                                   Y2
       7.999999999993 - 02    1.162319305093 + 00    1.748555709302 + 01

после четвертого обращения к подпрограмме -

                      X                               YX(1)                              YX(2)
       3.499999999997 - 02    1.445110875044 + 00    1.786699744989 + 01
                      H                               Y1                                   Y2
       9.999999999991 - 03    1.445110875235 + 00    1.786699744992 + 01

после пятого обращения к подпрограмме -

                      X                               YX(1)                              YX(2)
       4.999999999995 - 02    1.301167173347 + 00    1.772342131648 + 01
                      H                               Y1                                   Y2
       3.999999999996 - 02    1.301167173613 + 00    1.772342131642 + 01

после шестого обращения к подпрограмме -

                      X                               YX(1)                              YX(2)
       5.999999999995 - 02    1.225418558231 + 00    1.760987211324 + 01
                      H                               Y1                                   Y2
       5.999999999995 - 02    1.225418560085 + 00    1.760987211161 + 01

после седьмого обращения к подпрограмме -

                      X                               YX(1)                              YX(2)
       1.199999999999 - 01    9.673799790717 - 01    1.677102875494 + 01
                      H                               Y1                                   Y2
       1.199999999999 - 01    9.673800434666 - 01    1.677102869086 + 01