Текст подпрограммы и версий ( Фортран )
de30r.zip , de30d.zip 		Тексты тестовых примеров ( Фортран )
tde30r.zip , tde30d.zip
Текст подпрограммы и версий ( Си )
de30r_c.zip , de30d_c.zip 		Тексты тестовых примеров ( Си )
tde30r_c.zip , tde30d_c.zip
Текст подпрограммы и версий ( Паскаль )
de30r_p.zip , de30e_p.zip 		Тексты тестовых примеров ( Паскаль )
tde30r_p.zip , tde30e_p.zip

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

Назначение

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

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

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

          Y ' = AY ,

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

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

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

(1)    Y(XK) = [EXP( AH )]N*YN ,  H = ( XK-XN ) / N ,   N > 0.

Для аппроксимации матричной экспоненты EXP(AH) используется частичная сумма ряда 

                                     7
             EXP(AH)  =   ∑   ( AH )k / k!  ,
                                   k=0

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

(2)             || A H || ≤ 1

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

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

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

    SUBROUTINE  DE30R (M, XN, YN, XK, N, A, Y, E, R, R1, R2, IERR) 
   Параметры
M - количество уравнений в системе (тип: целый);
XN, YN - начальные значения аргумента и решения; YN представляет одномерный массив длины M (тип XN и YN: вещественный);
XK - значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования). XK может быть больше, меньше или pавно XN (тип: вещественный);
N - целое число равных частей, на которые делится отрезок интегрирования; при обращении к подпрограмме параметр N выбирается таким образом, чтобы || A (XK - XN) / N || ≤ 1 (см. замечания по использованию);
A - вещественный двумерный массив размера M*M, содержащий матрицу системы уравнений;
Y - искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента XK; Y представляет одномерный массив длины M. B случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: вещественный);
            E -
         R1, R2   		вещественные двумерные рабочие массивы размеpа M*M;
R - вещественный одномерный рабочий массив длины M;
IERR - целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если значение параметра N меньше 1. B этом случае интегрирование системы прекращается; при этом значения параметров M, XN, YN, XK, N, А не изменяются. 
   Версии
DE30D - вычисление решения задачи Коши для линейной однородной системы обыкновенных дифференциальных уравнений первого порядка с постоянными коэффициентами экспоненциальным методом с повышенной точностью. При этом параметры XN, YN, XK, A, Y, E, R1, R2, R должны иметь тип DOUBLE PRECISON. 
   Вызываемые подпрограммы
UTDE10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE30R.
UTDE11 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE30D. 
   Замечания по использованию
 

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

Значение параметра N, задаваемое при обращении к подпрограмме таким образом, что || A H || ≤ 1, есть, по - существу, число шагов, которые необходимо было бы выполнить, чтобы получить численную устойчивость приближенного решения данной задачи, если ее решать каким - нибудь методом типа Рунге - Кутта.

При работе подпрограммы значения параметров M, XN, YN, XK, N сохраняются. На месте матрицы А помещается матрица AH. Если после работы подпрограммы нет необходимости иметь начальное значение решения YN, то параметры YN и Y при обращении к ней можно совместить.

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

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

           y1'  =  -20 y1  +  y2
           y2'  =  -y1 - 20 y2
           y3'  =  -21 y1 - 19 y2 ,      0 < x < 10  ,
              
          y1(0)  =  10 ,     y2(0)  =  0 ,     y3(0)  =  0.

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

          DIMENSION  A(3, 3), YN(3), Y(3), E(3, 3), R(3), R1(3, 3), R2(3, 3)
          DATA  A /-20., -1., -21., 1., -20., -19., 0., 0., 0./
          M = 3
          XN = 0
          YN(1) = 10.
          YN(2) = 0.
          YN(3) = 0.
          XK = 10
          N = 200
          CALL  DE30R (M, XN, YN, XK, N, A, Y, E, R, R1, R2, IERR)

Результаты:

          Y(1)  =   0 ,  Y(2)  =   0 ,  Y(3)  =  -10

          IERR  =  0