Документ взят из кэша поисковой машины. Адрес оригинального документа : http://num-anal.srcc.msu.ru/lib_na/cat/de/de40r.htm
Дата изменения: Tue Dec 1 11:54:32 2015
Дата индексирования: Sun Apr 10 00:31:11 2016
Кодировка: Windows-1251
БЧА НИВЦ МГУ. DE40R. Задачи Коши для системы уравнений второго порядка
Текст подпрограммы и версий ( Фортран )
de40r.zip , de40d.zip
Тексты тестовых примеров ( Фортран )
tde40r.zip , tde40d.zip
Текст подпрограммы и версий ( Си )
de40r_c.zip , de40d_c.zip
Тексты тестовых примеров ( Си )
tde40r_c.zip , tde40d_c.zip
Текст подпрограммы и версий ( Паскаль )
de40r_p.zip , de40e_p.zip
Тексты тестовых примеров ( Паскаль )
tde40r_p.zip , tde40e_p.zip

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

Назначение

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

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

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

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

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

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

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

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

Жоголев E.A. Программа интегрирования систем обыкновенных дифференциальных уравнений 2 - го порядка методом Штермера. Сб. "Вычислительные методы и программирование", вып. 1 , Изд - во МГУ, 1962.

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

    SUBROUTINE  DE40R (F, M, JSTART, YX, X, DY, H, DELTY, DF,
                                             RFN, RF, YP, DYP, IERR) 

Параметры

F - имя подпрограммы вычисления значений правой части системы дифференциальных уpавнений. Первый оператор подпрограммы должен иметь вид:
SUBROUTINE  F (X, Y, DY, M).
Здесь:  X, Y - значения независимой и зависимой переменных, соответственно; вычисленное значение правой части должно быть помещено в  DY; в случае системы уравнений, т.е. когда  M ≠ 1, параметры  Y и  DY представляют одномерные массивы длиной  M (тип параметров  X, Y, DY: вещественный);
M - количество уравнений в системе (тип: целый);
JSTART - целый указатель режима использования подпрограммы, имеющий следующие значения:
0 - первое обращение к подпрограмме должно быть выполнено с нулевым значением   JSTART;
+1 - выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных и шага интегрирования, заданных параметрами  X, YX и  H, соответственно;
-1 - повторить последний шаг интегрирования с новым значением параметра  H;
  на выходе из подпрограмы  JSTART полагается равным + 1;
X, YX - заданные вещественные значения аргумента и решения в нем; в результате работы подпрограммы в  X получается новое значение аргумента, а в  YX - соответствующее значение решения; в случае системы уравнений, т.е. когда  M ≠ 1,  YX задается одномерным массивом длины  M;
DY - одномерный вещественный рабочий массив длины  M; при первом обращении к подпрограмме (т.е. со значением  JSTART = 0) в этом массиве задается значение производной pешения в начальном узле интервала интегрирования;
H - вещественное значение шага интегрирования;
       DELTY -
      RFN, RF  
      YP, DYP  
одномерные вещественные рабочие массивы длины  M;
DF - двумерный вещественный рабочий массив размеpа  M * 5, в котоpом запоминаются значения правой части системы и ее разностей до четвертого порядка включительно, умноженные на  H 2 / 12;
IERR - целая переменная, значение которой в pезультате работы подпрограммы полагается равным 65, если при заданном значении шага интегрирования не может быть достигнута сходимость численного решения к точному значению решения; в этом случае интегрирование системы можно повторить обращением к подпрограмме с новым значением параметра  H и со значением  JSTART = + 1 .

Версии

DE40D - выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений второго порядка методом Штермера без контроля точности, при этом все вычисления над действительными числами ведутся с удвоенным числом значащих цифр. В этом случае параметры   YX, X, DY, H, DELTY, DF, RFN, RF, YP, DYP и параметры  X, Y и  DY в подпрограмме  F должны иметь тип DOUBLE PRECISION.

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

DE35R - построение начальных значений при интегрировании системы обыкновенных дифференциальных уравнений второго порядка методом Штермера без контроля точности; вызывается при работе подпрограммы DE40R.
DE35D - построение начальных значений при интегрировании системы обыкновенных дифференциальных уравнений второго порядка методом Штермера без контроля точности, при этом все вычисления над действительными числами ведутся с удвоенным числом значащих цифр; вызывается при работе подпрограммы DE40D.
UTDE12 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE40R.
UTDE13 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE40D.
  Kpоме того, DE40R и DE40D используют рабочие подпрограммы DE42RP, DE28RS и DE42DP, DE28DS, соответственно.

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

 

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

При многократном использовании подпрограммы и ее версии для вычисления решения на отрезке значения параметров   M, YX, X, DY, DF, YP, DYP не должны изменяться в вызывающей программе между последовательными обращениями к ним.

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

     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  Y(2), DY(2), DELTY(2), DF(10), RFN(2), RF(2), YP(2), 
     *                       DYP (2)
      EXTERNAL  F
      M = 2
      JSTART = 0
      Y(1) = SQRT(2.) / 2.
      Y(2) = - Y(1)
      X = 0.75 * 3.14159265359
      DY(1) = Y(2)
      DY(2) = - Y(1)
      H = 0.01
      IH = 0
   6 CONTINUE
      CALL  DE40R (F, M, JSTART, Y, X, DY, H, DELTY, DF, RFN, RF,
     *                        YP, DYP, IERR)
      IH = IH + 1
      Y1 = SIN(X)
      Y2 = COS(X)
C    ПEЧATЬ TOЧHOГO И ПPИБЛИЖEHHOГO ЗHAЧEHИЙ PEШEHИЯ 
      PRINT 1, X, Y, H, Y1, Y2
      IF (IH .EQ. 1) GO TO 6
      IF (IH .EQ. 5) GO TO 20
      JSTART = - 1
      IF (IH .EQ. 3) GO TO 7
      IF (IH .EQ. 4) GO TO 8
      H = - 0.005
      GO TO 6
   7 H = 0.02
      GO TO 6
   8 H = 0.01
      GO TO 6
 20 STOP
   1 FORMAT (3E18.9)

Результаты:

после первого обращения к подпрограмме -

                  X                             Y(1)                                Y(2)
       2.366194490 + 00      7.000004762 - 01        - 7.141423761 - 01

                  H                             Y1                                  Y2
       1.000000000 - 02      7.000004762 - 01        - 7.141423761 - 01

после второго обращения к подпрограмме -
 
                  X                             Y(1)                                Y(2)
       2.376194490 + 00      6.928241717 - 01        - 7.211065574 - 01

                  H                             Y1                                  Y2
       1.000000000 - 02      6.928241717 - 01        - 7.211065574 - 01

после третьего обращения к подпрограмме -

                  X                             Y(1)                                Y(2)
       2.361194490 + 00      7.035624232 - 01        - 7.106334615 - 01

                  H                             Y1                                  Y2
     - 5.000000000 - 03      7.035624232 - 01        - 7.106334615 - 01

после четвертого обращения к подпрограмме -
  
                  X                             Y(1)                                 Y(2)
       2.386194490 + 00      6.855785854 - 01        - 7.279986286 - 01

                  H                             Y1                                   Y2
       2.000000000 - 02      6.855785854 - 01        - 7.279986286 - 01

после пятого обращения к подпрограмме -

                  X                             Y(1)                                 Y(2)
       2.376194490 + 00      6.928241717 - 01        - 7.211065574 - 01

                  H                             Y1                                   Y2
       1.000000000 - 02      6.928241717 - 01        - 7.211065574 - 01