Текст подпрограммы и версий ( Фортран )
de54r.zip , de54d.zip , de56r.zip , de56d.zip 		Тексты тестовых примеров ( Фортран )
tde54r.zip , tde54d.zip , tde56r.zip , tde56d.zip
Текст подпрограммы и версий ( Си )
de54r_c.zip , de54d_c.zip , de56r_c.zip , de56d_c.zip 		Тексты тестовых примеров ( Си )
tde54r_c.zip , tde54d_c.zip , tde56r_c.zip , tde56d_c.zip
Текст подпрограммы и версий ( Паскаль )
de54r_p.zip , de54e_p.zip , de56r_p.zip , de56e_p.zip 		Тексты тестовых примеров ( Паскаль )
tde54r_p.zip , tde54e_p.zip , tde56r_p.zip , tde56e_p.zip

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

Назначение

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

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

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

(1)       y '' + f(x) y ' + g(x) y = r(x)

c заданными на концах  xN и  xК отрезка интегрирования граничными условиями вида: 

(2)       aN y '(xN) + bN y(xN) = cN   ,      a2N + b2N ≠ 0

(3)       aK y '(xK) + bK y(xK) = cK   ,      a2K + b2K ≠ 0

- методом прогонки А.А.Абрамова. Решение вычисляется на заданной на отрезке интегрирования равномерной сетке с заданной мерой погрешности  EPS. Контроль точности по меpе погрешности заключается в следующем: если численное решение по абсолютной величине не меньше некоторой заданной константы  P, то контроль точности ведется по относительной погрешности, иначе - по абсолютной.

А.А.Абрамов. Вариант метода прогонки. Ж. вычисл. матем. и матем. физ., 1961, 1, N 2, 349 - 351.

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

    SUBROUTINE  DE54R (FF, FG, FR, XN, XK, CN, CK, NX, XO, HX,
                                           BULDY, HMIN, EPS, P, H, Y, DY, NX1, IERR) 
   Параметры
         FF, FG -
         FR   		имена вещественных подпрограмм - функций вычисления значений коэффициентов уравнения (1)  f (x), g (x) и  r (x), соответственно. Эти подпрограммы - функции должны иметь по одному паpаметpу, представляющему вещественное значение независимой переменной  x;
XN, XK - концы отрезка интегрирования, в которых заданы граничные условия (2) и (3), соответственно;  XK > XN или  XK < XN (тип: вещественный);
CN, CK - одномерные вещественные массивы длины 3, в которых помещаются коэффициенты граничных условий (2) и (3), т.е.
CN (1) = aN , CN (2) = bN , CN (3) = cN ,
CK (1) = aK , CK (2) = bK , CK (3) = cK ;
коэффициенты граничных условий должны удовлетворять неpавенствам:
a2N + b2N ≠ 0  и  a2K + b2K ≠ 0 ;
NX - число узлов равномерной сетки на отрезке интегрирования, в которых требуется вычислить решение задачи;  NX ≥ 1 (тип: целый);
XO - начальный узел равномерной сетки;  XO должен находиться между  XN и  XK или совпадать с каким-нибудь концом (тип: вещественный);
HX - шаг равномерной сетки;  HX > 0, если  XK > XN и  HX < 0, если  XK < XN (тип: вещественный);
BULDY - логическая переменная или константа, принимающая при входе в подпрограмму значение .TRUE., если вместе с решением задачи требуется вычислить и его производную, и .FALSE. - в противном случае, т.е. когда необходимо иметь только решение задачи;
HMIN - минимальное значение абсолютной величины шага интегрирования, котоpое разрешается использовать при численном решении воспомагательных задач Коши, к которым сводится данная краевая задача (тип: вещественный);
EPS - допустимая меpа погрешности, с которой тpебуется вычислить решение краевой задачи (тип: вещественный);
P - граница перехода, используемая при оценке меры погрешности решения (тип: вещественный);
H - вещественная переменная, содержащая начальное значение шага интегрирования, используемое при численном решении воспомагательных задач Коши;
NX1 - целая переменная, указывающая на выходе из подпрограммы число узлов сетки, в которых вычислено решение данной краевой задачи;  NX1 ≤ NX;
Y, DY - одномерные вещественные массивы длины не более  NX1, в которых запоминаются значения pешения и его производной, вычисленные в  NX1 узлах сетки, при этом элементы  Y (I) и  DY (I) содержат значения решения и производной в узле  XO + (I - 1) * H, причем производная вычисляется только тогда, когда  BULDY=.TRUE. (таким образом, предполагается, что узлы сетки занумерованы в направлении от XN до  XK);
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом:
IERR=1 - когда последний узел сетки  XO + (N - 1) * H выходит за конец  XK интервала интегрирования, а предпоследний узел  XO + (N - 2) * H является внутренним узлом, т.е. лежит внутри интервала интегрирования;
IERR=3 - когда несколько узлов сетки (больше одного) выходят за конец  XK интервала интегрирования, а самый последний узел, принадлежащий отрезку интегрирования, является внутренним узлом;
IERR=2 - когда один из узлов сетки совпадает с концом  XK интервала интегрирования, а все последующие узлы не попадают в интервал интегрирования;
  во всех указанных выше случаях решение и его производная вычисляются только в тех узлах сетки, которые принадлежат отрезку интегрирования, а в случае  IERR = 1 и  IERR = 3 они вычисляются также в конце  XK отрезка и запоминаются в  Y (NX1) и  DY (NX1), соответственно;
IERR=65 - когда в граничном условии (2) на конце  XN коэффициенты  aN = bN = 0;
IERR=66 - когда в граничном условии (3) на конце  XK коэффициенты  aK = bK = 0;
IERR=67 - когда начальный узел  XO сетки не принадлежит отрезку интегрирования;
IERR=68 и IERR=70 - когда решение краевой задачи не может быть вычислено с требуемой точностью  EPS при заданном минимальном значении  HMIN шага интегрирования, при этом  IERR = 68 указывает, что точность не достигается при прямой прогонке Абрамова, а IERR=70 - при обратной прогонке;
  при  IERR = 65, 66, 67, 68, 70 интегрирование уравнения прекращается; в случае  IERR = 68 и  IERR = 70 интегрирование можно повторить обращением к подпрограмме с новыми значениями  H и  HMIN;
IERR=69 - когда данная краевая задача не может быть решена методом прогонки А.А.Абрамова. 
   Версии
DE54D - вычисление решения двухточечной краевой задачи для линейного дифференциального уравнения второго порядка на равномерной сетке методом прогонки А.А.Абрамова с повышенной точностью. При этом параметры  FF, FG, FR, XN, XK, CN, CK, XO, HX, HMIN, EPS, P, H, Y, DY и формальный параметp в подпрограммах - функциях FF, FG и  FR должны иметь тип DOUBLE PRECISION.
DE56R - вычисление решения двухточечной краевой задачи для линейного дифференциального уравнения второго порядка на неравномерной сетке методом прогонки А.А.Абрамова с контролем точности. Первый оператор подпрограммы имеет вид:
 

SUBROUTINE  DE56R (FF, FG, FR, XN, XK, CN, CK, NX, X, BULDY, HMIN, EPS, P, H, Y, DY, NX1, IERR)

Здесь:  X - одномерный вещественный массив длины  NX, представляющий узлы неравномерной сетки, в которых требуется вычислить решение задачи, при этом узлы сетки располагаются в  X в направлении от  XN до  XK.

Остальные параметры подпрограммы DE56R имеют тот же смысл, что и одноименные параметры подпрограммы DE54R. При этом выход из подпрограммы DE56R со значением  IERR = 67 указывает на то, что первый узел сетки  X (1) не принадлежит интервалу интегрирования.
DE56D - вычисление решения двухточечной краевой задачи для линейного дифференциального уравнения второго порядка на неравномерной сетке методом прогонки А.А.Абрамова с повышенной точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE56R; при этом параметры  FF, FG, FR, XN, XK, CN, CK, X, HMIN, EPS, P, H, Y, DY и формальный параметр в подпрограммах - функциях  FF, FG, FR должны иметь тип DOUBLE PRECISION. 
   Вызываемые подпрограммы
UTDE12 - подпрограмма выдачи диагностических сообщений при работе подпрограмм DE54R и DE56R.
UTDE13 - подпрограмма выдачи диагностических сообщений при работе подпрограмм DE54D и DE56D.
  Kpоме того, подпрограммы DE54R, DE56R и подпрограммы DE54D, DE56D используют рабочие подпрограммы DE54RU и DE54DU, соответственно. 
   Замечания по использованию
 

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

Значения параметров  XN, XK, CN, CK, NX, XO, HX, BULDY, HMIN, EPS, P, X сохраняются.

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

    y '' + x2 y ' + xy = 10e- x*x (2x2 - 1) (2 - X) ,    - 1≤  x≤ 1

    y ' - 2y = 0   при  x = - 1 ,
    y ' + 2y = 0   при  x = 1

( точное решение задачи  y (x) = 10e- x*x )

      FUNCTION  F (X)
      F = X * X
      RETURN
      END

      FUNCTION  G (X)
      G = X
      RETURN
      END

      FUNCTION  R (X)
      X2 = X * X
      R = 10. * EXP(- X2) * (2. * X2 - 1.) * (2. - X)
      RETURN
      END

      DIMENSION  CN(3), CK(3), Y(11), DY(11)
      LOGICAL  BULDY 
      DATA  CN(1), CN(2), CN(3) /1., -2., 0./
      DATA  CK(1), CK(2), CK(3) /1., 2., 0./
      EXTERNAL  F, G, R
      XN = - 1.
      XK = 1.
      NX = 11 
      XO = XN
      HX = 2. / 199.
      BULDY = .TRUE.
      HMIN = 1.E - 10
      EPS = 0.00001
      P = 11.
      H = 0.01
      CALL  DE54R (F, G, R, XN, XK, CN, CK, NX, XO, HX, BULDY,
     *                         HMIN, EPS, P, H, Y, DY, NX1, IERR)

Результаты:

      IERR  =  0
      NX1   =  11

                    Y                                     DY
      4.45260617314 + 00       - 8.01021606015 + 00
      4.37238323319 + 00       - 7.95378250048 + 00
      4.29273839036 + 00       - 7.89518713214 + 00
      4.21369299696 + 00       - 7.83450451556 + 00
      4.13526765439 + 00       - 7.77180951660 + 00
      4.05748221040 + 00       - 7.70717722605 + 00
      3.98035575729 + 00       - 7.64068288024 + 00
      3.90390663078 + 00       - 7.57240178309 + 00
      3.82815240978 + 00       - 7.50240922954 + 00
      3.75310991678 + 00       - 7.43078043029 + 00
      3.67879521918 + 00       - 7.35759043838 + 00