Текст подпрограммы и версий ( Фортран )
de31r.zip , de31d.zip 		Тексты тестовых примеров ( Фортран )
tde31r.zip , tde31d.zip
Текст подпрограммы и версий ( Си )
de31r_c.zip , de31d_c.zip 		Тексты тестовых примеров ( Си )
tde31r_c.zip , tde31d_c.zip
Текст подпрограммы и версий ( Паскаль )
de31r_p.zip , de31e_p.zip 		Тексты тестовых примеров ( Паскаль )
tde31r_p.zip , tde31e_p.zip

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

Назначение

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

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

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

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

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

    SUBROUTINE  DE31R (FA, M, XN, YN, XK, N, NS, Y, A, E, 
                                            R, R1, R2, R3, IERR) 
   Параметры
FA - подпрограмма вычисления матрицы A(x) в любой точке X интервала интегрирования. Первый оператор подпрограммы должен иметь вид:
SUBROUTINE  FA(A, X, M).
Здесь 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.
  В каждом из этих случаев интегрирование системы прекращается. 
   Версии
DE31D - вычисление решения задачи Коши для линейной однородной системы обыкновенных дифференциальных уравнений первого порядка с переменными коэффициентами экспоненциальным методом с повышенной точностью. При этом параметры XN, YN, XK, Y, A, E, R, R1, R2, R3, а также параметры A, X в подпрограмме FA должны иметь тип DOUBLE PRECISION. 
   Вызываемые подпрограммы
UTDE10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE31R.
UTDE11 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE31D. 
   Замечания по использованию
 

Значение параметра 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  .

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

          SUBROUTINE  FA (A, X, M)
          DIMENSION   A(M, M)
          A(1, 1) = - (2. + X)/ (1. + X)
          A(M, M) = A(1, 1)
          A(M, 1) = X* (-20.0)
          A(1, M) = -A(M, 1)
          RETURN
          END
   
          REAL YN(2), Y(2), A(2, 2), E(2, 2), R(2), R1(2, 2), R2(2, 2), R3(2)
          EXTERNAL  FA
          M = 2
          XN = 0.
          YN(1) = 0.
          YN(2) = 1.
          XK = 6.
          N = 128
          IS1 = 128
          CALL  DE31R (FA, M, XN, YN, XK, N, IS1, Y, A, E, R, R1, R2, R3,
         *                        IERR)

Результаты:

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

          IERR  =  0