Текст подпрограммы и версий ( Фортран )
de51r.zip , de51d.zip 		Тексты тестовых примеров ( Фортран )
tde51r.zip , tde51d.zip
Текст подпрограммы и версий ( Си )
de51r_c.zip , de51d_c.zip 		Тексты тестовых примеров ( Си )
tde51r_c.zip , tde51d_c.zip
Текст подпрограммы и версий ( Паскаль )
de51r_p.zip , de51e_p.zip 		Тексты тестовых примеров ( Паскаль )
tde51r_p.zip , tde51e_p.zip

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

Назначение

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

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

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

(1)       y '' + f(x, y, y ') y ' + g(x, y, y ') y = r(x, y, y ')  ,     xN≤ x ≤ xK

с заданными на концах отрезка интегрирования  [ xN, xK ] граничными условиями вида 

(2)       aN y '(xN) + bN y(xN) = CN   ,

(3)       aK y '(xK) + bK y(xK) = CK   .

Решение вычисляется на заданной на отрезке  [ xN, xK ] равномерной сетке с заданной абсолютной погрешностью  EPS. Используется метод линеаризации исходного уравнения (1) с последующим применением метода конечных разностей для решения линейного уравнения.

L.Fox, Proc. Roy. Soc. A 190, 31 - 59, 1947.

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

    SUBROUTINE  DE51R (F, NX, XN, XK, CN, CK, EPS, IMAX, IDER,
                                            Y, IK, DELTY, RAB, IERR) 
   Параметры
F - имя подпрограммы вычисления значений коэффициентов  f (x, y, y'),  g (x, y, y'),  r (x, y, y') уравнения (1). Первый оператор подпрограммы должен иметь вид:
 

SUBROUTINE  F (CF, CG, CR, X, Y, I, RAB)
Здесь:

X, Y - представляют первые два аргумента функций  f (x, y, y'),  g (x, y, y') и   r (x, y, y');
RAB - одномерный массив длины 11 * NX (NX - число узлов сетки);
I - указывает компоненту  RAB(I) массива   RAB, которая представляет третий аpгумент функций  f (x, y, y'),  g (x, y, y'),   r (x, y, y');
CF, CG, CR - представляют вычисленные подпрограммой  F значения функций  f (x, y, y'),  g (x, y, y') и  r (x, y, y');
  параметры  CF, CG, CR, X, Y и  RAB имеют тип REAL, а  I - INTEGER;
NX - число узлов равномерной сетки на отрезке   [xN, xK], (включая  xN и  xK), в которых требуется вычислить решение задачи;  NX ≥ 3 (тип: целый);
XN, XK - концы отрезка интегрирования, в которых заданы граничные условия (2) и (3), соответственно;  XN < XK (тип: вещественный);
CN, CK - одномерные вещественные массивы длины 3, в которых помещаются коэффициенты граничных условий (2) и (3), соответственно, т.е.  CN (1) = aN,  CN (2) = bN,  CN (3) = CN,  CK (1) = aK, CK (2) = bK,  CK (3) = CK. Коэффициенты должны удовлетворять условию
a2N + b2N ≠ 0   и   a2K + b2K ≠ 0;
EPS - абсолютная погрешность, с которой требуется вычислить решение краевой задачи;  EPS > 0 (тип: вещественный);
IMAX - указывает максимальное допустимое число линеаризаций исходного уравнения (1) и число итераций, используемых при решении получающегося после линеаризации уравнения: число линеаризаций pавно абсолютному значению  | IMAX | параметра  IMAX, а число итераций pавно  IMAX, если  IMAX > 0 и нулю, если  IMAX < 0;  IMAX ≠ 0. Если  IMAX < 0, то в результате работы подпрограммы значение  IMAX заменяется на противоположное; в этом случае  IMAX зажается переменной или элементом массива (тип: целый);
IDER - указывает на зависимость коэффициентов  f (x, y, y'),  g (x, y, y'),  r (x, y, y') уравнения (1) от производной  y': если хотя бы один из этих коэффициентов зависит от  y', то  IDER > 0, иначе  IDER = 0 (тип: целый);
Y - одномерный вещественный массив длины  NX, который должен содержать при обращении к подпрограмме значения начального приближения решения краевой задачи (1), (2), (3) во всех узлах сетки. В результате работы подпрограммы в этом массиве получаются значения конечного приближения к решению (эти значения помещаются в  Y в любом случае, даже и тогда, когда они не могут быть вычислены с требуемой точностью);
IK - целая переменная, содержащая число линеаризаций уравнения (1), выполненных в действительности подпрограммой для достижения требуемой точности  EPS решения. В случае, когда данная точность решения нелинейного уравнения не достигается за  IMAX линеаризаций или когда решение соответствующего линейного уравнения не может быть вычислено с точностью  EPS за  IMAX итераций, то значение параметра  IK в результате работы подпрограммы полагается равным  - 1;
DELTY - одномерный вещественный массив длины  NX, значение которого в результате работы подпрограммы полагается равным разности между начальным и конечным приближениями решения того линейного уравнения, котоpое получается после выполнения самой последней линеаризации исходного уpавнения (1), причем эта разность вычисляется во всех узлах сетки. В случае, когда  IMAX < 0, все элементы массива  DELTY в результате работы подпрограммы полагаются равными 0;
RAB - одномерный вещественный рабочий массив длины 11 * NX;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом:
IERR=65 - если в граничном условии (2) коэффициенты  aN = bN = 0;
IERR=66 - если в граничном условии (3) коэффициенты  aK = bK = 0;
IERR=67 - если  IMAX = 0;
IERR=68 - если  IDER < 0;
IERR = 1 - если заданная точность решения нелинейного уравнения не достигается за  | IMAX | линеаризаций или если решение соответствующего линейного уpавнения не может быть вычислено с точностью  EPS за  IMAX итераций;
IERR = 3 - если заданная точность решения нелинейного уравнения не достигается за  | IMAX | линеаризаций или если решение соответствующего линейного уpавнения не может быть вычислено с точностью  EPS за  IMAX итераций, при этом при аппроксимации исходного уравнения линейным уравнением использовались производные, вычисленные с меньшей точностью, чем это необходимо;
  В каждом из выше указаных случаев решение задачи прекращается;
IERR = 2 - если при аппроксимации исходного уpавнения линейным уравнением использовались производные, вычисленные с меньшей точностью, чем это необходимо. При  IERR = 1, 2, 3 полученное приближенное решение запоминается в массиве  Y; при желании интегрирование можно повторить обращением к подпрограмме с новыми значениями параметров  NX, Y или  IMAX. 
   Версии
DE51D - вычисление решения двухточечной краевой задачи для нелинейного дифференциального уравнения второго порядка с повышенной точностью. При этом параметры  XN, XK, CN, CK, EPS, Y, DELTY, RAB должны иметь тип DOUBLE PRECISION. 
   Вызываемые подпрограммы
DE50R - вычисление решения двухточечной краевой задачи для линейного дифференциального уравнения второго порядка методом конечных разностей.
DE50D - вычисление решения двухточечной краевой задачи для линейного дифференциального уравнения второго порядка методом конечных разностей с повышенной точностью.
ID10R - подпрограмма вычисления первых и вторых производных с помощью центральных разностей.
ID10D - подпрограмма вычисления первых и вторых производных с помощью центральных разностей с двойной точностью.
  Подпрограммы DE50R и ID10R вызываются подпрограммой DE51R, а подпрограммы DE50D и ID10D - подпрограммой DE51D.
UTDE10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE51R.
UTDE11 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE51D. 
   Замечания по использованию
 

Число узлов сетки  NX, задаваемое пользователем при обращении к подпрограмме, должно быть достаточно большим для того, чтобы обеспечить сходимость конечно - разностных аппроксимаций, используемых для решения линейных уравнений, и последовательных линеаризаций исходного уравнения (1).

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

Пусть  hi ,  i = 1, ..., NX, представляет разность между текущим приближением к решению и приближением, полученным в результате предыдущей линеаризации уравнения (1). Тогда текущее приближение принимается за окончательное решение, если для всех  i (1 ≤ i ≤ NX) выполняется условие  | hi | ≤ EPS .

Значения параметров  NX, XN, XK, CN, CK, EPS, IMAX, в результате работы подпрограммы сохраняются.

Kpоме того, DE50R и DE50D используют рабочие подпрограммы ID10RP и ID10DP, соответственно.

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

Решается краевая задача:

      y '' + ( xy' ) y' + e -y y = ( xy' + 2e-y - 2 ) cos x - ( 1 + 2xy' - e-y ) sin x 
      2y - y' = 0         при x = 0
      y + y' = 1          при x = π / 2

(точное решение задачи  y = sin x + 2cos x).

      SUBROUTINE  FUNCTS (F, G, R, X, Y, I, RAB)
      DIMENSION  RAB(1)
      F = X * RAB(I)
      G = EXP(- Y)
      R =  (X * RAB(I) - 2. + 2. * EXP(- Y)) * COS(X) - (1. + 2. * X * RAB(I) -
     *        EXP(- Y)) * SIN(X)
      RETURN
      END

      DIMENSION  CN(3), CK(3), YX(150), DELTY(150), RAB(1650)
      EXTERNAL  FUNCTS
      NX = 150
      XN = 0.
      XK = 3.14159265359
      XK = XK / 2
      CN(1) = 2.
      CN(2) = - 1.
      CN(3) = 0.
      CK(1) = 1.
      CK(2) = 1.
      CK(3) = - 1.
      IMAX = 4
      EPS = 0.00001
      IDER = 1
      DO 1  I = 1, NX
   1 YX(I) = 1.5
      CALL  DE51R (FUNCTS, NX, XN, XK, CN, CK, EPS, IMAX, IDER,
     *                        YX, IK, DELTY, RAB, IERR)

Результаты:      (приводятся только для первых 10 узлов)

      YX(1)   =  1.999999086
      YX(2)   =  2.010430009
      YX(3)   =  2.020637497
      YX(4)   =  2.030620414
      YX(5)   =  2.040377652
      YX(6)   =  2.049908127
      YX(7)   =  2.059210778
      YX(8)   =  2.068284572
      YX(9)   =  2.077128500
      YX(10) =  2.085741581

      IK  =  3
      IERR  =  2