Текст подпрограммы и версий ( Фортран )
de05r.zip , de05d.zip 		Тексты тестовых примеров ( Фортран )
tde05r.zip , tde05d.zip
Текст подпрограммы и версий ( Си )
de05r_c.zip , de05d_c.zip 		Тексты тестовых примеров ( Си )
tde05r_c.zip , tde05d_c.zip
Текст подпрограммы и версий ( Паскаль )
de05r_p.zip , de05e_p.zip 		Тексты тестовых примеров ( Паскаль )
tde05r_p.zip , tde05e_p.zip

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

Назначение

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

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

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

(1)                         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   =   ( y1 0,..., yM 0 ) .

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

Метод Лоусона является одношаговым методом и заключается в следующем. Допустим, что искомое решение системы (1) уже вычислено в некоторой точке  x = xn  интервала интегрирования, т.е. известно  yn ≈ y (xn). Для отыскания решения Y(xn + 1) = Y(xn + H) в следующем узле  xn + 1 = xn + H выполняются такие действия. Исходная система уравнений с помощью замены искомой функции Y (x) на  xn ≤ x ≤ xn + H  по формуле 

                             Y(x)  =  exp [ ( x - xn ) A0 ] Z(x)  ,

где A0 - некоторая постоянная матрица, преобразуется в систему уравнений относительно новой неизвестной функции Z (X): 

(2)  Z ' (x)  =  A1(x) Z(x)  +  φ1(x)  =  exp [ - ( x - xn ) A0 ] { A(x) - A0 } *
                                  * exp [ ( x - xn ) A0 ] Z(x) + exp [ ( - ( x - xn ) A0 ] φ(x)

           xn ≤ x ≤ xn + H

Данное преобразование выполняется самой подпрограммой. В качестве матрицы A0 подпрограмма выбирает матрицу A0 = A (xn + H /2). Если шаг H достаточно мал, то преобразование позволяет уменьшить характеристические корни матрицы A1 (x) по сравнению с характеристическими корнями исходной матрицы A (x). Это приводит к уменьшению константы Липшица системы (2) по сравнению с константой Липшица системы (1). Для решения системы (2) применяются формулы классического метода Рунге - Кутта четвертого порядка точности, причем одновременно с решением (2) производится обратное преобразование от функции Z (x) к функции Y (x).

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

J.Douglas Louson. Generalized Runge - Kutta processes for stable systems with large Lipshitz constants, SIAM Journal on Numerical Analisys. Vol 4, No.3, 1967.

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

    SUBROUTINE  DE05R (FA, FI, M, XN, YN, XK, HMIN, EPS, P, H, Y,
                                             R, 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 - значение аргумента, при котором требуется вычислить решение задачи Коши (конец интервала интегрирования); XK может быть больше, меньше или равно XN (тип: вещественный);
HMIN - минимальное значение абсолютной величины шага, которое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
EPS - допустимая мера погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный);
P - граница перехода, используемая при оценке меры погрешности решения (тип: вещественный);
H - вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если XN < XK, отрицательным, если XN > XK, или без всякого учета в виде абсолютной величины; на выходе из подпрограммы содержит значение последнего шага интегрирования;
Y - искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента XK; для системы уравнений (когда M ≠ 1) задается одномерным массивом длины M. В случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: вещественный);
R - одномерный вещественный рабочий массив длины (5*M*M + 8*M + 1);
IERR - целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS; в этом случае интегрирование системы прекращается; при желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров HMIN и H. 
   Версии
DE05D - вычисление решения задачи Коши для линейной устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица в конце интервала интегрирования методом Лоусона с удвоенным числом значащих цифр. При этом параметры XN, YN, XK, HMIN, EPS, P, H, Y, R и параметры A, G, X в подпрограммах FA и FI должны иметь тип DOUBLE PRECISION. 
   Вызываемые подпрограммы
       DE04R -
       DE04D   		выполнение одного шага интегрирования линейной системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона.
      UTDE20 -
      UTDE21   		подпрограммы выдачи диагностических сообщений.
  Подпрограммы DE04R, UTDE20 вызываются при работе подпрограммы DE05R, а подпрограммы DE04D, UTDE21 - при работе DE05D. 
   Замечания по использованию
 

Данная подпрограмма предназначена для интегрирования линейных систем, имеющих малую неоднородность  φ (x).

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

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

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

Подпрограммы DE05R и DE05D предназначены также для решения задачи Коши для жестких дифференциальных уравнений (1).

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

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

       0 ≤ x ≤ 6

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

      SUBROUTINE  FA (A, X, M) 
      DIMENSION A(2, 2) 
      A(1, 1) = -(2. + X)/(1. + X) 
      A(1, 2) = 20.*X 
      A(2, 1) = -A(1, 2) 
      A(2, 2) = A(1, 1) 
      RETURN 
      END 

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

      DIMENSION Y(2), YN(2), R(37) 
      EXTERNAL FA, FI 
      M = 2 
      XN = 0. 
      YN(1) = 2. 
      YN(2) = 18. 
      HMIN = 1.E-10 
      EPS = 1.E-10 
      P = 100. 
      XK = 6. 
      H = 0.01 
      CALL  DE05R (FA, FI, M, XN, YN, XK, HMIN, EPS, P, H, Y, R, IERR) 
      END 
 
   Результаты: 

             Y(1)                                Y(2)                             H 
      5.911150077338-03    -2.487346326713-03     1.600000000001-01 

      IERR = 0 
                                    
   2)
        y1'  =  - 20 x y1 + ( (1 + 2x)/(1 + 3x) ) y2  +  ( 1 / 10 ) x2 ,
        y1(0) = 22 
        y2'  =  19 x y1  +  ( (2 + x)/(1 + x) ) y2  -  ( 9 / 10 ) x2 ,
        y2(0) = 18 
              0 ≤ x ≤ 3

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

      SUBROUTINE  FA (A, X, M) 
      DIMENSION A(2, 2) 
      A(1, 1) = -20.*X 
      A(1, 2) = (1. + 2.*X)/(1. + 3.*X) 
      A(2, 1) = 19.*X 
      A(2, 2) = -(2. + X)/(1. + X) 
      RETURN 
      END 

      SUBROUTINE  FI (R1, X, M) 
      DIMENSION R1(2) 
      R1(1) = X*X/10. 
      R1(2) = -9.*R1(1) 
      RETURN 
      END 

      DIMENSION Y(2), YN(2), R(37) 
      EXTERNAL FA, FI 
      M = 2 
      XN = 0. 
      YN(1) = 22. 
      YN(2) = 18. 
      HMIN = 1.E-10 
      EPS = 1.E-10 
      P = 100. 
      XK = 3. 
      H = 0.01 
      CALL  DE05R (FA, FI, M, XN, YN, XK, HMIN, EPS, P, H, Y, R, IERR) 
      END 

   Результаты: 

             Y(1)                                 Y(2)                             H 
      2.134285516536-02     4.227925779755-01     2.500000000001-03 

      IERR = 0