Подпрограмма:  DE52R (модуль DE52R_p)

Назначение

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

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

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

(1)          Y ' = A (X) Y + f (X)   ,

где  A (X) - квадратная матрица размера  M * M   ,
        f (X)  - М - мерная   вектор - функция   ,
с линейными краевыми условиями
(2)          BY (XN) = b   ,
(3)          CY (XK) = c   ,

где  B - прямоугольная матрица размера  (M - K) * M 
             (ранг pавен  M - K)   ,
      С - прямоугольная  матрица размера  K * M
             (ранг pавен  K)   ,
      b - (M - K) - мерный   вектоp   ,
      c - К - мерный   вектоp   , 

- методом ортогональной прогонки Годунова [1]. Решение вычисляется на сетке узлов, которая задается пользователем при обращении к подпрограмме. Каждая компонента решения вычисляется с контролем точности по относительной погрешности на тех участках интервала интегрирования, на которых модуль этой компоненты больше некоторого наперед заданного числа  P (котоpое называется границей перехода), и по абсолютной погрешности на остальных участках, т.е. там, где модуль проверяемой на точность компоненты меньше этого числа.

Реализованный в подпрограмме метод включает в себя в качестве подалгоритмов следующие задачи:

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

procedure DE52R(F :Proc_F_DE; M :Integer; XN :Real; NX :Integer;
                var X :Array of Real; HMIN :Real; EPS :Real;
                P :Real; NXORT :Integer; var XORT :Array of Real;
                K :Integer; BULODN :Boolean; var B :Array of Real;
                var C :Array of Real; var Y :Array of Real;
                var U :Array of Real; var IRAB :Array of Integer;
                var YR :Array of Real; var RK :Array of Real;
                var T :Array of Real; var IERR :Integer);

   Параметры 

   Версии 

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

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

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

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

(4)       U '''' - 24 * U ''' - 169 * U '' - 324 * U ' - 180 * U = 0 

Его частное решение
       U (x) = e-x - 2e-2x + e-3x
удовлетворяет начальным условиям
       U (0) = 0,   U ' (0) = 0,   U '' (0) = 2,   U ''' (0) = - 12.

При  x = 1 это решение удовлетворяет еще условиям:
       U (1) = 0.146996,    U ' (1) = 0.0241005.  

Находились численные значения  U (x) как решения уравнения (4) при следующих краевых условиях:

      U (0) = 0 ,   U ' (0) = 0  и  U (1) = 0.146996 ,   U ' (1) = 0.0241005 . 

Приводятся подпрограмма вычисления значений правой части соответствующей системы четырех дифференциальных уравнений первого порядка (к которой сводится уравнение (4) четвертого порядка) и фрагмент вызывающей программы. Так как эта система однородная, то в подпрограмме  FDE52R элемент BUL1 глобальной записи отсутствует.

Unit tde52r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE52R_p, DE52R_p;

function tde52r: String;

implementation

function tde52r: String;
var
K,M,NXS,NXO,I,_i,IERR :Integer;
XN,HMIN,EPS,P,H :Real;
BUL :Boolean;
IRAB :Array [0..3] of Integer;
XS :Array [0..1] of Real;
XO :Array [0..19] of Real;
U :Array [0..15] of Real;
YR :Array [0..11] of Real;
RK :Array [0..15] of Real;
T :Array [0..5] of Real;
Y :Array [0..7] of Real;
C :Array [0..9] of Real;
B :Array [0..9] of Real;
const
C0 :Array [0..9] of Real = ( 1.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.146996,
0.0241005 );
B0 :Array [0..9] of Real = ( 1.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0 );
label
_1;
begin
for I:=0 to 9 do
 begin
 C[I] := C0[I];
 B[I] := B0[I];
 end;
Result := '';  { результат функции }
K := 2;
M := 4;
XN := 0.0;
NXS := 2;
XS[0] := 1.0;
XS[1] := 0.0;
HMIN := 1.E-10;
EPS := 1.E-6;
P := 0.1;
NXO := 10;
H := 0.1;
XO[0] := H;
for I:=2 to NXO do
 begin
_1:
  XO[I-1] := XO[I-2]+H;
 end;
BUL := True;
DE52R(FDE52R,M,XN,NXS,XS,HMIN,EPS,P,NXO,XO,K,BUL,B,
     C,Y,U,IRAB,YR,RK,T,IERR);
Result := Result + #$0D#$0A;
Result := Result + '  XS' + #$0D#$0A + #$0D#$0A;
for _i:=0 to 1 do
 begin
  Result := Result + Format('%20.16f ',[XS[_i]]);
  if ( ((_i+1) mod 2)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + '  Y' + #$0D#$0A + #$0D#$0A;
for _i:=0 to 7 do
 begin
  Result := Result + Format('%20.16f ',[Y[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A + '  IERR=' + Format('%2d ',[IERR]);
UtRes('tde52r',Result);  { вывод результатов в файл tde52r.res }
exit;
end;

end.

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

procedure fde52r(X :Real; var Y :Array of Real; var Z :Array of Real;
                M :Integer);

implementation

procedure fde52r(X :Real; var Y :Array of Real; var Z :Array of Real;
                M :Integer);
begin
Z[0] := Y[1];
Z[1] := Y[2];
Z[2] := Y[3];
Z[3] := 24.0*Y[3]+169.0*Y[2]+324.0*Y[1]+180.0*Y[0];
end;

end.

Результаты:

       Y(1, 1)  =     0,1469960000000000
       Y(1, 2)  =    -0,0000002357020366
       Y(2, 1)  =     0,0241005000000000
       Y(2, 2)  =     0,0000009105752631
       Y(3, 1)  =    -0,2667186437063472
       Y(3, 2)  =     1,9999974802225444
       Y(4, 1)  =     0,4532534661842457
       Y(4, 2)  =   -11,9999934834509414

       IERR  =  0