Текст подпрограммы и версий
de31r_p.zip , de31e_p.zip 		Тексты тестовых примеров
tde31r_p.zip , tde31e_p.zip

Подпрограмма:  DE31R (модуль DE31R_p)

Назначение

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

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

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

(1)          Y ' = A(x) Y

здесь Y = ( y1, ..., yM ), А(x) - квадратная вещественная матрица порядка M, элементами которой являются трижды непрерывно дифференцируемые функции. Начальные значения заданы в точке XN : 

          Y(XN) = YN ,     YN = ( y10,..., yM0 ) ,

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

Интервал интегрирования разбивается на NS равных частей длины T, на каждой из которых исходная система уравнений (1) аппроксимируется системой уравнений с постоянными коэффициентами 

               ~         ~   ~         ~
(2)          Y '  =  A  Y ,      A = CONST.

Решение системы (2) вычисляется с помощью матричной экспоненты 

                         ~
(3)          EXP( A T )  =
                                               ~           N
                             =  [ EXP(  A  H ) ]      ,   H = T/N ,   N > 0

Значения N и NS задаются пользователем при обращении к подпрограмме, причем значение N выбирается таким образом, чтобы выполнялось условие 

(4)               max        || A(x) H ||  ≤  1
           XN ≤ x ≤ XK

(в качестве || A(x) H || можно взять максимальную сумму модулей элементов матрицы A (x) H по стpокам).

С.Ф.Залеткин, O численном решении задачи Коши для обыкновенных линейных однородных дифференциальных уравнений на больших отрезках интегрирования, "Вычислительные методы и программирование", вып.26, Изд-во МГУ, 1977.

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

procedure DE31R(FA :Proc_F3_DE; M :Integer; XN :Real;
                var YN :Array of Real; XK :Real; N :Integer;
                NS :Integer; var Y :Array of Real;
                var A :Array of Real; var E :Array of Real;
                var R :Array of Real; var R1 :Array of Real;
                var R2 :Array of Real; var R3 :Array of Real;
                var IERR :Integer);


   Параметры
FA - подпрограмма вычисления матрицы A(x) в любой точке X интервала интегрирования. Первый оператор подпрограммы должен иметь вид:
procedure FA (var A :Array of Real; X :Real; M :Integer);
Здесь A - двумерный массив размера M*M, в который помещается матрица системы, вычисленная при значении аргумента X, т.е. A (x) (тип параметров A и X: вещественный);
M - количество уравнений в системе (тип: целый);
XN, YN - начальные значения аргумента и решения; YN представляет одномерный массив длины M (тип: вещественный);
XK - значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования). XK может быть больше, меньше или pавно XN (тип: вещественный);
N - целое число равных частей, на которые делится отрезок интегрирования системы линейных дифференциальных уравнений с постоянными коэффициентами; при обращении к подпрограмме параметр N выбирается таким образом, чтобы
|| A(x)*( ((XK - XN)/NS)/N )|| ≤ 1 для XN ≤ x ≤ XK (см. замечания по использованию);
NS - число равных частей ,на которые разбивается интервал интегрирования. NS должно быть таким, чтобы система уравнений (2) достаточно хорошо аппроксимировала исходную систему уравнений с переменными коэффициентами; NS > 0 (тип: целый);
Y - искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента XK; Y представляет одномерный массив длины M. В случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: вещественный);
            A, E -
         R1, R2   		вещественные двумерные рабочие массивы размеpа M*M;
R, R3 - вещественные одномерные рабочие массивы длины M;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом:
IERR=65 - когда значение параметра N меньше 1;
IERR=66 - когда значение параметра NS меньше 1.
  В каждом из этих случаев интегрирование системы прекращается. 
   Версии
DE31E - вычисление решения задачи Коши для линейной однородной системы обыкновенных дифференциальных уравнений первого порядка с переменными коэффициентами экспоненциальным методом с расширенной (Extended) точностью. При этом параметры XN, YN, XK, Y, A, E, R, R1, R2, R3, а также параметры A, X в подпрограмме FA должны иметь тип Extended. 
   Вызываемые подпрограммы
UTDE10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE31R.
UTDE11 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE31E. 
   Замечания по использованию
 

Значение параметра N, задаваемое при обращении к подпрограмме таким образом, чтобы
max  || A(x)*(((XK-XN)/NS)/N) || ≤ 1 ,  XN ≤ x ≤ XK есть, по - существу, число шагов, котоpое необходимо было бы выполнить, чтобы получить численную устойчивость приближенного решения каждой линейной системы с постоянными коэффициентами, если ее интегрировать каким - нибудь методом типа Рунге - Кутта.

При работе подпрограммы значения параметров M, XN, YN, XK, N, NS сохраняются.

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

Использование подпрограммы иллюстрируется на примере: 

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

          y1(0)  =  0 ,     y2(0)  =  1 ,     0 ≤ x ≤ 6  .

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

Unit TDE31R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FADE31R_p, DE31R_p;

function TDE31R: String;

implementation

function TDE31R: String;
var
M,N,IS1,_i,IERR :Integer;
XN,ХК :Real;
YN :Array [0..1] of Real;
YK :Array [0..1] of Real;
A :Array [0..3] of Real;
E :Array [0..3] of Real;
R :Array [0..1] of Real;
R1 :Array [0..3] of Real;
R2 :Array [0..3] of Real;
R3 :Array [0..1] of Real;
Z :Array [0..1] of Real;
Y :Array [0..1] of Real;
begin
Result := '';  { результат функции }
M := 2;
XN := 0.0;
YN[0] := 0.0;
YN[1] := 1.0;
ХК := 6.0;
N := 128;
IS1 := 128;
DE31R(FADE31R,M,XN,YN,XK,N,IS1,Y,A,E,R,R1,R2,R3,IERR);
Result := Result + Format('%s',[' XK=']);
Result := Result + Format('%20.16f ',[XK]);
Result := Result + Format('%s',['      Y=']);
Result := Result + #$0D#$0A;
for _i:=0 to 1 do
 begin
  Result := Result + Format('%20.16f ',[Y[_i]]);
  if ( ((_i+1) mod 2)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TDE31R',Result);  { вывод результатов в файл TDE31R.res }
exit;
end;

end.

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

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

implementation

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

end.

Результаты:

          Y(1)  =   0.3395896401 * 10-3
          Y(2)  =  -0.1004661347 * 10-3

          IERR  =  0