Текст подпрограммы и версий ( Фортран )
de35r.zip , de35d.zip
Тексты тестовых примеров ( Фортран )
tde35r.zip , tde35d.zip
Текст подпрограммы и версий ( Си )
de35r_c.zip , de35d_c.zip
Тексты тестовых примеров ( Си )
tde35r_c.zip , tde35d_c.zip
Текст подпрограммы и версий ( Паскаль )
de35r_p.zip , de35e_p.zip
Тексты тестовых примеров ( Паскаль )
tde35r_p.zip , tde35e_p.zip

Подпрограмма:  DE35R

Назначение

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

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

Для системы  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, а другая является экстраполяционной формулой, применяемой для предсказания приближенного значения решения в используемом для интегрирования методе Штермера.

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

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

    SUBROUTINE  DE35R (F, M, IORDER, XN, YN, DYN, H,
                                            X, YX, DY, DF, RFN, RF, R) 

Параметры

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 (тип: вещественный);
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.

Версии

DE35D - построение начальных значений при интегрировании системы обыкновенных дифференциальных уpавнений второго порядка методом Штермера без контроля точности, при этом все вычисления с действительными числами выполняются с удвоенным числом значащих цифр; в этом случае параметры   XN, YN, DYN, H, X, YX, DY, DF, N, RF, R и параметры  X, Y и  DY в подпрограмме  F должны иметь тип DOUBLE PRECISION.

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

  Подпрограммы DE35R и DE35D используют рабочие подпрограммы DE28RS и DE28DS, соответственно.

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

 

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

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

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

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

      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)
      H = 0.01
      CALL  DE35R (F, M, IORDER, XN, YN, DYN, H, X,
     *                          Y, DY, DF, RFN, RF, DELTY)

Результаты:

       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