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

Назначение

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

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

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

          Y' = F (X, Y) ,

          Y = ( y1, ..., yM ) ,
          F = ( f1 (X, y1,..., yM),..., fM (X, y1,..., yM) )     методом Гира.

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

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

По заданному значению YX решения в узле X_n вычисляется значение этого решения в узле X_n + H. Считается, что решение в узле X_n + H вычислено с требуемой точностью EPS, если выполняется следующее соотношение:

                M   
             (  ∑ ( eI / YPM(I) )2 )1/2  ≤  EPS  ,
                I=1 

где e_I - вычисляемая подпрограммой оценка абсолютной погрешности приближенного значения I - й компоненты решения в узле X_n + H, YPM (I) - задается при обращении к подпрограмме и, по существу, позволяет использовать разный тип погрешности для каждой компоненты. Если, например, YPM(I) = 1, то для I - ой компоненты YX(I) решения будет использоваться абсолютная погрешность; если YPM (I) = | YX (I) | ≠ 0 , где YX(I) - значение I - й компоненты YX(I) в узле X_n, то для этой компоненты берется отношение абсолютной погрешности приближенного ее значения в узле X_n + H к абсолютной величине приближенного значения этой компоненты в предыдущем узле Х_n. Значение H может быть меньше или pавно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.

Gear C.W., The automatic integration of ordinary differential equations. Communications 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  DE21R (F, M, ISTIFJ, IORDER, JSTART, HMIN, HMAX,
                                           EPS, YX, X, H, BUL, YPM, DELTY,
                                           RAB, YP, IERR) 

   Параметры 

   Версии 

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

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

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

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

        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), DELTY(4), YPM(4), YP(32)
          LOGICAL  BUL
          EXTERNAL  F
          M = 4
          ISTIFJ = 0
          IORDER = 4
          JSTART = 0
          HMIN = 1.E-18
          HMAX = 0.1
          EPS = 0.00001
          Y(1) = 0.09411764706
          Y(2) = 0.
          Y(3) = 0.
          Y(4) = 4.5
          X = 0.
          H = 0.01
          BUL = .FALSE.
          DO 10  I = 1, M
    10  YPM(I) = 1.
          CALL  DE21R (F, M, ISTIFJ, IORDER, JSTART, HMIN, HMAX,
         *                        EPS, Y, X, H, BUL, YPM, DELTY, RAB, YP, IERR)

Результаты:

          X  =   4.057708603-05
          H  =   4.057708603-05

          Y(1)  =   9.411746119-02
          Y(2)  =   1.825965265-04
          Y(3)  =  -4.580764462-03
          Y(4)  =   4.499991113

          IERR  =   0