БЧА НИВЦ МГУ. DE94R (DE96R, DE90R, DE92R). Задача Коши для жестких систем и систем с большой константой Липшица

Текст подпрограммы и версий ( Фортран ) de94r.zip, de94d.zip, de96r.zip, de96d.zip, de90r.zip, de90d.zip, de92r.zip, de92d.zip	Тексты тестовых примеров ( Фортран ) tde94r.zip, tde94d.zip, tde96r.zip, tde96d.zip , tde90r.zip, tde90d.zip, tde92r.zip, tde92d.zip
Текст подпрограммы и версий ( Си ) de94r_c.zip, de94d_c.zip, de96r_c.zip, de96d_c.zip, de90r_c.zip, de90d_c.zip, de92r_c.zip, de92d_c.zip	Тексты тестовых примеров ( Си ) tde94r_c.zip, tde94d_c.zip, tde96r_c.zip, tde96d_c.zip , tde90r_c.zip, tde90d_c.zip, tde92r_c.zip, tde92d_c.zip
Текст подпрограммы и версий ( Паскаль ) de94r_p.zip, de94e_p.zip, de96r_p.zip, de96e_p.zip, de90r_p.zip, de90e_p.zip, de92r_p.zip, de92e_p.zip	Тексты тестовых примеров ( Паскаль ) tde94r_p.zip, tde94e_p.zip, tde96r_p.zip, tde96e_p.zip , tde90r_p.zip, tde90e_p.zip, tde92r_p.zip, tde92e_p.zip

Подпрограмма: DE94R (версии: DE96R, DE90R, DE92R)

Назначение

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

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

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

                                Y ' =  F( X, Y ) ,

          Y = ( y₁, ..., y_M ) ,   F = ( f₁ ( X, y₁, ..., y_M ), ..., f_M ( X, y₁, ..., y_M ) )

методом типа Розенброка четвертого порядка.

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

Артемьев С.С., Демидов Г.В. A - устойчивый метод типа Розенброка четвертого порядка точности решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений// Некоторые проблемы вычислительной и прикладной математики. Новосибирск: "Наука", 1975.

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

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

    SUBROUTINE  DE94R (F, FJ, FX, M, JSTART, HMIN, EPS, P, YX, X,
                                            H, BUL, XP, YP, A, A1, A2, RAB,
                                            R1, R2, R3, R4, YD, IERR)

   Параметры

F -	имя подпрограммы вычисления значений правой части дифференциальных уравнений. Первый оператор подпрограммы должен иметь вид: SUBROUTINE F (X, Y, DY, M). Здесь X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY; в случае системы уравнений, т.е. когда M ≠ 1, параметры Y и DY представляют одномерные массивы длиной M (тип параметров X, Y и DY: вещественный);
FJ -	имя подпрограммы вычисления значений элементов матрицы Якоби ∂f / ∂y правой части системы. Первый оператор подпрограммы должен иметь вид: SUBROUTINE FJ (X, Y, Z, M). Здесь X, Y - значения независимой и зависимой переменных, соответственно. В случае системы уравнений, т.е. когда M ≠ 1, параметр Y представляет собой одномерный массив длины M, а параметр Z - двумерный массив размера M * M. Значения элементов матрицы Якоби ∂f / ∂y правой части системы должны быть помещены в массив Z, при этом частная производная от правой части I - го уравнения по J - ой переменной Y ( J ) запоминается в элементе Z ( I, J )(тип параметров X, Y и Z: вещественный);
FX -	имя подпрограммы вычисления значений частных производных по X правой части системы ∂f / ∂x. Первый оператор подпрограммы должен иметь вид: SUBROUTINE FX (X, Y, Z1, M). Здесь: X, Y - значения независимой и зависимой переменных, соответственно. В случае системы уравнений, т.е. когда M ≠ 1, параметры Y и Z1 представляют собой одномерные массивы длиной M. Значения частных производных по X правой части системы ∂f / ∂x должны быть помещены в массив Z1; при этом частная производная от правой части I - го уравнения запоминается в элементе Z1 ( I ) (тип параметров X, Y и Z1: вещественный);
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 - A1, A2	вещественные двумерные рабочие массивы размера M * M;
RAB -	вещественный одномерный рабочий массив длины 5 * M; на выходе из подпрограммы в первых M элементах этого массива содержится погрешность вычисленного приближенного решения, полученная по правилу Рунге;
R1, R2 - R3, R4 YD	вещественные одномерные рабочие массивы длины M;
IERR -	целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров H, HMIN и значением JSTART = - 1 .

   Версии

DE96R -	выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка четвертого порядка. Отличие подпрограммы DE96R от подпрограммы DE94R состоит в том, что при обращении к подпрограмме DE96R не требуется задавать подпрограмму вычисления матрицы Якоби системы и подпрограмму вычисления частных производных правой части системы по x . Все частные производные от правой части вычисляются в подпрограмме DE96R с помощью разностных отношений. Первый оператор подпрограммы DE96R имеет вид: SUBROUTINE DE96R (F, M, JSTART, HMIN, EPS, P, YX, X, H, BUL, XP, YP, A, A1, A2, RAB, R1, R2, R3, R4, YD, IERR) Список формальных параметров подпрограммы DE96R отличается от списка параметров подпрограммы DE94R отсутствием параметров FJ и FX. Параметры подпрограммы DE96R имеют тот же смысл, что и одноименные параметры подпрограммы DE94R.
DE90R -	выполнение одного шага численного интегрирования автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка четвертого порядка. Данная подпрограмма предназначена для интегрирования систем уравнений, в правые части которых не входит независимая переменная x , т.е. систем вида y' = f (y) . Первый оператор подпрограммы DE90R имеет вид: SUBROUTINE DE90R (F, FJ, M, JSTART, HMIN, EPS, P, YX, X, H, BUL, XP, YP, A, A1, A2, RAB, R1, R2, YD, IERR) Список формальных параметров подпрограммы DE90R отличается от списка параметров DE94R отсутствием параметров FX, R3, R4. Параметры подпрограммы DE90R имеют тот же смысл, что и одноименные параметры подпрограммы DE94R, кроме параметра RAB. В подпрограмме DE90R параметр RAB представляет одномерный вещественный рабочий массив длины 4 * M.
DE92R -	выполнение одного шага численного интегрирования жесткой автономной системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка четвертого порядка. Данная подпрограмма предназначена для интегрирования систем уравнений, в правые части которых не входит независимая переменная X, т.е. систем вида y' = f (y) . Отличие подпрограммы DE92R от подпрограммы DE94R состоит в том, что при обращении к подпрограмме DE92R не требуется задавать подпрограмму вычисления матрицы Якоби системы ∂f / ∂y и подпрограмму вычисления частных производных ∂f / ∂x . Первый оператор подпрограммы DE92R имеет вид: SUBROUTINE DE92R (F, M, JSTART, HMIN, EPS, P, YX, X, H, BUL, XP, YP, A, A1, A2, RAB, R1, R2, YD, IERR) Список формальных параметров подпрограммы DE92R отличается от списка параметров подпрограммы DE94R отсутствием параметров FJ, FX, R3, R4. Параметры подпрограммы имеют тот же смысл, что и в подпрограмме DE94R, кроме параметра RAB. В подпрограмме DE92R параметр RAB представляет одномерный вещественный рабочий массив длины 4 * M.
DE94D -	выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка четвертого порядка с повышенной точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE94R; при этом параметры HMIN, EPS, P, YX, X, H, XP, YP, A, A1, A2, RAB, R1, R2, R3, R4, YD и параметры X, Y, DY, Z, Z1 в подпрограммах F, FJ, FX должны иметь тип DOUBLE PRECISION.
DE96D -	выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка четвертого порядка с повышенной точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE96R; при этом параметры HMIN, EPS, P, YX, X, H, XP, YP, A, A1, A2, RAB, R1, R2, R3, R4, YD и параметры X, Y, DY в подпрограмме F должны иметь тип DOUBLE PRECISION.
DE90D -	выполнение одного шага численного интегрирования автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка четвертого порядка с повышенной точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE90R; при этом параметры HMIN, EPS, P, YX, X, H, XP, YP, A, A1, A2, RAB, R1, R2, YD и параметры X, Y, DY, Z в подпрограммах F, FJ должны иметь тип DOUBLE PRECISION.
DE92D -	выполнение одного шага численного интегрирования жесткой автономной системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка четвертого порядка с повышенной точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE92R; при этом параметры HMIN, EPS, P, YX, X, H, XP, YP, A, A1, A2, RAB, R1, R2, YD и параметры X, Y, DY в подпрограмме F должны иметь тип DOUBLE PRECISION.

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

AIG2R -	обращение вещественной матрицы методом Жордана с выбором ведущего элемента по всей матрице. Вызывается при работе подпрограмм DE90R, DE92R, DE94R, DE96R.
AIG2D -	обращение вещественной матрицы, заданной с удвоенной точностью, методом Жордана с выбором ведущего элемента по всей матрице. Вызывается при работе подпрограмм DE90D, DE92D, DE94D, DE96D.
UTDE20 -	подпрограмма выдачи диагностических сообщений при работе подпрограмм DE90R, DE92R, DE94R, DE96R.
UTDE21 -	подпрограмма выдачи диагностических сообщений при работе подпрограмм DE90D, DE92D, DE94D, DE96D.

Кроме того, при работе подпрограммы DE90R вызываются рабочие подпрограммы DE92RS, DE92RK; при работе подпрограммы DE92R вызываются рабочие подпрограммы DE92RS, DE92RK, DE84RJ; при работе подпрограммы DE94R выэываются рабочие подпрограммы DE96RS, DE92RK; при работе подпрограммы DE96R вызываются рабочие подпрограммы DE96RS, DE92RK, DE84RJ, DE96RX; при работе подпрограммы DE90D вызываются рабочие подпрограммы DE92DS, DE92DK; при работе подпрограммы DE92D вызываются подпрограммы DE92DS, DE92DK, DE84DJ; при работе подпрограммы DE94D вызываются подпрограммы DE94DS, DE92DK; при работе подпрограммы DE96D вызываются подпрограммы DE96DS, DE92DK, DE84DJ, DE96DX.

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

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

При работе подпрограммы и ее версий значения параметров M, HMIN, EPS, P сохраняются. При работе подпрограмм F, FJ, FX значения параметров M, X, Y не должны изменяться.

При выполнении из точки x_n одного шага H метода типа Розенброка в подпрограмме DE94R и ее версиях вычисляются четыре значения правой части системы дифференциальных уравнений. Среди этих значений есть одно, которое должно вычисляться при значении независимой переменной x, равном x_n - H. Это значение аргумента x_n находится слева от точки x_n при H > 0 и справа от точки x_n при H < 0. Это следует учитывать при составлении подпрограммы F вычисления правой части системы. Если правая часть системы дифференциальных уравнений не определена для x < x_n при интегрировании слева направо (H > 0) или для x < x_n при интегрировании справа налево (H < 0), то попытка вычислить правую часть для указанного значения аргумента x может привести к авосту.

Если при решении системы дифференциальных уравнений не задаются подпрограммы FJ вычисления матрицы Якоби и FX вычисления частных производных по x правой части системы (т.е. производится обращение к подпрограммам DE92R, DE96R, DE92D, DE96D), то все частные производные по y и по x в этих подпрограммах аппроксимируются центральными разностными отношениями. Замена точных значений производных разностными аппроксимациями может привести к росту погрешности приближенного решения системы дифференциальных уравнений по сравнению с тем, что было бы, если бы все частные производные вычислялись точно с помощью подпрограмм FJ и FX. Даже если будет мала погрешность аппроксимации частных производных (например, если правая часть системы является линейной функцией своих аргументов, то погрешность аппроксимации частных производных равна нулю), может оказаться значительной вычислительная погрешность, особенно, когда вычисления выполняются с одинарной точностью. При этом вычислительная погрешность может даже превосходить верхний предел погрешности приближенного решения, заданный при обращении к подпрограмме параметром EPS. В этом случае целесообразно использовать не подпрограммы DE92R, DE96R, а их версии, выполняющие вычисления с удвоенным числом значащих цифр, т.е. подпрограммы DE92D, DE96D.

Если частные производные по x от правой части системы вычисляются с помощью центральных разностных отношений (т.е. когда производится обращение к подпрограммам DE96R, DE96D), то значения правой части системы, необходимые для аппроксимации этой производной, будут вычисляться с помощью подпрограммы F при значениях аргумента X, находящихся по разные стороны от точки X_n, а именно, при X_n - Δ, X_n + Δ где Δ - значение, задаваемое параметром EPS. Если правая часть системы дифференциальных уравнений не определена для x < x_n при интегрировании слева направо или для x > x_n при интегрировании справа налево, то попытка вычислить правую часть для указанного значения аргумента x может привести к авосту.

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

   Пример 1.
   Решается задача Коши

      y₁' = - 100y₁ 
      y₂' = - 100y₁  - 2y₂  + 20e^{- 100x} + 2e^{- x} cos x 
      y₃' = - 100y₁  + 9998y₂ - 9990y₃ - 10y₄  + 20e^{- 100x} + 2e^{- x} cos x 
      y₄' = - 100y₁  + 9988y₂  + 20y₃ - 10010y₄  + 20e^{- 100x} + 2e^{- x} cos x 

      x  =  0 ,  y₁(0) = 10 ,  y₂(0) = 11 ,  y₃(0) = 111 ,  y₄(0) = 111 . 

   Точное решение задачи имеет вид:

      y₁  = 10e^{- 100x} 
      y₂  = 10e^{- 100x} + e^{- x} cos x + e^{- x} sin x 
      y₃  = y₂  + 100e^{- 10000x} cos 10x 
      y₄  = y₃  + 100e^{- 10000x} sin 10x 

   Матрица Якоби правой части системы имеет вид:

              - 100        0             0              0      
              - 100      - 2             0              0      
              - 100        9998     - 9990     - 10     
              - 100        9988       20         - 10010   

   Частные производные по  x  правой части системы имеют вид:

                        0                 
     - 2000e^{- 100x} - 2e^{- x} (cos x + sin x)   
     - 2000e^{- 100x} - 2e^{- x} (cos x + sin x)   
     - 2000e^{- 100x} - 2e^{- x} (cos x + sin x)

Ниже приводятся фрагмент вызывающей программы для DE94R, подпрограммы F, FJ, FX и результаты счета, полученные после нескольких обращений к подпрограмме DE94R.

       DIMENSION  YX(4), YP(4), A(4, 4), A1(4, 4), A2(4, 4), RAB(20), 
      *                       R1(4), R2(4), R3(4), R4(4), YD(4)
       LOGICAL  BUL
       EXTERNAL  F, FJ, FX 
       M = 4 
       X = 0.
       YX(1) = 10.
       YX(2) = 11.
       YX(3) = 111.
       YX(4) = 111.
       HMIN = 1.E - 10
       EPS = 1.E - 8
       P = 1.E3
       JSTART = 0
       H = 0.01
       IH = 0
  10 IH = IH + 1
       CALL  DE94R (F, FJ, FX, M, JSTART, HMIN, EPS, P, YX, X, H, BUL,
      *                        XP, YP, A, A1, A2, RAB, R1, R2, R3, R4, YD, IERR)
C     BЫЧИCЛEHИE TOЧHOГO PEШEHИЯ CИCTEMЫ: 
       Y1 = 10.*EXP(- 100.*X)
       Y2 = Y1 + EXP(- X)*(COS(X) + SIN(X))
       Y3 = Y2 + 100.*EXP(- 1.E4*X)*COS(10.*X)
       Y4 = Y3 + 100.*EXP(- 1.E4*X)*SIN(10.*X)
       GO TO (11, 12, 13, 14, 15, 16), IH
  11 H = 1.E - 8
       GO TO 10
  12 H = - 1.E - 8
       GO TO 10
  13 JSTART = - 1
       GO TO 10
  14 JSTART = - 1
       H = - 1.E - 8
       GO TO 10
  15 H = 0.01
       GO TO 10
  16 CONTINUE
   
       SUBROUTINE  F (X, Y, DY, M) 
       DIMENSION  Y(4), DY(4) 
       DY(1) = - 100.*Y(1) 
       T1 = 20.*EXP(- 100.*X) + 2.*EXP(- X)*COS(X) 
       DY(2) = DY(1) - 2.*Y(2) + T1 
       DY(3) = DY(1) + 9998.*Y(2) - 9990.*Y(3) - 10.*Y(4) + T1 
       DY(4) = DY(1) + 9988.*Y(2) + 20.*Y(3) - 10010.*Y(4) + T1 
       RETURN 
       END 

       SUBROUTINE  FJ (X, Y, DF, M) 
       DIMENSION  Y(4), DF(4, 4) 
       DF(1, 1) = - 100. 
       DF(1, 2) = 0. 
       DF(1, 3) = 0. 
       DF(1, 4) = 0. 
       DF(2, 1) = - 100. 
       DF(2, 2) = - 2. 
       DF(2, 3) = 0. 
       DF(2, 4) = 0. 
       DF(3, 1) = - 100. 
       DF(3, 2) = 9998. 
       DF(3, 3) = - 9990. 
       DF(3, 4) = - 10. 
       DF(4, 1) = - 100. 
       DF(4, 2) = 9988. 
       DF(4, 3) = 20. 
       DF(4, 4) = - 10010. 
       RETURN 
       END 

       SUBROUTINE  FX (X, Y, DX, M) 
       DIMENSION  Y(4), DX(4) 
       DX(1) = 0. 
       DX(2) = - 2.E3*EXP(- 1.E2*X) - 2.*EXP(- X)*(COS(X) + SIN(X)) 
       DX(3) = DX(2) 
       DX(4) = DX(2) 
       RETURN 
       END 

Результаты:

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

              Y1                                   Y2 
      9.997558891759 + 00     1.099755889174 + 01
              Y3                                   Y4 
      1.085857138761 + 02     1.085880963993 + 02
             YX(1)                               YX(2)
      9.997558891700 + 00     1.099755889169 + 01
             YX(3)                                YX(4)
      1.085857138677 + 02     1.085880963909 + 02
              X                                     H
      2.441406250001 - 06     2.441406250001 - 06

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
    - 9.701276818918 - 13     - 2.910383045670 - 12
            RAB(3)                           RAB(4)
      7.823109626764 - 09       7.784304519497 - 09

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

              Y1                                   Y2 
      9.997548894156 + 00     1.099754889414 + 01
              Y3                                   Y4 
      1.085759455506 + 02     1.085783375935 + 02
             YX(1)                               YX(2)
      9.997548894098 + 00     1.099754889408 + 01
             YX(3)                               YX(4)
      1.085759455420 + 02     1.085783375847 + 02
              X                                     H
      2.451406249999 - 06     2.000000000004 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
      0.000000000000 + 00     0.000000000000 + 00
            RAB(3)                           RAB(4)
    - 7.761021455135 - 12   - 1.552204291027 - 11

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

              Y1                                   Y2 
      9.997558891759 + 00     1.099755889174 + 01
              Y3                                   Y4 
      1.085857138761 + 02     1.085880963993 + 02
             YX(1)                               YX(2)
      9.997558891613 + 00     1.099755889160 + 01
             YX(3)                               YX(4)
      1.085857138667 + 02     1.085880963901 + 02
              X                                     H
      2.441406249998 - 06   - 2.000000000004 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
    - 9.701276818918 - 13   - 9.701276818918 - 13
            RAB(3)                           RAB(4)
    - 2.328306436536 - 11   - 1.552204291027 - 11

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

              Y1                                   Y2 
      9.997568889303 + 00     1.099756888929 + 01
              Y3                                   Y4 
      1.085954831772 + 02     1.085978561790 + 02
             YX(1)                               YX(2)
      9.997568889157 + 00     1.099756888914 + 01
             YX(3)                               YX(4)
      1.085954831680 + 02     1.085978561698 + 02
              X                                     H
      2.431406249996 - 06   - 4.000000000008 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
    - 1.940255363783 - 12   - 1.940255363783 - 12
            RAB(3)                           RAB(4)
    - 1.552204291027 - 11   - 1.552204291027 - 11

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

              Y1                                   Y2 
      9.997558891759 + 00     1.099755889174 + 01
              Y3                                   Y4 
      1.085857138761 + 02     1.085880963993 + 02
             YX(1)                               YX(2)
      9.997558891613 + 00     1.099755889160 + 01
             YX(3)                               YX(4)
      1.085857138667 + 02     1.085880963901 + 02
              X                                     H
      2.441406249998 - 06   - 2.000000000004 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
    - 9.701276818918 - 13   - 9.701276818918 - 13
            RAB(3)                           RAB(4)
    - 2.328306436536 - 11   - 1.552204291027 - 11

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

              Y1                                   Y2 
      9.995118379389 + 00     1.099511837936 + 01
              Y3                                   Y4 
      1.062995982552 + 02     1.062342483763 + 02
             YX(1)                               YX(2)
      9.995118379200 + 00     1.099511837916 + 01
             YX(3)                               YX(4)
      1.062995982379 + 02     1.062342483588 + 02
              X                                     H
      4.882812499996 - 06     2.441406250001 - 06

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
    - 1.940255363783 - 12   - 2.910387045670 - 12
            RAB(3)                           RAB(4)
      7.644606133297 - 09     7.590278983110 - 09

Пример 2.
Решается задача Коши из примера 1. Приводятся фрагмент вызывающей программы для DE96R, подпрограмма F и результаты счета, полученные после нескольких обращений к подпрограмме DE96R.

      DIMENSION  YX(4), YP(4), A(4, 4), A1(4, 4), A2(4, 4), RAB(20), 
     *                       R1(4), R2(4), R3(4), R4(4), YD(4)
      LOGICAL BUL
      EXTERNAL F 
      M = 4 
      X = 0.
      YX(1) = 10.
      YX(2) = 11.
      YX(3) = 111.
      YX(4) = 111.
      HMIN = 1.E - 10
      EPS = 1.E - 8
      P = 1.E3 
      JSTART = 0
      H = 0.01
      IH = 0
 10 IH = IH + 1
      CALL  DE96R (F, M, JSTART, HMIN, EPS, P, YX, X, H, BUL,
     *                        XP, YP, A, A1, A2, RAB, R1, R2, R3, R4, YD, IERR)
C     BЫЧИCЛEHИE TOЧHOГO PEШEHИЯ CИCTEMЫ:
      Y1 = 10.*EXP(- 100.*X)
      Y2 = Y1 + EXP(-X)*(COS(X) + SIN(X))
      Y3 = Y2 + 100.*EXP(- 1.E4*X)*COS(10.*X)
      Y4 = Y3 + 100.*EXP(- 1.E4*X)*SIN(10.*X)
       GO TO (11, 12, 13, 14, 15, 16), IH
 11  H = 1.E - 8
       GO TO 10
 12  H = - 1.E - 8
       GO TO 10
 13  JSTART = - 1
       GO TO 10
 14  JSTART = - 1
      H = - 1.E - 8
       GO TO 10
 15  H = 0.01
       GO TO 10
 16  CONTINUE

      SUBROUTINE  F (X, Y, DY, M)
      DIMENSION  Y(4), DY(4) 
      DY(1) = - 100.*Y(1)
      T1 = 20.*EXP(- 100.*X) + 2.*EXP(- X)*COS(X)
      DY(2) = DY(1) - 2.*Y(2) + T1 
      DY(3) = DY(1) + 9998.*Y(2) - 9990.*Y(3) - 10.*Y(4) + T1
      DY(4) = DY(1) + 9988.*Y(2) + 20.*Y(3) - 10010.*Y(4) + T1
      RETURN 
      END
 
Результаты:

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

              Y1                                   Y2 
      9.999961853100 + 00     1.099996185309 + 01
              Y3                                   Y4 
      1.109618221549 + 02     1.109618602874 + 02
             YX(1)                               YX(2)
      9.999961853042 + 00     1.099996185304 + 01
             YX(3)                               YX(4)
      1.109618221468 + 02     1.109618602339 + 02
              X                                     H
      3.814697265626 - 08     3.814697265626 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
    - 1.940255363783 - 12   - 1.940255363783 - 12
            RAB(3)                           RAB(4)
    - 1.777273913224 - 09     3.523503740628 - 09

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

              Y1                                   Y2
      9.999951853111 + 00     1.099995185310 + 01
              Y3                                   Y4
      1.109518164690 + 02     1.109518645927 + 02
             YX(1)                               YX(2)
      9.999951853039 + 00     1.099995185304 + 01
             YX(3)                               YX(4)
      1.109518164597 + 02     1.109518645362 + 02
              X                                     H
      4.814697265626 - 08     2.000000000004 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
    - 9.701276818918 - 13   - 9.701276818918 - 13
            RAB(3)                           RAB(4)
      1.862645149229 - 10   - 1.552204291027 - 11

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

              Y1                                   Y2
      9.999961853100 + 00     1.099996185309 + 01
              Y3                                   Y4
      1.109618221549 + 02     1.109618602874 + 02
             YX(1)                               YX(2)
      9.999961852940 + 00     1.099996185294 + 01
             YX(3)                               YX(4)
      1.109618221460 + 02     1.109618602293 + 02
              X                                     H
      3.814697265621 - 08   - 2.000000000004 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
    - 1.940255363783 - 12   - 1.940255363783 - 12
            RAB(3)                           RAB(4)
    - 6.208817164108 - 11     6.984919309620 - 11

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

              Y1                                   Y2
      9.999971853074 + 00     1.099997185306 + 01
              Y3                                   Y4
      1.109718288408 + 02     1.109718569798 + 02
             YX(1)                               YX(2)
      9.999971852914 + 00     1.099997185291 + 01
             YX(3)                               YX(4)
      1.109718288386 + 02     1.109718569178 + 02
              X                                     H
      2.814697265622 - 08   - 2.000000000004 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                            RAB(2)
      0.000000000000 + 00      0.000000000000 + 00
            RAB(3)                            RAB(4)
      1.396983861924 - 10      3.182018796603 - 10

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

              Y1                                   Y2
      9.999961853100 + 00      1.099996185309 + 01
              Y3                                   Y4
      1.109618221549 + 02      1.109618602874 + 02
             YX(1)                                YX(2)
      9.999961852940 + 00      1.099996185294 + 01
             YX(3)                                YX(4)
      1.109618221460 + 02      1.109618602293 + 02
              X                                      H
      3.814697265621 - 08    - 2.000000000004 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                            RAB(2)
    - 1.940255363783 - 12    - 1.940255363783 - 12
            RAB(3)                            RAB(4)
    - 6.208817164108 - 11      6.984919309620 - 11

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

              Y1                                   Y2
      9.999923706346 + 00      1.099992370633 + 01
              Y3                                   Y4
      1.109236588570 + 02      1.109237350927 + 02
             YX(1)                                YX(2)
      9.999923706127 + 00      1.099992370613 + 01
             YX(3)                                YX(4)
      1.109236588565 + 02      1.109237350506 + 02
              X                                      H
      7.629394531243 - 08      3.814697265626 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                            RAB(2)
    - 1.940255363783 - 12    - 1.940255363783 - 12
            RAB(3)                            RAB(4)
      4.020209113755 - 09      3.554547826448 - 09

Пример 3.
Решается задача Коши

      y₁' = - 10000y₁  + 100y₂ - 10y₃  + y₄ 
      y₂' = - 1000y₂  + 10y₃ -10y₄ 
      y₃' = - y₃  + 10y₄ 
      y₄' = - 0.1y₄ 

      x  =  0 ,  y₁(0) = 1 ,  y₂(0) = 1 ,  y₃(0) = 1 ,  y₄(0) = 1 

   Точное решение задачи имеет вид:
      y₁  = (1 - d - f - g) t₃  + d t₂  + f t₁  + g y₄ 
      y₂  = (1 - a - b) t₂  + a t₁  + b y₄ 
      y₃  = (- 9.1 / 0.9) t₁ + (10 / 0.9) y₄
      y₄  = e^{- 0.1 x} 
 
  где t = e^{- x}  ,    t = e^{- 1000x}  ,    t = e^{- 10000x}  ,

  a  =  - 91 / (0.9*999) ,    b  =  91 / (0.9*999.9)  ,     d  =  (1 - a - b) / 90  , 
              
  f  =  [100 a + (91/0.9)] / 9999  ,    g  =  [100 b - (100/0.9)  + 1] / 9999.9  .

   Матрица Якоби правой части системы имеет вид:

       - 10000     100    - 10      1     
         0          - 1000     10   - 10    
         0             0        - 1       10    
         0             0          0     - 0.1

Ниже приводятся фрагмент вызывающей программы для DE90R, подпрограммы F и FJ и результаты счета, полученные после нескольких обращений к подпрограмме DE90R.

       DIMENSION  YX(4), YP(4), A(4, 4), A1(4, 4), A2(4, 4), RAB(16), 
       *                    R1(4), R2(4), YD(4)
        LOGICAL  BUL
        EXTERNAL  F, FJ 
        M = 4 
        X = 0.
        YX(1) = 1.
        YX(2) = 1.
        YX(3) = 1.
        YX(4) = 1.
        HMIN = 1.E - 10
        EPS = 1.E - 8
        P = 1.E3
        JSTART = 0
        H = 0.01
        IH = 0
  10  IH = IH + 1
        CALL  DE90R (F, FJ, M, JSTART, HMIN, EPS, P, YX, X, H, BUL, 
       *                        XP, YP, A, A1, A2, RAB, R1, R2, YD, IERR)
C     BЫЧИCЛEHИE TOЧHOГO PEШEHИЯ CИCTEMЫ
        Y4 = EXP(- 0.1*X)
        T1 = EXP(- X)
        Y3 = -(9.1/0.9)*T1 + (10./0.9)*Y4
        T2 = EXP(- 1.E3*X)
        AA = -91./(0.9*999.)
        BB = 91./(0.9*999.9)
        C = 1. - AA - BB
        Y2 = C*T2 + AA*T1 + BB*Y4
        T3 = EXP(- 1.E4*X)
        DD = C/90.
        FF = (100.*AA + 91./0.9)/9999.
        GG = (100.*BB - 100./0.9 + 1.)/9999.9
        QQ = 1. - DD - FF - GG 
        Y1 = QQ*T3 + DD*T2 + FF*T1 + GG*Y4 
        GO TO (11, 12, 13, 14, 15, 16), IH
   11  H = 1.E - 8
        GO TO 10
   12  EPS = 1.E - 2
        GO TO 10
   13  JSTART = - 1
        EPS = 1.E - 8
        GO TO 10
   14  JSTART = - 1
        H = - 1.E - 8
        GO TO 10
   15  H = 0.01
        GO TO 10
   16  CONTINUE

        SUBROUTINE  F (X, Y, DY, M) 
        DIMENSION  Y(4), DY(4) 
        DY(1) = - 1.E4*Y(1) + 1.E2*Y(2) - 10.*Y(3) + Y(4)
        DY(2) = - 1.E3*Y(2) + 10.*Y(3) - 10.*Y(4)
        DY(3) = -Y(3) + 10.*Y(4)
        DY(4) = - 0.1*Y(4)
        RETURN 
        END
        
        SUBROUTINE  FJ (X, Y, DF, M) 
        DIMENSION  Y(4), DF(4, 4) 
        DF(1, 1) = - 1.E4 
        DF(1, 2) = 1.E2 
        DF(1, 3) = - 10.
        DF(1, 4) = 1. 
        DF(2, 1) = 0.
        DF(2, 2) = - 1.E3 
        DF(2, 3) = 10.
        DF(2, 4) = - 10.
        DF(3, 1) = 0.
        DF(3, 2) = 0.
        DF(3, 3) = - 1. 
        DF(3, 4) = 10.
        DF(4, 1) = 0.
        DF(4, 2) = 0.
        DF(4, 3) = 0.
        DF(4, 4) = - 0.1
        RETURN 
        END
    
Результаты:

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

              Y1                               Y2
      9.527772901474 - 01     9.951290901317 - 01
              Y3                               Y4
      1.000043945169 + 00     9.999995117187 - 01
             YX(1)                           YX(2)
      9.527772877191 - 01     9.951290901299 - 01
             YX(3)                           YX(4)
      1.000043945183 + 00     9.999995117150 - 01
              X                                 H
      4.882812500002 - 06     4.882812500002 - 06

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
      2.270887004367 - 09   - 1.212659602364 - 13
            RAB(3)                           RAB(4)
     -3.637978807088 - 13   - 1.818989403544 - 13

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

              Y1                               Y2
      9.526829222541 - 01     9.951191388946 - 01
              Y3                               Y4
      1.000044035131 + 00     9.999995107200 - 01
             YX(1)                           YX(2)
      9.526829198239 - 01     9.951191388900 - 01
             YX(3)                           YX(4)
      1.000044035176 + 00     9.999995107109 - 01
              X                                 H
      4.892812499997 - 06     2.000000000004 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
    - 1.818989403544 - 13   - 1.212659602364 - 13
            RAB(3)                           RAB(4)
    - 2.425319204729 - 13   - 1.212659602364 - 13

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

              Y1                               Y2
      9.524942147509 - 01     9.950992367185 - 01
              Y3                               Y4
      1.000044215194 + 00     9.999995087192 - 01
             YX(1)                           YX(2)
      9.524942123180 - 01     9.950992367112 - 01
             YX(3)                           YX(4)
      1.000044215169 + 00     9.999995087064 - 01
              X                                 H
      4.912812499994 - 06     4.000000000008 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
    - 1.212659602364 - 13   - 1.212659602364 - 13
            RAB(3)                           RAB(4)
    - 1.212659602364 - 13   - 1.818989403544 - 13

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

              Y1                               Y2
      9.523055449436 - 01     9.950793349426 - 01
              Y3                               Y4
      1.000044395201 + 00     9.999995067164 - 01
             YX(1)                           YX(2)
      9.523055425125 - 01     9.950793349353 - 01
             YX(3)                           YX(4)
      1.000044395165 + 00     9.999995067074 - 01
              X                                 H
      4.932812499997 - 06     8.000000000017 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
    - 6.063298011824 - 14   - 6.063298011824 - 14
            RAB(3)                           RAB(4)
    - 3.637978807088 - 13   - 6.063298011824 - 14

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

              Y1                               Y2
      9.527772901474 - 01     9.951290901317 - 01
              Y3                               Y4
      1.000043945169 + 00     9.999995117187 - 01
             YX(1)                           YX(2)
      9.527772877100 - 01     9.951290901245 - 01
             YX(3)                           YX(4)
      1.000043945169 + 00     9.999995117078 - 01
              X                                 H
      4.882812499996 - 06   - 2.000000000004 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
    - 1.818989403544 - 13   - 6.063298011824 - 14
            RAB(3)                           RAB(4)
    - 2.425319204729 - 13   - 1.212659602364 - 13

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

              Y1                               Y2
      9.078026703492 - 01     9.902819081972 - 01
              Y3                               Y4
      1.000087890148 + 00     9.999990234382 - 01
             YX(1)                           YX(2)
      9.078026657144 - 01     9.902819081881 - 01
             YX(3)                           YX(4)
      1.000087890114 + 00     9.999990234237 - 01
              X                                 H
      9.765624999991 - 06     4.882812500002 - 06

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
      2.162717767836 - 09     0.000000000000 + 00
            RAB(3)                           RAB(4)
    - 3.637978807088 - 13   - 6.063298011824 - 14

Пример 4.
Решается задача Коши из примера 3. Приводятся фрагмент вызывающей программы для DE92R, подпрограмма F и результаты счета, полученные после нескольких обращений к подпрограмме DE92R.

       DIMENSION  YX(4), YP(4), A(4, 4), A1(4, 4), A2(4, 4), RAB(16), 
      *                       R1(4), R2(4), YD(4)
       LOGICAL  BUL
       EXTERNAL  F 
       M = 4 
       X = 0.
       YX(1) = 1.
       YX(2) = 1.
       YX(3) = 1.
       YX(4) = 1.
       HMIN = 1.E - 10
       EPS = 1.E - 8 
       P = 1.E3
       JSTART = 0
       H = 0.01
       IH = 0
 10  IH = IH + 1
       CALL  DE92R (F, M, JSTART, HMIN, EPS, P, YX, X, H, BUL,
      *                        XP, YP, A, A1, A2, RAB, R1, R2, YD, IERR)
C     BЫЧИCЛEHИE TOЧHOГO PEШEHИЯ CИCTEMЫ:
       Y4 = EXP(- 0.1*X)
       T1 = EXP(- X)
       Y3 = -(9.1/0.9)*T1 + (10./0.9)*Y4
       T2 = EXP(- 1.E3*X)
       AA = - 91./(0.9*999.9)
       BB = 91./(0.9*999.9)
       C = 1. - AA - BB
       Y2 = C*T2 + AA*T1 + BB*Y4
       T3 = EXP(- 1.E4*X)
       DD = C/90.
       FF = (100.*AA + 91./0.9)/9999.
       GG = (100.*BB - 100./0.9 + 1.)/9999.9
       QQ = 1. - DD - FF - GG 
       Y1 = QQ*T3 + DD*T2 + FF*T1 + GG*Y4 
       GO TO (11, 12, 13, 14, 15, 16), IH
  11  H = 1.E - 8
       GO TO 100
  12  EPS = 1.E - 2
       GO TO 100
  13  JSTART = - 1
       EPS = 1.E - 8
       GO TO 100
  14  JSTART = - 1
       H = - 1.E - 8
       GO TO 100
  15  H = 0.01
       GO TO 100
  16  CONTINUE

       SUBROUTINE  F (X, Y, DY, M) 
       DIMENSION  Y(4), DY(4) 
       DY(1) = - 1.E4*Y(1) + 1.E2*Y(2) - 10.*Y(3) + Y(4)
       DY(2) = - 1.E3*Y(2) + 10.*Y(3) - 10.*Y(4)
       DY(3) = -Y(3) + 10.*Y(4)
       DY(4) = - 0.1*Y(4)
       RETURN 
       END

Результаты:

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

              Y1                               Y2 
      9.527772901474 - 01     9.951298901317 - 01
              Y3                               Y4 
      1.000043945169 + 00     9.999995117187 - 01
             YX(1)                           YX(2)
      9.527772906931 - 01     9.951290903537 - 01
             YX(3)                           YX(4)
      1.000043945183 + 00     9.999995117150 - 01
              X                                 H
      4.882812500002 - 06     4.882812500002 - 06

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
    - 7.571846557159 - 10   - 1.412748436753 - 11
            RAB(3)                           RAB(4)
    - 3.637978807088 - 13   - 1.818989403544 - 13

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

              Y1                               Y2
      9.526829222541 - 01     9.951191388946 - 01
              Y3                               Y4
      1.000044035231 + 00     9.999995107200 - 01
             YX(1)                           YX(2)
      9.526829227989 - 01     9.951191391137 - 01
             YX(3)                           YX(4)
      1.000044035176 + 00     9.999995107109 - 01
              X                                 H
      4.892812499997 - 06     2.000000000004 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
    - 1.212259602364 - 13   - 1.212659602364 - 13
            RAB(3)                           RAB(4)
    - 2.425319204729 - 13   - 1.212659602364 - 13

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

              Y1                               Y2
      9.524942147509 - 01     9.950992367185 - 01
              Y3                               Y4
      1.000044215194 + 00     9.999995087192 - 01
             YX(1)                           YX(2)
      9.524942152921 - 01     9.950992369349 - 01
             YX(3)                           YX(4)
      1.000044215169 + 00     9.999995087064 - 01
              X                                 H
      4.912812499994 - 06     4.000000000008 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
    - 1.212659602364 - 13   - 1.212659602364 - 13
            RAB(3)                           RAB(4)
    - 1.212659602364 - 13   - 1.818989403544 - 13

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

              Y1                               Y2
      9.523055449436 - 01     9.950793349426 - 01
              Y3                               Y4
      1.000044395201 + 01     9.999995067164 - 01
             YX(1)                           YX(2)
      9.523055454865 - 01     9.950793351591 - 01
             YX(3)                           YX(4)
      1.000044395165 + 00     9.999995067074 - 01
              X                                 H
      4.932812499997 - 06     8.000000000017 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
    - 6.063298011814 - 13     0.000000000000 + 00
            RAB(3)                           RAB(4)
    - 3.637978807088 - 13   - 6.063298011824 - 14

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

              Y1                               Y2
      9.527772901474 - 01     9.951290901317 - 01
              Y3                               Y4
      1.000043945169 + 00     9.999995117187 - 01
             YX(1)                           YX(2)
      9.527772906868 - 01     9.951290903482 - 01
             YX(3)                           YX(4)
      1.000043945169 + 00     9.999995117078 - 01
              X                                 H
      4.882812499996 - 06   - 2.000000000004 - 08

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
    - 1.212659602364 - 13   - 6.063298011824 - 14
            RAB(3)                           RAB(4)
    - 2.425319204729 - 13   - 1.212659602364 - 13

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

              Y1                               Y2
      9.078026703492 - 01     9.902819081972 - 01
              Y3                               Y4
      1.000087890148 + 00     9.999990234382 - 01
             YX(1)                           YX(2)
      9.078026512007 - 01     9.902819086356 - 01
             YX(3)                           YX(4)
      1.000087890114 + 00     9.999990234237 - 01
              X                                 H
      9.765624999991 - 06     4.882812500002 - 06

Распечатка контрольного члена погрешности:

            RAB(1)                           RAB(2)
      3.313774262397 - 09   - 1.552204291027 - 11
            RAB(3)                           RAB(4)
    - 3.637978807088 - 13   - 6.063298011824 - 14