Текст подпрограммы и версий ( Фортран )
de04r.zip , de04d.zip
Тексты тестовых примеров ( Фортран )
tde04r.zip , tde04d.zip
Текст подпрограммы и версий ( Си )
de04r_c.zip , de04d_c.zip
Тексты тестовых примеров ( Си )
tde04r_c.zip , tde04d_c.zip
Текст подпрограммы и версий ( Паскаль )
de04r_p.zip , de04e_p.zip
Тексты тестовых примеров ( Паскаль )
tde04r_p.zip , tde04e_p.zip

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

Назначение

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

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

Выполняется один шаг численного интегрирования линейной устойчивой системы 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) ) 

Предполагается, что среди характеристических корней матрицы A(X) имеются большие по модулю корни, а функция  φ(x) является достаточно малой. По заданному значению решения YX в узле  xn  вычисляется значение этого решения в узле  xn + H  . Вычисление производится по методу Лоусона.

Метод Лоусона заключается в следующем. Исходная система уравнений с помощью замены искомой функции Y (x) на  [ xn, 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) 

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

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

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

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

    SUBROUTINE  DE04R (FA, FI, M, JSTART, HMIN, EPS, P, YX, X, H,
                                            BUL, XP, YP, A, E1, E2, E3, W1,
                                            R1, R2, R3, R4, R5, R6, YD, 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 - количество уравнений в системе (тип: целый);
JSTART - целый указатель режима использования подпрограммы, имеющий следующие значения:
0,+1 - выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных и шага интегрирования, заданных параметрами X, YX, и H, соответственно;
-1 - повторить последний шаг интегрирования с новыми значениями параметров H и/или HMIN; на выходе из подпрограммы JSTART полагается равным  + 1;
HMIN - минимальное значение абсолютной величины шага, которое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
EPS - допустимая мера погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный);
P - граница перехода, используемая при оценке меры погрешности решения (тип: вещественный);
YX, X - заданные вещественные значение решения и соответствующее ему значение аргумента; в результате работы подпрограммы в X получается новое значение аргумента, в YX - соответствующее значение решения; в случае системы уравнений, т.е. когда M ≠ 1, YX задается одномерным массивом длины M;
H - вещественная переменная, содержащая значение шага интегрирования; если для этого значения шага точность приближенного решения достигается, то именно он и реализуется подпрограммой, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность EPS; на выходе из подпрограммы H содержит рекомендуемое подпрограммой значение следующего шага интегрирования, определяемое ею с целью достижения более экономного способа интегрирования;
BUL - логическая переменная, значение которой при обращении к подпрограмме полагается равным .TRUE., если заданный в H шаг выводит в конец интервала интегрирования, и .FALSE. в противном случае; в результате работы подпрограммы BUL равно .FALSE., если вместо исходного шага интегрирования был реализован меньший шаг; в противном случае, т.е. когда был выполнен именно заданный при обращении в H шаг, значение параметры BUL не меняется;
XP, YP - вещественная рабочая переменная и одномерный рабочий массив длины M, соответственно; значения параметров XP, YP на выходе из подпрограммы равны тем значениям, которые имели параметры X, YX при входе в нее (т.е. предыдущий узел и решение в нем);
       A, E1 -
         E2, E3, 
         W1  
вещественные двумерные рабочие массивы размера M*M;
       R1, R2 -
       R3, R4, 
         R5, R6, 
         YD  
вещественные одномерные рабочие массивы длины M;
IERR - целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров H, HMIN и значением JSTART = - 1 .

Версии

DE04D - выполнение одного шага численного интегрирования линейной устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона с удвоенным числом значащих цифр. При этом параметры HMIN, EPS, P, YX, X, H, XP, YP, A, E1, E2, E3, W1, R1, R2, R3, R4, R5, R6, YD и параметры A, G, X в подпрограммах FA и FI должны иметь тип DOUBLE PRECISION.

Вызываемые подпрограммы

AME2R - подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы DE04R.
AME2D - подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы DE04D.
UTDE20 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE04R.
UTDE21 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE04D.
 

Кроме того при работе подпрограмм DE04R и DE04D вызываются подпрограммы DE04RD, DE04RM, DE04RP, DE04RQ и DE04DD, DE04DM, DE04DP, DE04DQ соответственно.

Замечания по использованию

 

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

В общем случае заданная точность не гарантируется. При работе подпрограммы значения параметров M, HMIN, EPS, P сохраняются.

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

При обращении к подпрограмме со значением JSTART = - 1 в качестве исходных значений аргумента и решения принимаются значения параметров XP, YP, соответственно, т.е. те значения, которые эти параметры получили после самого последнего обращения к подпрограмме с неотрицательным значением JSTART.

После работы подпрограммы в массиве R1 содержится значение оценки погрешности на шаге, вычисленной по правилу Рунге.

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

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

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

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

      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 A(2, 2), E1(2, 2), E2(2, 2), E3(2, 2), W1(2, 2), R1(2),  
     *                      R2(2), R3(2), R4(2), R5(2), R6(2), YD(2), YP(2), YX(2) 
      LOGICAL BUL 
      EXTERNAL FA, FI 
      M = 2 
      X = 0. 
      YX(1) = 0. 
      YX(2) = 1. 
      HMIN = 1.E-10 
      EPS = 1.E-8 
      P = 100. 
      JSTART = 1 
      H = 0.01 
      IH = 0 
 10 IH = IH + 1 
      CALL  DE04R (FA, FI, M, JSTART, HMIN, EPS, P, YX, X, H, BUL,
     *                         XP, YP, A, E1, E2, E3, W1, R1, R2, R3,
     *                         R4, R5, R6, YD, IERR) 
      GO TO (11, 12, 13, 14, 15, 16, 17), IH 
 11 EPS = 1.E-10 
      GO TO 10 
 12 EPS = 1.0 
      GO TO 10 
 13 JSTART = -1 
      H = 0.01 
      EPS = 1.E-5 
      GO TO 10 
 14 JSTART = -1 
      H = -0.0001 
      EPS = 1.E-2 
      GO TO 10 
 15 H = 1. 
      EPS = 1.E-10 
      GO TO 10 
 16 EPS = 1.E-4 
      GO TO 10 
 17 CONTINUE 

  Результаты: 

   После первого обращения к подпрограмме
                 X                                   YX(1)                           YX(2) 
         1.000000000001-02     9.802471967628-04     9.802468700200-01 
                 H                                   R1(1)                            R1(2) 
         2.000000000001-02     0.000000000000+00     6.063298011814-14 

   после второго обращения к подпрограмме
                 X                                   YX(1)                           YX(2) 
         3.000000000000-02     8.479506692396-03     9.421419716109-01 
                 H                                   R1(1)                            R1(2) 
         4.000000000002-02     1.231607408649-14     1.333925562599-12 

   после третьего обращения к подпрограмме
                 X                                   YX(1)                           YX(2) 
         6.999999999994-02     4.268132414620-02     8.703501915807-01 
                 H                                   R1(1)                            R1(2) 
         8.000000000004-02     1.970571853839-12     3.565219230945-11 

   после четвертого обращения к подпрограмме
                 X                                   YX(1)                           YX(2) 
         3.999999999996-02     1.478074532271-02     9.237177506848-01 
                 H                                   R1(1)                            R1(2) 
         2.000000000001-02     0.000000000000+00     6.063298011814-14 

   после пятого обращения к подпрограмме
                 X                                   YX(1)                           YX(2) 
         2.990000000000-02     8.424732657545-03     9.423281849631-01 
                 H                                   R1(1)                            R1(2) 
        -2.000000000004-04     0.000000000000+00     6.063298011814-14 

   после шестого обращения к подпрограмме
                 X                                   YX(1)                           YX(2) 
         6.115000000000-02     3.314040588361-02     8.858545422663-01 
                 H                                   R1(1)                            R1(2) 
         3.125000000000-02     4.395891058563-13     1.097456940137-11 

   после седьмого обращения к подпрограмме
                 X                                   YX(1)                           YX(2) 
         9.240000000000-02     7.117143016819-02     8.315812963056-01 
                 H                                   R1(1)                            R1(2) 
         6.250000000000-02     7.80649190203-13      8.852415097246-12