Текст подпрограммы и версий ( Фортран )
de23r.zip , de23d.zip , de25r.zip , de25d.zip 		Тексты тестовых примеров ( Фортран )
tde23r.zip , tde23d.zip , tde25r.zip , tde25d.zip
Текст подпрограммы и версий ( Си )
de23r_c.zip , de23d_c.zip , de25r_c.zip , de25d_c.zip 		Тексты тестовых примеров ( Си )
tde23r_c.zip , tde23d_c.zip , tde25r_c.zip , tde25d_c.zip
Текст подпрограммы и версий ( Паскаль )
de23r_p.zip , de23e_p.zip , de25r_p.zip , de25e_p.zip 		Тексты тестовых примеров ( Паскаль )
tde23r_p.zip , tde23e_p.zip , tde25r_p.zip , tde25e_p.zip

Подпрограмма:  DE23R (версия DE25R)

Назначение

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

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

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

          Y' = F (X, Y) ,

          Y = ( y1, ..., yM ) ,
          F = ( f1 (X, y1,..., yM),..., fM (X, y1,..., yM) )
 с начальными условиями, заданными в точке XN:
          Y(XN) = YN ,     YN = ( y10,..., yM0 ) ,

методом Гира. Решение вычисляется в одной точке XK, которая является концом интервала интегрирования.

Метод Гира для нежесткой системы является многошаговым предсказывающе - исправляющим методом Адамса, записанным в форме Нордсика, при этом предсказание и исправление имеют один и тот же порядок.

B случае, когда система уравнений является жесткой, интегрирование осуществляется специальным методом, основанном на методе типа Адамса и использующим якобиан ( ∂F/∂Y ) системы, который вычисляется подпрограммой по формулам численного дифференцирования. При интегрировании данной системы уравнений численное решение проверяется на точность; считается что значение решения в узле  xn вычислено с требуемой точностью ЕРS, если выполняется следующее соотношение: 

                M   
             (  ∑ δI2 )1/2  ≤  EPS  ,
                I=1

где δI - одна из погрешностей следующих типов: абсолютная, относительная или стандартная. При этом под относительной погрешностью приближенного значения I - й компоненты решения в узле  xn подразумевается отношение абсолютной погрешности  eI этого значения в узле  xn к абсолютной величине значения I - й компоненты в предыдущем узле  xn - 1, т.е.  eI / | yIn - 1|, а под стандартной погрешностью - отношение  eI / YPM (I), где 

           YPM(I)  =  max ( | yI0 |, | yI1 |,..., | yIn-1 | )

Tип погрешности специфицируется пользователем при обращении к подпрограмме.

Gear C.W. The automatic integration of ordinary differential equations. Communicatuons of the ACM, 14, 3 (March 1971), 176-179.

Gear C.W., Numerical Initial Value Problems in Ordinary Differential Equations, Prentice - Hall, Englewood Cliffs, N.J., 1971.

Gear C.W., The automatic integration of stiff ordinary differential equations. Information Processing 68, A.J.H.

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

    SUBROUTINE  DE23R (F, M, XN, YN, XK, HMIN, HMAX, EPS, ISTIFJ, 
                                            IORDER, IU, H, Y, YPM, DELTY, RAB,
                                            YP, IERR) 
   Параметры
F - имя подпрограммы вычисления значений правой части дифференциального уравнения. Первый оператоp подпрограммы должен иметь вид:
SUBROUTINE  F (X, Y , DY, M).
Здесь: X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY. B случае системы уравнений, т.е. когда M ≠ 1, параметры Y и DY представляют массивы длины M (тип параметров X, Y и DY: вещественный);
M - количество уравнений в системе (тип: целый);
XN, YN - начальные значения аргумента и решения; в случае системы уравнений (т.е. M ≠ 1) YN представляет массив длины M (тип: вещественный);
XK - значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования); XK может быть больше, меньше, или pавно XN (тип: вещественный);
HMIN - минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений; это значение должно быть много меньше среднего ожидаемого шага интегрирования, задаваемого параметром H (тип: вещественный);
HMAX - максимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
EPS - допустимая погрешность, с которой требуется вычислить все компоненты решения; тип погрешности специфицируется с помощью параметpа IU (тип: вещественный);
ISTIFJ - целый указатель метода численного интегрирования:
ISTIFJ=0 - интегрирование системы ведется методом Адамса;
ISTIFJ=1 - интегрирование ведется специальным методом, предназначенным для жестких систем;
IORDER - целая переменная, указывающая максимальный допустимый порядок метода; IORDER должен быть не больше 7 для метода Адамса и не больше 6 для метода интегрирования жестких систем;
IU - целый указатель типа погрешности численного решения:
IU = 1 - для стандартной погрешности;
IU = 2 - для относительной погрешности;
IU = 3 - для абсолютной погрешности;
H - вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления итегрирования, т.е. положительным, если XK > XN, отрицательным, если XK < XN, или без такого учета в виде абсолютной величины;
Y - искомое решение задачи Коши, вычисленное подпрограммой для значения аргумента XK; для системы уравнений (когда M ≠ 1) задается массивом длиной M; в случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: вещественный);
         YPM -
       DELTY   		одномерные вещественные рабочие массивы длиной M;
RAB - одномерный вещественный рабочий массив; при интегрировании нежесткой системы уравнений RAB имеет размер 17*M, при интегрировании жесткой системы - M*(M + 17);
YP - двумерный вещественный рабочий массив размеpа 8*M;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом:
IERR= 1 - когда неправильно задан параметр IORDER, а именно, когда IORDER превосходит максимальный допустимый порядок метода; в этом случае интегрирование системы ведется методом Гира порядка не выше 7 для нежесткой системы, и не выше 6 для жесткой;
IERR=65 - когда решение системы не может быть вычислено с требуемой точностью EPS при заданных начальном шаге H, его минимальном значении HMIN и порядке метода IORDER;
IERR=66 - когда приближенное значение решения не может быть вычислено, т.к. итерационный процесс его определения не сходится для шагов интегрирования H, больших заданного минимального значения HMIN;
IERR=67 - когда требуемая точность EPS вычисления приближенного решения меньше той, которая может быть достигнута для данной задачи при тех размерах шага интегрирования, начальное значение которого задано параметром H;
  при IERR = 65, 66, 67 интегрирование системы прекращается; при желании интегрирование можно повторить обращением к подпрограмме с новыми значениями параметров H, HMIN и IORDER;
IERR=68 - когда приближенное значение решения для жесткой системы не может быть вычислено с заданной точностью; для достижения тебуемой точности следует воспользоваться версиями подпрограммы DE23D, DE25R или DE25D. 
   Версии
DE23D - вычисление решения задачи Коши для нежесткой и жесткой систем обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Гира с повышенной точностью. При этом параметры XN, YN, XK, HMIN, HMAX, EPS, H, Y, YPM, DELTY, RAB, YP и параметры X, Y и DY в подпрограмме F должны иметь тип DOUBLE PRECISION.
DE25R - вычисление решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Гира с автоматическим выбором шага. Первый оператор подпрограммы имеет вид:
SUBROUTINE  DE25R (F, FJ, M, XN, YN, XK, HMIN, HMAX, EPS, IORDER, IU, H, Y, YPM, DELTY, RAB, YP, IERR)
Здесь: FJ - имя подпрограммы вычисления якобиана правой части системы; первый оператор этой подпрограммы имеет вид:
SUBROUTINE  FJ (X, Y, Z, M)
Здесь: X, Y - значения независимой и зависимой переменных, соответственно, причем Y представляет одномерный массив длины M; вычисленное значение якобиана должно быть помещено в двумерный массив Z размера M*M, при этом частная производная от правой части I - ого уравнения по J - ой переменной Y (J) запоминается в элементе Z (I, J) (тип параметров X, Y и Z: вещественный). Остальные параметры подпрограммы DE25R имеют тот же смысл, что и одноименные параметры подпрограммы DE23R.
DE25D - вычисление решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Гира с повышенной точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE25R; при этом параметры XN, YN, XK, HMIN, HMAX, EPS, H, Y, YPM, DELTY, RAB, YP и параметры X, Y, DY и Z в подпрограммах F и FJ должны иметь тип DOUBLE PRECISION. 
   Вызываемые подпрограммы
DE21R - выполнение одного шага численного интегрирования нежесткой и жесткой систем обыкновенных дифференциальных уравнений первого порядка методом Гира с контролем точности; вызывается при работе подпрограммы DE23R.
DE21D - выполнение одного шага численного интегрирования нежесткой и жесткой систем обыкновенных дифференциальных уравнений первого порядка методом Гира с повышенной точностью; вызывается при работе подпрограммы DE23D.
DE24R - выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом Гира с контролем точности; вызывается при работе подпрограммы DE25R.
DE24D - выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом Гира с повышенной точностью; вызывается при работе подпрограммы DE25D.
UTDE12 - подпрограмма выдачи диагностических сообщений при работе подпрограмм DE21R, DE23R, DE24R, DE25R.
UTDE13 - подпрограмма выдачи диагностических сообщений при работе подпрограмм DE21D, DE23D, DE24D, DE25D. 
   Замечания по использованию
 

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

При работе подпрограммы и ее версий значения параметров M, XN, YN, XK, HMIN, HMAX, EPS, ISTIFJ, IORDER, IU сохраняются. Если после работы подпрограммы нет необходимости иметь начальное значение решения YN, то параметры YN и Y при обращении к ней можно совместить.

Значение HMIN должно быть много меньше среднего ожидаемого шага интегрирования, задаваемого параметром H при обращении к подпрограмме и ее версиям, т.к. интегрирование системы начинается с метода первого порядка.

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

1)        y1'  =   y3
           y2'  =   y4
           y3'  =  - ( 1 - 0.002 x ) y1 / ( y12  +  y22 )3/2
           y4'  =  - ( 1 - 0.002 x ) y2 / ( y12  +  y22 )3/2  ,  0 ≤ x ≤ 6π ,

          y1(0)  =  0.09411764706 ,     y2(0)  =  0 ,     y3(0)  =  0 ,     y4(0)  =  4.5

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

          SUBROUTINE  F (X, Y, DY, M)
          DIMENSION  Y(4), DY(4)
          R2 = Y(1)*Y(1) + Y(2)*Y(2)
          R2 = R2*SQRT(R2)
          R = - (1. - 0.002*X)/R2
          DY(1) = Y(3)
          DY(2) = Y(4)
          DY(3) = R*Y(1)
          DY(4) = R*Y(2)
          RETURN
          END

          DIMENSION  RAB(68), Y(4), YN(4), DELTY(4), YPM(4), YP(32)
          EXTERNAL  F
          M = 4
          XN = 0.
          YN(1) = 0.09411764706
          YN(2) = 0.
          YN(3) = 0.
          YN(4) = 4.5
          XK = 6.*3.14159265359
          HMIN = 1.E-18
          HMAX = 0.1
          EPS = 1.E-9
          ISTIFJ = 0
          IORDER = 7
          IU = 3
          H = 0.01
          CALL  DE23R (F, M, XN, YN, XK, HMIN, HMAX, EPS, ISTIFJ,
         *                        IORDER, IU, H, Y, YPM, DELTY, RAB, YP, IERR)

Результаты:

          IERR  =  0

          Y(1)  =  -9.776495884-01
          Y(2)  =  -4.398340941-01
          Y(3)  =   9.302044292-01
          Y(4)  =  -1.472281399-02

2)         y1'  =  -20 y1  +  y2
            y2'  =   19 y1 - 2 y2 ,     0 ≤ x ≤ 1

           y1(0)  =  2 ,      y2(0)  =  18
 
 Собственные значения якобиевой матрицы системы
   λ1  =  -21 ,    λ2  =  -1 .
 Поэтому эту систему можно рассматривать как жесткую.
    
          SUBROUTINE  F (X, Y, DY, M)
          DIMENSION  Y(2), DY(2)
          DY(1) = -20.*Y(1) + Y(2)
          DY(2) = 19.*Y(1)-2.*Y(2)
          RETURN
          END
 
          DIMENSION  RAB(38), Y(2), YN(2), DELTY(2), YPM(2), YP(16)
          EXTERNAL  F
          M = 2
          XN = 0.
          YN(1) = 2.
          YN(2) = 18.
          XK = 1.
          HMIN = 1.E-10
          HMAX = 0.1
          EPS = 0.00001
          ISTIFJ = 1
          IORDER = 6
          IU = 3
          H = 0.01
          CALL  DE23R (F, M, XN, YN, XK, HMIN, HMAX, EPS, ISTIFJ,
         *                        IORDER, IU, H, Y, YPM, DELTY, RAB, YP, IERR)

Результаты:

          IERR  =  0

          Y(1)  =  3.678796482-01
          Y(2)  =  6.989709436

3) Pешается пример № 2 с помощью подпрограммы DE25R. Приводится только подпрограмма вычисления якобиана и фрагмент вызывающей программы. 

          SUBROUTINE  FJ (X, Y, PW, M)
          DIMENSION  Y(2), PW(2, 2)
          PW(1, 1) = -20.
          PW(1, 2) = 1.
          PW(2, 1) = 19.
          PW(2, 2) = -2.
          RETURN
          END
    
          DIMENSION  RAB(38), Y(2), YN(2), DELTY(2), YPM(2), YP(16)
          EXTERNAL  F, FJ
          M = 2
          XN = 0.
          YN(1) = 2.
          YN(2) = 18.
          XK = 1.
          HMIN = 1.E-10
          HMAX = 0.1
          EPS = 0.00001
          IORDER = 6
          IU = 3
          H = 0.01
          CALL  DE25R (F, FJ, M, XN, YN, XK, HMIN, HMAX, EPS, IORDER,
         *                        IU, H, Y, YPM, DELTY, RAB, YP, IERR)

Результаты:

          IERR  =  0

          Y(1)  =  3.678796482-01
          Y(2)  =  6.989709436