Текст подпрограммы и версий ( Фортран )
de21r.zip , de21d.zip , de24r.zip , de24d.zip 		Тексты тестовых примеров ( Фортран )
tde21r.zip , tde21d.zip , tde24r.zip , tde24d.zip
Текст подпрограммы и версий ( Си )
de21r_c.zip , de21d_c.zip , de24r_c.zip , de24d_c.zip 		Тексты тестовых примеров ( Си )
tde21r_c.zip , tde21d_c.zip , tde24r_c.zip , tde24d_c.zip
Текст подпрограммы и версий ( Паскаль )
de21r_p.zip , de21e_p.zip , de24r_p.zip , de24e_p.zip 		Тексты тестовых примеров ( Паскаль )
tde21r_p.zip , tde21e_p.zip , tde24r_p.zip , tde24e_p.zip

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

Назначение

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

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

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

          Y' = F (X, Y) ,

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

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

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

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

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

где eI - вычисляемая подпрограммой оценка абсолютной погрешности приближенного значения I - й компоненты решения в узле Xn + H, YPM (I) - задается при обращении к подпрограмме и, по существу, позволяет использовать разный тип погрешности для каждой компоненты. Если, например, YPM(I) = 1, то для I - ой компоненты YX(I) решения будет использоваться абсолютная погрешность; если YPM (I) = | YX (I) | ≠ 0 , где YX(I) - значение I - й компоненты YX(I) в узле Xn, то для этой компоненты берется отношение абсолютной погрешности приближенного ее значения в узле Xn + 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) 
   Параметры
F - имя подпрограммы вычисления значений правой части дифференциального уравнения. Первый оператоp подпрограммы должен иметь вид:
SUBROUTINE  F (X, Y , DY, M).
Здесь: X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY. B случае системы уравнений, т.е. когда M ≠ 1, параметры Y и DY представляют массивы длины M (тип параметров X, Y и DY: вещественный);
M - количество уравнений в системе (тип: целый);
ISTIFJ - целый указатель метода численного интегрирования:
ISTIFJ=0 - интегрирование системы ведется методом Адамса;
ISTIFJ=1 - интегрирование ведется специальным методом, предназначенным для жестких систем;
IORDER - целая переменная, указывающая максимальный допустимый порядок метода; IORDER должен быть не больше 7 для метода Адамса и не больше 6 для метода интегрирования жестких систем;
JSTART - целый указатель режима использования подпрограммы, имеющий следующие значения:
0 - первое обращение к подпрограмме должно быть выполнено с нулевым значением JSTART;
+1 - выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных и шага интегрирования, заданных параметрами X, YX и H ,соответственно;
-1 - повторить последний шаг интегрирования с новыми значениями параметров H и/или HMIN;
  на выходе из подпрограммы JSTART pавен текущему порядку метода;
HMIN - минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений; это значение должно быть много меньше среднего ожидаемого шага интегрирования при первом обращении к подпрограмме с JSTART = 0 (тип: вещественный);
HMAX - максимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
EPS - требуемая точность вычисления решения; размер шага интегрирования и порядок метода выбираются в подпрограмме автоматически таким образом, чтобы вычисляемые ею оценки абсолютных погрешностей всех компонент решения, деленные на YPM (I), были не больше EPS в евклидовой ноpме (тип: вещественный);
X, YX - начальные вещественные значения аргумента и решения уравнений в нем; в результате работы подпрограммы в X получается новое значение аpгумента, а в YX - соответствующее значение pешения; в случае системы уравнений, т.е. когда M ≠ 1 , YX задается одномерным массивом длины M;
H - вещественная переменная, содержащая значение шага интегрирования; если для этого значения шага точность приближенного решения достигается, то именно он и реализуется подпрограммой, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность; на выходе из подпрограммы H содержит рекомендуемое подпрограммой значение следующего шага интегрирования, определяемое ею с целью достижения наиболее экономного способа интегрирования;
BUL - логическая переменная, значение которой при обращении к подпрограмме полагается равным .TRUE., если заданный в H шаг выводит в конец интервала интегрирования, и .FALSE. в противном случае; в результате работы подпрограммы BUL pавно .FALSE., если вместо исходного шага интегрирования был реализован меньший шаг; в противном случае, т.е. когда был выполнен именно заданный при обращении в H шаг, значение параметра BUL не меняется;
YPM - одномерный вещественный массив длины M; значение YPM, котоpое он имеет на входе в подпрограмму, используется при вычислении погрешности приближенного решения; считается, что приближенное решение достигает требуемой точности, если евклидова ноpма вектоpа, составленного из абсолютных погрешностей вычисленных значений всех компонент решения, деленных на соответствующие элементы массива YPM (т.е. абсолютная погрешность I - й компоненты делится на I - й элемент YPM(I)), не превосходит EPS; в pезультате работы подпрограммы значение I - го элемента YPM (I) сохраняется, если вновь вычисленное значение I - й компоненты YX (I) решения не превосходит по абсолютной величине исходного значения YPM (I), и заменяется на абсолютную величину нового значения, если она больше первоначального значения YPM (I); если на входе в подпрограмму YPM (I) = 1, то для I - ой компоненты решения YX (I) будет использоваться абсолютная погрешность; если на входе YPM (I) = | YX (I) | ≠ 0, то для I - й компоненты берется отношение абсолютной погрешности приближенного ее значения в узле X + H к абсолютной величине значения этой компоненты в предыдущем узле X;
DELTY - одномерный вещественный рабочий массив длины M;
RAB - одномерный вещественный рабочий массив; при интегрировании нежесткой системы уравнений RAB имеет размер 17*M, при интегрировании жесткой системы - M*(M + 17);
YP - двумерный вещественный рабочий массив размера 8*M;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом:
IERR= 1 - когда неправильно задан параметр IORDER, а именно, когда IORDER превосходит максимальный допустимый порядок метода; в этом случае интегрирование ведется методом Гира порядка не выше 7 для нежесткой системы и не выше 6 для жесткой;
IERR=65 - когда интегрирование системы выполнено с заданным минимальным шагом, но требуемая точность полученного при этом значения YX решения не достигнута; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров H и HMIN и со значением JSTART = - 1;
IERR=66 - когда приближенное значение решения не может быть вычислено, т.к. итерационный процесс его определения не сходится для шагов интегрирования H, больших заданного минимального значения HMIN;
IERR=67 - когда требуемая точность EPS вычисления приближенного решения меньше той, которая может быть достигнута для данной задачи при тех размерах шага интегрирования, начальное значение которого задано параметром H;
  при IERR = 66 и 67 значения параметров X и YX сохраняются, а в H засылается то значение шага интегрирования, котоpое было на выходе в предыдущем обращении к подпрограмме; в этом случае интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров H и HMIN, если IERR = 66, и параметра H, если IERR = 67;
IERR=68 - когда приближенное значение решения для жесткой системы не может быть вычислено с заданной точностью; в этом случае значения X и YX сохраняются, а в H засылается то значение шага интегрирования, котоpое было на выходе в предыдущем обращении к подпрограмме; для достижения требуемой точности следует воспользоваться версиями подпрограммы DE21D, DE24R, или DE24D. 
   Версии
DE21D - выполнение одного шага численного интегрирования нежесткой и жесткой систем обыкновенных дифференциальных уравнений первого порядка методом Гира с повышенной точностью. При этом параметры HMIN, HMAX, EPS, YX, X, H, YPM, DELTY, RAB, YP и параметры X, Y и DY в подпрограмме F должны иметь тип DOUBLE PRECISION.
DE24R - выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом Гира с контролем точности. Первый оператор подпрограммы имеет вид:
SUBROUTINE  DE24R ( F, FJ, M, IORDER, JSTART, HMIN, HMAX, EPS, YX, X, H, BUL, 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: вещественный). Остальные параметры подпрограммы DE24R имеют тот же смысл, что и одноименные параметры подпрограммы DE21R.
DE24D - выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом Гира с повышенной точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE24R; при этом паметры HMIN, HMAX, EPS, YX, X, H, YPM, DELTY, RAB, YP и параметры X, Y, DY и Z в подпрограммах F и FJ должны иметь тип DOUBLE PRECISION. 
   Вызываемые подпрограммы
UTDE12 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE21R и DE24R.
UTDE13 - подпрограмма выдачи диагностических сообщений при работе подпрограмм DE21D и DE24D.
  Kpоме того, DE21R, DE24R и DE21D, DE24D используют pабочие подпрограммы DE21RU, DE21RP и DE21DU, DE21DP, соответственно. 
   Замечания по использованию
 

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

При работе подпрограммы и ее версий значения параметров M, ISTIFJ, IORDER, HMIN, HMAX, EPS сохраняются.

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

При многократном использовании подпрограммы или ее веpсий для вычисления решения системы уравнений на отрезке значения параметров M, ISTIFJ, EPS, YX, X и значения рабочих массивов, задаваемых параметрами YPM, DELTY, RAB, YP, не должны изменяться в вызывающей программе между последовательными обращениями к подпрограмме DE21R или ее версиям.

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

         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