Текст подпрограммы и версий ( Фортран )
de37r.zip , de37d.zip 		Тексты тестовых примеров ( Фортран )
tde37r.zip , tde37d.zip
Текст подпрограммы и версий ( Си )
de37r_c.zip , de37d_c.zip 		Тексты тестовых примеров ( Си )
tde37r_c.zip , tde37d_c.zip
Текст подпрограммы и версий ( Паскаль )
de37r_p.zip , de37e_p.zip 		Тексты тестовых примеров ( Паскаль )
tde37r_p.zip , tde37e_p.zip

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

Назначение

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

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

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

        Y ' (X) = A(X) * Y(X) + φ(X)   ,
        Y = ( y1, ..., yM )   ,
        A(X) = ( ai j(X) ),        i, j = 1, ..., M   ,
        φ(X) = ( φ1(X), ..., φM(X) )   .

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

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

трехстадийным  А - устойчивым неявным методом Рунге - Кутта шестого порядка точности. Решение вычисляется в одной точке  ХК, которая является концом интервала интегрирования.

Bсе компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы  P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Абсолютная погрешность приближенного решения оценивается по правилу Рунге.

Современные численные методы решения обыкновенных дифференциальных уравнений. Ред. Дж.Холл и Дж.Уатт. "Мир", M., 1979.

Butcher J.C. Implicit Runge - Kutta processes. Math. Comp., 18, 50 - 64, 1964.

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

    SUBROUTINE  DE37R (FA, FI, M, XN, YN, XK, HMIN, EPS, P,
                                            H, Y, RAB, IR, IERR) 
   Параметры
FA - подпрограмма вычисления матрицы системы  A (X) в любой точке  X. Первый оператор подпрограммы должен иметь вид:

SUBROUTINE  FA (A, X, M)

Здесь:  A - двумерный массив размера  M * M, в который помещается матрица системы, вычисленная при значении аргумента  X (тип параметров  A, X: вещественный);
FI - подпрограмма вычисления неоднородности правой части системы  φ (X) в любой точке  X. Первый оператор подпрограммы должен иметь вид:

SUBROUTINE  FI (G, X, M)

Здесь  G - одномерный массив длины  M, в который помещается неоднородность правой части системы, вычисленная при значении аргумента  X (тип параметров  G, X: вещественный);
M - количество уравнений в системе (тип: целый);
XN, YN - начальные значения аргумента и решения; в случае системы уравнений (т.е.  M ≠ 1)  YN представляет одномерный массив длины  M (тип: вещественный);
XK - значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования);  XK может быть больше, меньше или pавно  XN (тип: вещественный);
HMIN - минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
EPS - допустимая меpа погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный);
P - граница перехода, используемая при оценке погрешности решения (тип: вещественный);
H - вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если  XN < XK, отрицательным, если   XN > XK, или без всякого учета в виде абсолютной величины;
Y - искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента  XK; для системы уравнений (когда  M ≠ 1) задается одномерным массивом длины  M. В случае совпадения значений параметров  XN и  XK значение  Y полагается равным начальному значению  YN (тип: вещественный);
RAB - одномерный рабочий массив вещественного типа длины (10*M*M + 9*M + 1);
IR - целый одномерный рабочий массив длины 3*M;
IERR - целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая-нибудь компонента решения не может быть вычислена с требуемой точностью  EPS; в этом случае интегрирование системы прекращается; при желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров  HMIN и  H. 
   Версии
DE37D - вычисление решения задачи Коши для жесткой линейной системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования неявным методом Рунге - Кутта с удвоенным числом значащих цифр. При этом параметры  XN, YN, XK, HMIN, EPS, P, H, Y, RAB и параметры  A, G, X в подпрограммах  FA и  FI должны иметь тип DOUBLE PRECISION. 
   Вызываемые подпрограммы
       DE36R -
       DE36D   		выполнение одного шага численного интегрирования жесткой линейной системы обыкновенных дифференциальных уравнений первого порядка неявным методом Рунге - Kутта.
      UTDE16 -
      UTDE17   		подпрограммы выдачи диагностических сообщений.
  Подпрограммы DE36R, UTDE16 вызываются при работе подпрограммы DE37R, а подпрограммы DE36D, UTDE17 - при pаботе DE37D. 
   Замечания по использованию
 

B общем случае заданая точность не гарантируется.

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

При работе подпрограмм  FA и  FI значения параметров  X и M не должны изменяться.

Если после работы подпрограммы нет необходимости иметь начальное значение решения  YN, то параметры  YN и  Y при обращении к ней можно совместить. При этом следует иметь в виду, что в случае аварийного выхода из подпрограммы, т.е. со значением  IERR = 65, значение параметра  YN будет испорчено.

Tак как при интегрировании уравнений с помощью подпрограмм DE37R и DE37D используются общие блоки с именами COM36R и COM36D, соответственно, то пользователю не рекомендуется использовать для своих целей общие блоки с указанными именами.

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

      y1' = - 20y1 + y2   ,   y1 = 2   ,
      y2' =  19y1 - 2y2   ,   y2 = 18   ,   0 ≤ x ≤ 5   .

Точное решение системы:

    y1 = e - x + e - 21x   ,
    y2 = 19e - x - e - 21x   .

      SUBROUTINE  FBUT (A, X, M)
      DIMENSION  A(2, 2)
      A(1, 1) = - 20.
      A(1, 2) = 1.
      A(2, 1) = 19.
      A(2, 2) = - 2.
      RETURN
      END

      SUBROUTINE  FBUTFI (FI, X, M)
      DIMENSION  FI(2)
      FI(1) = 0.
      FI(2) = 0.
      RETURN
      END

      DIMENSION  YN(2), Y(2), RAB(59), IR(6)
      EXTERNAL  FBUT, FBUTFI
      M = 2
      XN = 0.
      YN(1) = 2.
      YN(2) = 18.
      XK = 5.
      HMIN = 1.E - 12
      EPS = 1.E - 5
      P = 100.
      DO 5 I = 1, 2
      H = 0.01
      CALL  DE37R (FBUT, FBUTFI, M, XN, YN, XK, HMIN, EPS,
     *                         P, H, Y, RAB, IR, IERR)
C      BЫЧИCЛEHИE TOЧHЫX ЗHAЧEHИЙ PEШEHИЯ
      Y1 = EXP(- XK) + EXP (- 21. * XK)
      Y2 = 19. * EXP(- XK) - EXP(- 21. * XK)
      PRINT 1, XK, Y, Y1, Y2, H
      EPS = 1.E - 7
   5 CONTINUE
      
Результаты:

после первого обращения к подпрограмме -

                     XK                            Y(1)                               Y(2)
      5.000000000000 + 00    6.737934370626 - 03    1.280207530417 - 01

                     H                               Y1                                 Y2
      1.280000000001 + 00    6.737946999117 - 03    1.280209929830 - 01

после второго обращения к подпрограмме -

                     XK                            Y(1)                                Y(2)
      5.000000000000 + 00    6.737946721380 - 03    1.280209877059 - 01

                     H                               Y1                                  Y2
      6.400000000003 - 01    6.737946999117 - 03    1.280209929830 - 01