Текст подпрограммы и версий ( Фортран )
de29r.zip , de29d.zip 		Тексты тестовых примеров ( Фортран )
tde29r.zip , tde29d.zip
Текст подпрограммы и версий ( Си )
de29r_c.zip , de29d_c.zip 		Тексты тестовых примеров ( Си )
tde29r_c.zip , tde29d_c.zip
Текст подпрограммы и версий ( Паскаль )
de29r_p.zip , de29e_p.zip 		Тексты тестовых примеров ( Паскаль )
tde29r_p.zip , tde29e_p.zip

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

Назначение

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

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

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

          Y ' = F (X, Y) ,

          Y = ( y1, ..., yM ) ,
          F = ( f1 (X, y1,..., yM), ... , fM (X, y1,..., yM) )
 с начальными условиями, заданными в точке XN:
          Y(XN) = YN ,     YN = ( y10,..., yM0 ) ,

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

Решение вычисляется в одной точке XK, которая является концом интервала интегрирования. Для всех компонент решения осуществляется контроль точности по меpе погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P (называемой границей перехода), то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.

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

Беленький В.З. Стандартная программа для интегрирования системы дифференциальных уравнений методом Адамса. Сб. "Вычислительные методы и программирование", вып. 3, Изд-во МГУ, 1965

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

    SUBROUTINE  DE29R (F, M, XN, YN, XK, HMIN, HMAX, EPS, P, H,
                                           Y, DELTY, DF, RF, YP, RY, RFN, IERR) 
   Параметры
F - имя подпрограммы вычисления значений правой части дифференциального уравнения. Первый оператоp подпрограммы должен иметь вид:
SUBROUTINE  F (X, Y , DY, M).
Здесь: X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY. B случае системы уравнений, т.е. когда M ≠ 1, параметры Y и DY представляют массивы длины M (тип параметров X, Y и DY: вещественный); имя подпрограммы вычисления значений правой части системы.
M - количество уравнений в системе (тип: целый);
XN, YN - начальные значения аргумента и решения; в случае системы уравнений (т.е. M ≠ 1) YN представляет массив длины M (тип: вещественный);
XK - значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования); XK может быть больше, меньше или pавно XN (тип: вещественный);
HMIN - минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
HMAX - максимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
EPS - допустимая меpа погрешности, с которой тpебуется вычислить все компоненты решения (тип: вещественный);
P - граница перехода, используемая при оценке меры погрешности решения (тип: вещественный);
H - вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если X > XN, отрицательным, если XK < XN, или без такого учета в виде абсолютной величины;
Y - искомое решение задачи Коши, вычисленное подпрограммой для значения аргумента XK; для системы уравнений (когда M ≠ 1) задается одномерным массивом длины M; в случае совпадения значений параметров XK и XN значение Y полагается равным начальному значению YN (тип: вещественный);
         DELTY -
         RF, YP  
         RFN   		одномерные вещественные рабочие массивы длины M;
DF, RY - двумерные вещественные рабочие массивы размеpа M*5;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом:
           IERR=65 -
         IERR=66   		когда какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS при заданных начальном шаге H и его минимальном значении HMIN; при этом IERR = 66 указывает, что требуемая точность не может быть достигнута при разгоне, а IERR = 65 - после разгона; при желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров H и HMIN. 
   Версии
DE29D - вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Адамса с повышенной точностью. При этом параметры XN, YN, XK, HMIN, HMAX, EPS, P, H, Y, DELTY, DF, RF, YP, RY, RFN и параметры X, Y и DY в подпрограмме F должны иметь тип DOUBLE PRECISION. 
   Вызываемые подпрограммы
DE28R - выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом Адамса с контролем точности.
DE28D - выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом Адамса с повышенной точностью.
      UTDE12 -       UTDE13   подпрограммы выдачи диагностических сообщений.
  Подпрограммы DE28R и UTDE12 вызываются при работе подпрограммы DE29R, а подпрограммы DE28D и UTDE13 - при работе подпрограммы DE29D. 
   Замечания по использованию
 

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

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

Если после работы подпрограммы нет необходимости иметь начальное значение решения YN, то параметры YN и Y при обращении к ней можно совместить. При этом надо иметь в виду, что в случае аварийного выхода из подпрограммы, т.е. со значением IERR = 65 или 66, значение параметра YN будет испорчено.

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

           y1'  =   y2
           y2'  =  -y1 ,     3/4 π ≤ x ≤ π

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

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

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

          DIMENSION  YN(2), Y(2), DELTY(2), DF(10), RF(2),
         *                        YP(2), RY(10), RFN(2)
          EXTERNAL  F
          M = 2
          XN = 0.75*3.14159265359
          YN(1) = SQRT(2.)/2.
          YN(2) = -YN(1)
          XK = 3.14159265359
          HMIN = 1.E-18
          HMAX = 0.1
          EPS = 0.00001
          P = 100.
          H = 0.01
          CALL  DE29R (F, M, XN, YN, XK, HMIN, HMAX, EPS, P, H, Y,
         *                        DELTY, DF, RF, YP, RY, RFN, IERR)

Результаты:

          IERR  =  0
          
          Y(1)  =   3.010855175-07
          Y(2)  =  -1.000001503+00