Текст подпрограммы и версий ( Фортран )
de07r.zip , de07d.zip 		Тексты тестовых примеров ( Фортран )
tde07r.zip , tde07d.zip
Текст подпрограммы и версий ( Си )
de07r_c.zip , de07d_c.zip 		Тексты тестовых примеров ( Си )
tde07r_c.zip , tde07d_c.zip
Текст подпрограммы и версий ( Паскаль )
de07r_p.zip , de07e_p.zip 		Тексты тестовых примеров ( Паскаль )
tde07r_p.zip , tde07e_p.zip

Подпрограмма:  DE07R

Назначение

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

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

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

                           Y ' (x)   =   A * Y(x) + φ(x)  ,

Здесь    Y   =   ( y1(x), ... , yM(x) )  ,
             A   -   квадратная вещественная матрица порядка M ,
            φ(x)   =   ( α1 e b1, ... , αM e bM )  ,
где         bi    =   βi x  ,     i = 1, ... , M   .

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

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

Решение вычисляется в заданной точке XK = XN + T по следующей формуле 

     Y(XK)  =  Y(XN+T)  =  e AT Y(XN)  +
                                                                    M
                                                                +  ∑   αm e  bm  g m ( T )   , 
                                                                   m=1
 где    bm  =  βm XK  , 
        g m(T)  -  m-ый столбец матрицы
                                 T
            Φm (T)   =    ∫  e λm dt  ,        λm  =  ( A - βm T )( T - t )
                               0

Матричная экспонента  e AT и интегралы от матричной экспоненты  g m (T) вычисляются с помощью рекуррентных соотношений, которым удовлетворяют матрица  e AH и векторы  g m (H) для значений скалярного аргумента H = H j , где 

        H j  =  ( T / K )K  ,   K  =  2 j  ,      1 ≤ j ≤ N  .

Значение N выбирается таким, чтобы выполнялось неpавенство 

                    || A  HN || ≤ 1

(в качестве || A HN || можно взять максимальную сумму модулей элементов матрицы A HN по стpокам). Это значение N задается пользователем при обращении к подпрограмме.

Залеткин С.Ф. Численное решение системы линейных обыкновенных дифференциальных уравнений с постоянными коэффициентами на больших промежутках интегрирования. - В сб.: Вопросы конструирования библиотек программ. M.: Изд - во МГУ, 1985.

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

    SUBROUTINE  DE07R (M, XN, YN, XK, N, Q, A, ALP, BET, JSTART, 
                                             Y, E, G, A1, R, R1, IERR) 
   Параметры
M - количество уравнений в системе (тип: целый);
XN, YN - начальные значения аргумента и решения; YN представляет собой одномерный массив длины M (тип XN и YN: вещественный);
XK - значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования); XK может быть больше, меньше или pавно XN (тип: вещественный);
N - показатель степени 2N, равный числу частей, на которые делится отрезок интегрирования [ XN, XK ]; при обращении к подпрограмме параметp N задается таким образом, чтобы || A (XK - XN) / 2N  || ≤ 1 (тип: целый);
Q - число членов ряда для аппроксимации матричной экспоненты 
                   Q
   eAH   =   ∑   ( AH ) i / i !  ,    H = ( XK - XN ) / 2N   ;
                  i=0
задается при обращении к подпрограмме (тип: целый);
A - вещественный двумерный массив порядка M на M, содержащий матрицу коэффициентов системы уpавнений;
ALP - вещественный одномерный массив длины M, содержащий коэффициенты  α i : ALP (I) = α i;
BET - вещественный одномерный массив длины M, содержащий показатели степени  β i ; BET (I) = β i;
JSTART - целый показатель режима использования подпрограммы, принимающий перед обращением к подпрограмме значение 0; в том случае, когда требуется вычислить решение задачи Коши на равномерной сетке узлов  xj = XN + j T,  j = 1, 2, ..., то первое обращение к подпрограмме для вычисления решения в первом узле сетки  x1 = XN + T производится со значением JSTART = 0, а для получения решения в каждом следующем узле  x j ,  j ≥ 2, рекомендуется повтоpно обращаться к подпрограмме со значением JSTART = 1, при этом необходимо положить XN = xj - 1, YN = y ( xj - 1 ) (см. также замечания по использованию);
Y - вещественный одномерный массив длины M, содержащий вычисленное решение задачи в точке XK;
            E, G -
            A1   		вещественные двумерные рабочие массивы размеpа M на M;
R, R1 - вещественные одномерные рабочие массивы длины M;
IERR - целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если N < 1, и 66, если Q < 1. В этом случае решение задачи прекращается; при этом значения параметров M, XN, YN, XK, N, Q, A, ALP, BET не изменяются. 
   Версии
DE07D - то же, что и DE07R, но решение задачи вычисляется с удвоенной значностью. При этом параметры XN, YN, XK, A, ALP, BET, Y, E, G, A1, R, R1 должны иметь тип DOUBLE PRECISION. 
   Вызываемые подпрограммы
UTDE18 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE07R.
UTDE19 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE07D. 
   Замечания по использованию
 

Значения параметров M, XN, YN, XK, N, Q, A, ALP, BET при работе подпрограмм сохраняются.

При многократном использовании подпрограммы для вычисления решения задачи Коши на равномерной сетке узлов со значением JSTART = 1 значения рабочих массивов, задаваемых параметрами A1 и G, не должны изменяться в вызывающей программе между последовательными обращениями к подпрограмме.

Если матрица A ( XK - XN ) имеет собственные значения с положительными вещественными частями, то работа программы при вычислении матричной экспоненты  eA (XK - XN) при больших значениях | XK - XN | может привести к переполнению арифметического устройства.

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

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

                                     
        y1'  =  - y1  +  2 y2  - 2 e - k  ,     k = 2x  ,     y1(0) = 5 
        y2'  =  - 3 y1  -  2 y2  +  15 e - x ,                  y2(0) = 1.

Решение вычисляется в двух точках  x1 = 5,  x2 = 10 
(точное решение задачи  y1 = 5 e- x,  y2 = e - k,  k = 2 x). 

      PROGRAM T07R
      REAL YN(2), A(2, 2), ALP(2), BET(2), Y(2), E(2, 2), G(2, 2), 
     *          A1(2, 2), R(2), R1(2)
      INTEGER Q
      M = 2
      XN = 0.
      YN(1) = 5.
      YN(2) = 1.
      XK = 5.
      N = 7
      Q = 10
      A(1, 1) = -1.
      A(1, 2) = 2.
      A(2, 1) = -3.
      A(2, 2) = -2.
      ALP(1) = -2.
      ALP(2) = 15.
      BET(1) = -2.
      BET(2) = -1.
      JSTART = 0
      CALL  DE07R (M, XN, YN, XK, N, Q, A, ALP, BET, JSTART, Y, E, G,
     *                         A1, R, R1, IERR)
      YT1 = 5.*EXP(-XK)
      YT2 = EXP(-2.*XK)
      YN(1) = Y(1)
      YN(2) = Y(2)
      XN = XK
      XK = XK + XK
      JSTART = 1
      CALL  DE07R (M, XN, YN, XK, N, Q, A, ALP, BET, JSTART, Y, E, G,
     *                         A1, R, R1, IERR)
      YT1 = 5.*EXP(-XK)
      YT2 = EXP(-2.*XK)
      END

Результаты:

      IERR = 0

      XK  =  5.00000000000 + 00          Y(1)  =  3.36897349947-02
      XK  =  5.00000000000 + 00          YT1  =  3.36897349956-02
                                                             Y(2)  =  4.53999302508-05
                                                             YT2  =  4.53999297625-05


      IERR = 0

      XK  =  1.00000000000 + 01          Y(1)  =  2.26999648806-04
      XK  =  1.00000000000 + 01          YT1  =  2.26999648813-04
                                                             Y(2)  =  2.06115688405-09
                                                             YT2  =  2.06115362244-09