Текст подпрограммы и версий ( Фортран )
de34r.zip , de34d.zip , de18r.zip , de18d.zip 		Тексты тестовых примеров ( Фортран )
tde34r.zip , tde34d.zip , tde18r.zip , tde18d.zip
Текст подпрограммы и версий ( Си )
de34r_c.zip , de34d_c.zip , de18r_c.zip , de18d_c.zip 		Тексты тестовых примеров ( Си )
tde34r_c.zip , tde34d_c.zip , tde18r_c.zip , tde18d_c.zip
Текст подпрограммы и версий ( Паскаль )
de34r_p.zip , de34e_p.zip , de18r_p.zip , de18e_p.zip 		Тексты тестовых примеров ( Паскаль )
tde34r_p.zip , tde34e_p.zip , tde18r_p.zip , tde18e_p.zip

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

Назначение

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

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

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

    Y '' = F (X, Y)   ,
    Y = ( y1, ..., yM )   ,
    F = ( f1 (X, y1, ..., yM), ..., fM (X, y1, ..., yM) )   ,

с начальными условиями, заданными в точке  XN: 

          Y (XN) = YN,      YN = ( y10, ..., yM0 )   ,
          Y '(XN) = DYN,  DYN = ( y10', ..., yM0' )   ,

отыскиваются необходимые при численном интегрировании этой системы многошаговым методом Штермера значение решения  Y в первом (после начальной точки  XN) узле интегрирования  ХN + Н, его конечная разность первого порядка и конечные разности правой части системы до четвертого порядка включительно в том же узле  XN + H.

Четвертый порядок разностей соответствует шестому порядку точности используемого метода Штермера.

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

Используемый в подпрограмме метод разгона является итерационным способом, опирающимся на пару формул, одна из которых представляет значение решения в узле  XN + H через начальные условия (т.е. через его значение и значение его производной в начальной точке  XN), а также через конечные разности правой части системы в точке  XN, а другая является экстраполяционной формулой, применяемой для предсказания приближенного значения решения в используемом для интегрирования методе Штермера.

B качестве критерия точности при разгоне принята оценка первого отброшенного члена экстраполяционной формулы Адамса - Штермера, при этом все компоненты решения проверяются на точность по меpе погрешности, который заключается в следующем: если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы  P, называемой границей перехода, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.

Березин И.С., Жидков Н.П. Методы вычислений, т.2. Физматгиз, М., 1960.

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

    SUBROUTINE  DE34R (F, M, IORDER, XN, YN, DYN, HMIN,
                                           EPS, P, H, X, YX, DY, DF, RFN, RF, R, IERR) 
   Параметры
F - имя подпрограммы вычисления значений правой части системы. Первый оператор подпрограммы должен иметь вид:
 

SUBROUTINE  F (X, Y, DY, M)

Здесь: X, Y - значения независимой и зависимой переменных, соответственно; вычисленное значение правой части должно быть помещено в  DY. B случае системы уравнений, т.е. когда  M ≠ 1, параметры  Y и  DY представляют одномерные массивы длиной  M (тип параметров X, Y и DY: вещественный);
M - количество уравнений в системе (тип: целый);
IORDER - порядок точности без единицы того метода Штермера, для которого выполняется разгон и который будет использоваться при интегрировании данной системы уравнений (тип: целый);
      XN, YN -          DYN   начальные значения аргумента, решения и его производной; в случае системы уравнений (т.е.  M ≠ 1)  YN и  DYN представляют одномерные массивы длиной  M (тип: вещественный);
HMIN - минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
EPS - допустимая меpа погрешности, с которой тpебуется вычислить все компоненты решения (тип: вещественный);
P - граница перехода, используемая при оценке погрешности решения (тип: вещественный);
H - вещественная переменная, содержащая значение шага интегрирования. Если для этого значения шага точность при разгоне достигается, то именно он и реализуется на разгоне, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность;
X - вещественная переменная, значение которой на выходе из подпрограммы представляет первый (после начальной точки  XN) узел интегрирования  XN + H, в котоpом вычислены необходимые для интегрирования данной системы уравнений начальные значения;
YX, DY - одномерные вещественные массивы длиной  M, в которых запоминаются значения решения и его первой разности, вычисленные в узле  X, при этом погрешность решения и разности имеет (IОRDЕR + 1) - й порядок по  H;
DF - двумерный вещественный массив размера M * IORDER, в котоpом запоминаются значения правой части системы и ее разностей до (IОRDЕR - 1) - го порядка включительно, вычисленные в узле  X и умноженные на коэффициент   H2 / 12, при этом элемент  DF (I, 1) этого массива содержит значение правой части  I - го уравнения системы, а  DF (I, J + 1) - ее  J - ю разность, погрешность которой имеет (IОRDЕR + 1) - й порядок по  H;
         RFN -
       RF, R   		одномерные вещественные рабочие массивы длины  M;
IERR - целая переменная, значение которой в pезультате работы подпрограммы полагается равным 65, если на разгоне не может быть достигнута требуемая точность  EPS; в этом случае разгон можно начать сначала обращением к подпрограмме с новыми значениями паpаметpов  H и  HMIN. 
   Версии
DE34D - построение начальных значений с повышенной точностью при интегрировании системы обыкновенных дифференциальных уравнений второго порядка методом Штермера с контролем точности. При этом параметры  XN, YN, DYN, HMIN, EPS, P, H, X, YX, DY, DF, RFN, RF, R и параметры  X, Y и  DY в подпрограмме должны иметь тип DOUBLE PRECISION.
DE18R - назначение этой подпрограммы такое же, как и подпрограммы DE34R, но схема организации вычислений в ней несколько иная, чем в DE34R, а именно: значение решения вычисляется через значение некоторой промежуточной переменной, введение которой позволяет уменьшить вычислительную погрешность в приближенном значении решения. Все формальные параметры DE18R имеют тот же смысл, что и соответствующие им параметры DE34R с той только разницей, что в массив  DY заносится разделенная разность решения, вычисленная с порядком  IORDER по  H, а не конечная разность как в DE34R, а в массив  DF - конечные разности правой части системы, умноженные на коэффициент  H / 12 и вычисленные с порядком  IORDER по  H.
DE18D - то же, что и DE18R, но все вычисления ведутся с удвоенным числом значащих цифр; при этом параметры  XN, YN, DYN, HMIN, EPS, P, H, X, YX, DY, DF, RFN, RF, R и параметры  X, Y и  DY в подпрограмме  F должны иметь тип DOUBLE PRECISION. 
   Вызываемые подпрограммы
UTDE12 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE34R.
UTDE13 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE34D.
UTDE16 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE18R.
UTDE17 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE18D.
  Kpоме того, подпрограммы DE34R и DE34D используют рабочие подпрограммы DE28RS и DE28DS, соответственно. 
   Замечания по использованию
 

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

Значение праметра  H, задаваемое при обращении к подпрограмме, должно быть таким, чтобы узел  XN + (IORDER - 1) * H не выходил за конец интервала интегрирования.

Значение параметра  IORDER на входе в подпрограмму полагается равным 5. B дальнейшем предполагается расширить допустимое множество значений этого параметра.

При работе подпрограммы и ее версии значения параметров   M, IORDER, XN, YN, DYN, HMIN, EPS, P сохраняются.

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

       y1'' = - y1   , 
       y2'' = - y2   ,
       y1 (3/4 * π) = √2 / 2,        y2 (3/4 * π) = - √2 / 2   ,
       y1' (3/4 * π) = - √2 / 2,     y2' (3/4 * π) = - √2 / 2   .

      SUBROUTINE  F (X, Y, DY, M)
      DIMENSION  Y(2), DY(2)
      DY(1) = - Y(1)
      DY(2) = - Y(2)
      RETURN
      END

      DIMENSION  YN(2), DYN(2), Y(2), DY(2), DELTY(2), DF(2, 5), 
     *                        RFN(2), RF(2)
      EXTERNAL  F
      M = 2
      IORDER = 5
      XN = 0.75 * 3.14159265359
      YN(1) = SQRT(2.) / 2.
      YN(2) = - YN(1)
      DYN(1) = YN(2)
      DYN(2) = - YN(1)
      HMIN = 1.E - 10
      EPS = 0.0001
      P = 100.
      H = 0.01
      CALL  DE34R (F, M, IORDER, XN, YN, DYN, HMIN,
     *                        EPS, P, H, X, Y, DY, DF, RFN, RF, DELTY, IERR)

Результаты:

       IERR  =  0

       X  =  2.366194490 + 00

       Y(1)   =   7.000004762 - 01        Y(2)   =  - 7.141423761 - 01

       DY(1)  =  - 7.106305006 - 03        DY(2)  =  - 7.035594917 - 03

       DF(1, 1)  =  - 5.83333730151 - 06       DF(1, 2)  =     5.92192083634 - 08
       DF(1, 3)  =    5.89251453187 - 10        DF(1, 4)  =  - 5.86490578325 - 12
       DF(1, 5)  =  - 5.76760861293 - 14

       DF(2, 1)  =    5.95118646749 - 06       DF(2, 2)  =     5.86299576250 - 08
       DF(2, 3)  =  - 5.89250072347 - 10       DF(2, 4)  =  - 5.92374760355 - 12
       DF(2, 5)  =    6.01324545713 - 14