БЧА НИВЦ МГУ. DE98D. Один шаг интегрирования системы с правой частью, зависящей от произврдной, методом рядов Чебышёва

Текст подпрограммы и версий ( Фортран ) de98d.zip	Тексты тестовых примеров ( Фортран ) tde98d1.zip , tde98d2.zip , tde98d3.zip , tde98d4.zip
Текст подпрограммы и версий ( Си ) de98d_c.zip	Тексты тестовых примеров ( Си ) tde98d1_c.zip , tde98d2_c.zip , tde98d3_c.zip , tde98d4_c.zip
Текст подпрограммы и версий ( Паскаль ) de98e_p.zip	Тексты тестовых примеров ( Паскаль ) tde98e1_p.zip , tde98e2_p.zip , tde98e3_p.zip , tde98e4_p.zip

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

Назначение

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

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

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

(1)       Y '' = F(X, Y, Y'),
          Y  = (y₁, ... , y_M),  Y'  = (y'₁, ... , y'_M),   
          F = (f₁ (X, y₁, ... , y_M, y'₁, ... , y'_M), ... , f_M (X, y₁, ... , y_M, y'₁, ... , y'_M)),    XN ≤ X ≤ XK,

с начальными условиями, заданными в точке XN :

(2)       Y(XN) = YN ,     YN = (y₁₀, ... , y_M0),     Y'(XN) = DYN,     DYN = (y'₁₀, ... , y'_M0),

при условии, что правая часть системы (1) имеет непрерывные ограниченные частные производные по переменным X, Y, Y'. Предполагается, что на отрезке [XN, XK] задача (1),(2) имеет единственное решение. Тогда решение задачи Коши Y(XN + αΔ) и его производные Y'(XN + αΔ), Y''(XN + αΔ) = F(XN + αΔ, Y(XN + αΔ), Y'(XN + αΔ)) = Φ(α) разлагаются на промежутке интегрирования [XN, XK] в равномерно сходящиеся ряды по смещенным многочленам Чебышёва первого рода

                                  _∞                                                                   ₁
(3)     Y(XN + αΔ) = ∑^'a_i^*[Y]T_i^*(α),     0 ≤ α ≤ 1,     a_i^*[Y] = 2/π ∫ Y(XN + αΔ)T_i^*(α) / √(α(1 - α))dα,
                                 ⁱ⁼⁰                                                                  ⁰
 
                                   _∞                                                                    ₁
(4)     Y'(XN + αΔ) = ∑^'a_i^*[Y']T_i^*(α),     0 ≤ α ≤ 1,     a_i^*[Y'] = 2/π ∫ Y'(XN + αΔ)T_i^*(α) / √(α(1 - α))dα,
                                  ⁱ⁼⁰                                                                   ⁰

                      _∞                                                                   ₁
(5)     Φ(α) = ∑^'a_i^*[Φ]T_i^*(α),     0 ≤ α ≤ 1,     a_i^*[Φ] = 2/π ∫ Φ(α)T_i^*(α) / √(α(1 - α))dα.
                     ⁱ⁼⁰                                                                  ⁰

Здесь: штрих у знака суммы означает, что слагаемое с индексом 0 берется с дополнительным множителем 1/2; T_i^* (α) - смещенный многочлен Чебышёва первого рода на [0, 1]: T_i^* (α) = T_i (2α-1) ; T_i (t) - многочлен Чебышёва первого рода i-го порядка на [-1, 1]; Δ = XK - XN. Если ряд Чебышёва (3) (и ряды (4), (5)) является быстросходящимся, то его сумма на [XN, XK] (и суммы рядов (4), (5)) хорошо приближается частичной суммой некоторого порядка. Эта частичная сумма принимается в качестве приближенного аналитического решения задачи (1), (2) на промежутке [XN, XK]. В противном случае, т.е. при медленной сходимости ряда (3) на интервале [XN, XK], получение аналитического решения в виде одной частичной суммы на всем отрезке интегрирования [XN, XK] может быть затруднено. Поэтому целесообразно использовать разбиение промежутка интегрирования [XN, XK] на такие элементарные сегменты некоторой, вообще говоря, переменной длины H: [x_s, x_s+H], x₀ = XN, s = 0, 1, ... , на каждом из которых ряды Чебышёва для решения Y(X) и его производных Y'(X), Y''(X) будут сходиться значительно быстрее. На каждом подобном сегменте решение исходной задачи Коши приближенно представляется в виде (K + 2) - й частичной суммы смещенного ряда Чебышёва

                                                            _K+2                                                                        ₁
(6)     Y(x_s + αH) ≈ U_K+2(x_s + αH) = ∑^'a_i^*[U_K+2]T_i^*(α),   0 ≤ α ≤ 1,   a_i^*[U_K+2] = 2/π ∫ U_K+2(x_s+αH)T_i^*(α) / √(α(1-α))dα,
                                                            ⁱ⁼⁰                                                                         ⁰

а его производные - в виде частичных сумм (K + 1)-го и К-го порядков

                                                              _K+1                                                                         ₁
(7)     Y'(x_s + αH) ≈ U'_K+2(x_s + αH) = ∑^'a_i^*[U'_K+2]T_i^*(α),   0 ≤ α ≤ 1,   a_i^*[U'_K+2] = 2/π ∫ U'_K+2(x_s + αH)T_i^*(α) / √(α(1-α))dα,
                                                              ⁱ⁼⁰                                                                          ⁰

                                                                _K                                                                             ₁         
(8)     Y''(x_s + αH) ≈ U''_K+2(x_s + αH) = ∑^'a_i^*[U''_K+2]T_i^*(α),   0 ≤ α ≤ 1,   a_i^*[U''_K+2] = 2/π ∫ U''_K+2(x_s + αH)T_i^*(α) / √(α(1-α))dα.
                                                               ⁱ⁼⁰                                                                            ⁰

В этом случае аналитическое решение задачи (1), (2) состоит из совокупности частичных сумм рядов Чебышёва, построенных на этих элементарных сегментах. Порядок этих частичных сумм задается пользователем при обращении к подпрограмме.

При обращении к подпрограмме DE98D задаются начало элементарного сегмента и его длина. По заданным значениям решения и первой производной в начале X = x_s элементарного сегмента [x_s, x_s + H] подпрограмма DE98D вычисляет значения решения и первой производной в конце элементарного сегмента, т.е. в узле X + H = x_s + H. Одновременно вычисляются коэффициенты Чебышёва a_i^*[U_K+2] (i = 0, 1, ... , K + 2) на элементарном сегменте [X, X + H] = [x_s, x_s + H] для решения задачи Коши Y(X + αH) ≈ U_K+2(X + αH), 0 ≤ α ≤ 1, коэффициенты Чебышёва a_i^*[U'_K+2] (i = 0, 1, ... , K + 1) его первой производной U'_K+2(X + αH) и коэффициерты Чебышёва a_i^*[U''_K+2] (i = 0, 1, ... , K) второй производной U''_K+2(X + αH) = Φ(α). Значение H может изменяться от сегмента к сегменту, в общем случае H = H_s.

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

При разбиении промежутка интегрирования на элементарные сегменты решение задачи сводится к определению нескольких наборов коэффициентов a_i^*[U_K+2], i = 0, 1, ... , K + 2. Коэффициенты a_i^*[U_K+2] ряда Чебышёва для решения и коэффициенты ряда Чебышёва для U'_K+2(x_s + αH) на сегменте [x_s, x_s + H] выражаются через коэффициенты a_i^*[Φ] ряда Чебышёва второй производной Φ(α) = F(x_s + αH, U_K+2(x_s + αH), U'_K+2(x_s + αH)), 0 ≤ α ≤ 1, на [x_s, x_s + H], которые, в свою очередь, вычисляются приближенно итерационным способом, исходя из некоторого начального приближения. Вычисления выполняются с помощью квадратурной формулы Маркова на [x_s, x_s + H] с K + 1 узлом. При этом один из узлов квадратурной формулы совпадает с x_s, а остальные K узлов лежат внутри интервала (x_s, x_s + H). Количество итераций, которое предписывается выполнить в этом итерационном процессе, задается при обращении к подпрограмме и может меняться от сегмента к сегменту. Если при выбранном H ряды Чебышёва для U_K+2(x_s + αH), 0 ≤ α ≤ 1 , и его производных U'_K+2(X + αH), U''_K+2(X + αH) на элементарном сегменте [X, X + H] = [x_s, x_s + H] быстро сходятся, то для того, чтобы приближенное решение в конце одного такого сегмента имело максимальный порядок точности относительно H, необходимо выполнить не менее K итераций; при этом погрешность приближенного решения в конце элементарного сегмента является величиной порядка O(H^K + 3) при H --> 0, а погрешность приближенного значения производной Y' - величиной порядка O(H^K + 2) при H --> 0. Если H подобрано достаточно малым (или, вернее сказать, выбрано довольно удачным), то хорошая точность приближенного решения может быть получена и при меньшем числе итераций.

Начальное приближение коэффициентов a_i^*[Φ] ряда Чебышёва для второй производной на сегменте [x_s, x_s + H] может вычисляться двумя способами. В первом способе начальное приближение определяется только с использованием значения решения U_K+2(X) и его первой производной U'_K+2(X) в узле X = x_s. При этом погрешность начального приближения для всех коэффициентов a₀^*[Φ], a₁^*[Φ], ..., a_K^*[Φ] является величиной O(H²) при H --> 0. Во втором способе начальное приближение определяется через коэффициенты ряда Чебышёва производной Φ(α) на предыдущем элементарном сегменте [x_s - 1, x_s] = [X - H', X], где H' - длина предыдущего элементарного сегмента. В этом случае погрешности начального приближения для коэффициентов a₀^*[Φ], a₁^*[Φ], ... , a_K^*[Φ] имеют, соответственно, порядки O(H), O(H²), ..., O(H^K + 1). Заметим, что длина H текущего элементарного сегмента [X, X + H] может быть больше или меньше длины H' предыдущего элементарного сегмента [X - H', X] или равна ей. Второй способ определения начального приближения в некоторых случаях может привести к более быстрой сходимости итерационного процесса и, тем самым, к меньшему числу выполняемых итераций. Второй способ может быть применен только начиная со второго элементарного сегмента [x₀ + H₀, x₁ + H], H₀ - длина начального (первого) элементарного сегмента, x₁ = x₀ + H₀. На начальном элементарном сегменте [x₀, x₀ + H] всегда применяется исключительно первый способ. Способ выбора начального приближения задается пользователем при обращении к подпрограмме и может меняться от сегмента к сегменту.

На заданном при обращении к подпрограмме DE98D частичном сегменте [X, X + H]=[x_s, x_s + H] для оценки погрешности решения U_K+2(x_s + αH) и его производной U'_K+2(x_s + αH), полученных при некотором заданном значении K = K₁, вычисляются на этом же сегменте специальным образом второе приближенное решение вместе с его производной при K = K₂ > K₁, которые также представляются указанными выше частичными суммами при K = K₂. Второе решение U_K₂+2, как имеющее более высокий порядок точности O(H^K₂ + 3), используется для оценки погрешности первого решения U_K₁+2 на элементарном сегменте [X, X + H]. Производная второго решения U'_K₂+2, имеющая порядок точности O(H^K₂ + 2), используется для оценки погрешности производной первого решения U'_K₁+2. Значения K₁ и K₂ задаются при обращении к подпрограмме и могут изменяться от сегмента к сегменту. Второе приближенное решение на сегменте [X, X + H] вычисляется аналогично первому решению итерационным способом. Единственным отличием при этом является то, что начальным приближением в этом дополнителном итерационном процессе служит уже вычисленное на отрезке [X, X + H] первое решение U_K₁+2(x_s + αH). Количество итераций, которое потребуется выполнить в этом дополнительном итерационном процессе, задается при обращении к подпрограмме и может меняться от сегмента к сегменту.

Для системы уравнений (M > 1) с проверкой на точность могут вычисляться либо все компоненты решения, либо некоторые из них (в частности, одна компонента). В последнем случае номера проверяемых на точность компонент задаются при обращении к подпрограмме. Совершенно аналогично вычисляются с проверкой на точность компоненты производной решения.

Для оценки уклонения приближенного решения от точного (т.е. для оценки абсолютной погрешности приближенного решения) в подпрограмме DE98D используются два способа, или две формулы. Первый способ опирается на первую формулу для абсолютной погрешности приближенного решения. Она является асимптотической (т.е. справедливой при H --> 0) и представляет собою разность двух частичных сумм смещенного ряда Чебышёва на сегменте [X, X + H] с приближенными коэффициентами, а именно разность приближенных частичных сумм (K₂+2)-го и (K₁+2)-го порядков в точке X + H. Второй способ опирается на несколько завышенную оценку абсолютной погрешности, которая представляет сумму абсолютных величин (модулей) разностей приближенных коэффициентов этих частичных сумм и поэтому несколько завышена по сравнению с первой. Назовем эту формулу второй. Более того, при использовании этой, второй, формулы в том случае, когда точность решения оценивается по относительной погрешности или по мере погрешности, то дополнительно используется также завышенная оценка сверху для относительной погрешности и завышенная оценка сверху для меры погрешности. Второй способ оценки погрешности, хотя и является более точным, накладывает более сильное ограничение на размер элементарного сегмента и приводит к меньшей длине сегментов. Используемый способ оценки абсолютной погрешности задается при обращении к подпрограмме. Совершенно аналогичный прием применяется также для оценки уклонения производной приближенного решения от производной точного решения.

Если для заданной при обращении к подпрограмме DE98D длины H сегмента [X, X + H] точность приближенного решения U_K₁+2 или его производной U'_K₁+2 не достигается, то данный сегмент исключается из рассмотрения и программа сокращает его длину H до тех пор, пока на новом, сокращенном, сегменте не будет достигнута заданная точность вновь вычисленного приближенного решения U_K₁+2 и его производной. При выполнении условия достижения на сегменте заданной точности приближенного решения U_K₁+2 и его производной этот сегмент (заданной, первоначальной, длины или сокращенной длины) принимается в качестве элементарного (частичного) сегмента, на который продлено решение дифференциального уравнения. В качестве такого решения на сегменте [X, X + H] принимается частичная сумма порядка K₁ + 2 с коэффициентами, равными коэффициентам Чебышёва a_i^*[U_K₂+2], i = 0, 1, ... , K₁ + 2, второго приближенного решения U_K₂+2, поскольку эти коэффициенты имеют более высокий порядок точности по сравнению с коэффициентами Чебышёва a_i^*[U_K₁+2] первого приближенного решения. В качестве значения решения в конце X + H сегмента [X, X + H] принимается значение второго приближенного решения в точке X + H, т.е. U_K₂+2(X + H), как имеющее более высокий порядок точности O(H^K₂+3) по сравнению с первым приближенным решением U_K₁+2(X + H). В качестве значения производной решения в конце X + H сегмента [X, X + H] принимается значение производной второго приближенного решения в точке X + H, т.е. U'_K₂+2(X + H), поскольку оно имеет более высокий порядок точности O(H^K₂+2) по сравнению с производной первого приближенного решения U'_K₁+2(X + H). В качестве коэффициентов Чебышёва производной решения на сегменте [X, X + H] принимаются по той же причине коэффициенты Чебышёва производной второго приближенного решения a_i^*[U'_K₂+2], i = 0, 1, ... , K₁ + 1, а в качестве коэффициентов Чебышёва второй производной решения принимаются коэффициенты Чебышёва второй производной второго приближенного решения a_i^*[U''_K₂+2], i = 0, 1, ... , K₁.

Таким образом, один шаг приближенного интегрирования по методу рядов Чебышёва, выполняемый подпрограммой DE98D, заключается в следующем. По заданным значениям решения YX и его производной DYX в начале X сегмента [X, X + H] вычисляются с контролем точности значения этого решения и его производной в конце X + H сегмента, а также определяются значения коэффициентов Чебышёва решения, значения коэффициентов Чебышёва его производной и значения коэффициентов Чебышёва его второй производной на сегменте [X, X + H]; другими словами, находятся ортогональное разложение решения и ортогональные разложения его первой и второй производных на сегменте [X, X + H].

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

Залеткин С.Ф. Численное интегрирование обыкновенных дифференциальных уравнений с использованием ортогональных разложений // Математическое моделирование. 2010. 22. № 1. 69 - 85.

Арушанян О.Б., Волченскова Н.И., Залеткин С.Ф. О применении ортогональных разложений для приближенного интегрирования обыкновенных дифференциальных уравнений // Вестник Московского университета. Серия 1. Математика. Механика. 2010. № 4. 40 - 43.

Арушанян О.Б., Волченскова Н.И., Залеткин С.Ф. Вычисление коэффициентов разложения решения задачи Коши в ряд по многочленам Чебышёва // Вестник Московского университета. Серия 1. Математика. Механика. 2012. № 5. 24 - 30.

Арушанян О.Б., Залеткин С.Ф. Обоснование одного подхода к применению ортогональных разложений для приближенного интегрирования канонических систем обыкновенных дифференциальных уравнений второго порядка // Вестник Московского университета. Серия 1. Математика. Механика. 2018. № 3. 29 - 33.

Арушанян О.Б., Залеткин С.Ф. К теории вычисления ортогонального разложения решения задачи Коши для обыкновенных дифференциальных уравнений второго порядка // Вычислительные методы и программирование. 2018. 19. 178 - 184.

Арушанян О.Б., Залеткин С.Ф. Приближенное интегрирование канонических систем обыкновенных дифференциальных уравнений второго порядка методом рядов Чебышёва с оценкой погрешности решения и его производной // Вестник Московского университета. Серия 1. Математика. Механика. 2022. № 4. 27 - 34.

Залеткин С.Ф. Приближенное интегрирование обыкновенных дифференциальных уравнений методом рядов Чебышёва с контролем точности // Математическое моделирование. 2022. 34. № 6. 53 - 74.

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

    SUBROUTINE  DE98D(F, M, K, INIAPR, IMAX, JSTART, YX, DYX, X, IU, EPS, THRESH, H, HMIN, HMAX, BUL, XP, YP, DYP, MEXACT, NUMBES,
                      DELTY, RO, MESTER, K2, IMAX2, YX2, DYX2, NATTEM, AU, ADU, AJK, AU2, ADU2, AJK2, XJ0, XJ02, AJKP, U, DU,
                      U2, DU2, ZFI, ZFI2, HD4I, RABC, RABC1, RABC2, RABC3, ALP, ALPN, CMAR, P, S, ALP2, ALPN2, CMAR2, IERR)

   Параметры

F -	имя подпрограммы вычисления значений правой части дифференциального уравнения. Первый оператор подпрограммы должен иметь вид: SUBROUTINE F (X, Y, DY, D2Y, M). Здесь: X, Y, DY - значения независимой, зависимой переменных и производной решения соответственно. Вычисленное значение правой части должно быть помещено в D2Y. B случае системы уравнений, т.е. когда M ≠ 1, параметры Y, DY, D2Y представляют одномерные массивы длины M (тип параметров X, Y, DY и D2Y: с двойной точностью);
M -	количество уравнений в системе (тип: целый);
K -	порядок частичной суммы смещенного ряда Чебышёва, с помощью которой аппроксимируется вторая производная первого приближенного решения задачи Коши на каждом элементарном сегменте разбиения интервала интегрирования; при этом само первое приближенное решение задачи Коши приближается на каждом элементарном сегменте частичной суммой (K + 2) - го порядка, а его производная - частичной суммой (K + 1) - го порядка; K≥2. Этот параметр относится только к первому приближенному решению, имеющему меньший порядок точности O(H^K+3), и его производной (см. "Математическое описание" и "Замечания по использованию"; тип: целый);
INIAPR -	целый указатель способа выбора начального приближения коэффициентов Чебышёва для второй производной первого приближенного решения на каждом элементарном сегменте:

INIAPR=1 -

для первого способа, когда начальное приближение определяется только с использованием значения решения и его первой производной в начале каждого элементарного сегмента;

INIAPR=2 -

для второго способа, когда начальное приближение коэффициентов Чебышёва на текущем элементарном сегменте (начиная со второго) определяется через коэффициенты Чебышёва, вычисленные на предыдущем элементарном сегменте, т.е. путем экстраполяции коэффициентов с предыдущего сегмента на следующий (см. "Математическое описание");

IMAX -

целая переменная, задающая количество итераций, которое предполагается выполнить в итерационном процессе вычисления коэффициентов Чебышёва для второй производной первого приближенного решения задачи Коши на элементарном сегменте [X, X + H], исходя из некоторого начального приближения, способ определения которого задается параметром INIAPR; IMAX ≥ 1. Для получения максимального порядка точности первого приближенного решения необходимо выполнить не менее K итераций. Этот параметр относится только к первому приближенному решению, имеющему меньший порядок точности O(H^K + 3), и его производной. Если правая часть дифференциального уравнения не зависит от переменных Y и Y', т.е. дифференциальное уравнение имеет вид Y'' = F(X), то число итераций при вычислении первого приближенного решения можно положить равным 1. В этом случае параметр IMAX = 1 (см. "Математическое описание" и "Замечания по использованию");

JSTART -

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

0 -

первое обращение к подпрограмме должно быть исполнено с нулевым значением JSTART: выполнить первый (начальный) шаг интегрирования для значений независимой переменной X, зависимой переменной YX, производной DYX и длины элементарного сегмента H. При данном значении параметра JSTART будет применен исключительно первый способ определения начального приближения для коэффициентов Чебышёва второй производной первого приближенного решения независимо от значения параметра INIAPR. Нулевое значение параметра JSTART также может означать, что необходимо выполнить очередной шаг интегрирования с измененным значением параметра K или с измененным значением параметра K2;

1 -

выполнить следующий (очередной) шаг интегрирования системы дифференциальных уравнений для значений независимой перемннной X, зависимой переменной YX, производной DYX и длины элементарного сегмента H. При данном значении параметра JSTART способ определения начального приближения для коэффициентов Чебышёва второй производной решения определяется параметром INIAPR. При обращении к подпрограмме DE98D со значением JSTART = 1 значения параметра K (см. выше) и параметра K2 (см. ниже) должны совпадать с их значениями при предыдущем обращении к подпрограмме.

	На выходе из подпрограммы DE98D параметр JSTART всегда принимает значение, равное 1 (см. "Примеры использования");
X, YX, DYX -	начальные значения аргумента, решения и первой производной (начало элементарного сегмента X, решение YX и производная DYX в нем); в результате работы подпрограммы в X получается новое значение аргумента, отстоящее от начального значения на величину выбранного подпрограммой элементарного сегмента, т.е. такого сегмента, на который продлено решение дифференциального уравнения (конец элементарного сегмента), а в YX и DYX - соответствующие значения решения и производной; в качестве такого решения принимается второе приближенное, оценивающее, решение, поскольку оно имеет более высокий порядок точности по сравнению с первым приближенным решением. По этой же причине в качестве производной принимается производная второго приближенного (оценивающего) решения. В случае системы уравнений, т.е. когда M ≠ 1, YX и DYX задаются одномерными массивами длины M (тип параметров X, YX, DYX: с двойной точностью);
IU -	целый указатель типа погрешности численного решения и его производной; задается одномерным массивом длины 2. Переменная с индексом IU(1) задает тип погрешности производной решения, а переменная с индексом IU(2) задает тип погрешности решения:

IU(1)=1 -	компоненты производной решения проверяются на точность по абсолютной погрешности;
IU(1)=2 -	компоненты производной решения проверяются на точность по относительной погрешности;
IU(1)=3 -	компоненты производной решения проверяются на точность по мере погрешности;
IU(2)=1 -	компоненты решения проверяются на точность по абсолютной погрешности;
IU(2)=2 -	компоненты решения проверяются на точность по относительной погрешности;
IU(2)=3 -	компоненты решения проверяются на точность по мере погрешности.

	Контроль точности по мере погрешности заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше определенной наперед заданной положительной константы THRESH (называемой границей перехода), то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Таким образом, контроль точности по мере погрешности состоит в том, что на тех участках интервала интегрирования, где абсолютная величина компоненты решения меньше значения THRESH, контроль точности ведется для нее по абсолютной погрешности, а там, где абсолютная величина этой компоненты решения равна значению THRESH или превосходит значение THRESH, контроль точности для нее ведется по относительной погрешности; аналогично выполняется контроль точности по мере погрешности и для производной решения (см. "Замечания по использованию");
EPS -	допустимая погрешность, с которой требуется вычислить проверяемые на точность компоненты решения и его производной; тип погрешности специфицируется с помощью параметра IU; задается одномерным массивом длины 2. Переменная с индексом EPS(1) задает точность для производной решения, а переменная с индексом EPS(2) задает точность для решения (см. "Замечания по использованию"; тип: с двойной точностью);
THRESH -	граница перехода, используемая при оценке меры погрешности решения и его производной; задается одномерным массивом длины 2. Переменная с индексом THRESH(1) задает границу перехода для производной решения, а переменная THRESH(2) задает границу перехода для решения (см. "Замечания по использованию"; тип: с двойной точностью);
H -	переменная с двойной точностью, содержащая начальное значение длины первого элементарного сегмента (аналог шага интегрирования для разностных методов). Может задаваться с учетом направления интегрирования, т.е. положительным, если XK > XN, отрицательным, если XK < XN , или без такого учета в виде абсолютной величины;
HMIN, HMAX -	соответственно минимальное и максимальное значение длины элементарного сегмента, которое разрешается использовать при разбиении интервала интегрирования на элементарные (частичные) сегменты (тип: с двойной точностью);
BUL -	логическая переменная, значение которой при обращении к подпрограмме полагается равным .TRUE., если заданная в H длина сегмента выводит в конец интервала интегрирования, и .FALSE. в противном случае; в результате работы подпрограммы BUL равно .FALSE., если вместо исходного сегмента был выбран меньший сегмент, в противном случае, т.е. когда был выбран именно заданный при обращении в H сегмент, значение параметра BUL не меняется.
XP, YP, DYP -	переменная XP и одномерные массивы YP, DYP длины M. На выходе из подпрограммы содержат начальные значения аргумента, решения и производной, т.е. те значения, которые имели параметры X, YX, DYX на входе в подпрограмму. Таким образом, на выходе из подпрограммы DE98D границы элементарного сегмента, к которому относятся вычисленные подпрограммой коэффициенты Чебышёва решения и коэффициенты Чебышёва его производных и содержащиеся в массивах AU, ADU и AJK, показываются параметрами XP и X, а именно: параметр XP содержит начало сегмента, а параметр X содержит конец этого сегмента (тип: с двойной точностью);
MEXACT -	количество компонент численного решения и его производной, проверяемых на точность; задается одномерным массивом длины 2. Переменная с индексом MEXACT(1) показывает, сколько компонент производной надо проверять на точность; переменная с индексом MEXACT(2) показывает, сколько компонент решения надо проверять на точность. Нулевое значение переменной с индексом MEXACT(1) или MEXACT(2) сообщает о том, что соответственно производная или решение на точность не проверяется. Например, элементы массива MEXACT/O, M/ сообщают о том, что производная на точность вообще не проверяется, а проверяются на точность все компоненты решения (здесь M - количество уравнений в системе). Элементы массива MEXACT/M, O/ указывают, что на точность надо проверять все компоненты производной, а решение на точность проверять не надо. Элементы массива MEXACT/M, M/ сообщают, что все компоненты производной и все компоненты решения надо проверять на точность. Элементы массива MEXACT/1, 1/ уведомляют, что на точность проверяется какая-то одна компонента производной и какая-то одна компонента решения (тип: целый);
NUMBES -	одномерный массив, содержащий номера проверяемых на точность компонент производной и номера проверяемых на точность компонент решения. Первые MEXACT(1) элементов массива NUMBES содержат номера проверяемых на точность компонент производной, а следующие MEXACT(2) элементов этого массива NUMBES содержат номера проверяемых на точность компонент решения. Следовательно, NUMBES представляет одномерный массив длины MEXACT(1) + MEXACT(2). Если MEXACT(1) = 0, то массив NUMBES содержит только номера проверяемых на точность компонент решения; если MEXACT(2) = 0, то массив NUMBES содержит только номера проверяемых на точность компонент производной. В том случае, если все компоненты производной и все компоненты решения проверяются на точность (этому случаю соответствуют значения MEXACT(1) = MEXACT(2) = M), этот массив в подпрограмме не используется. Таким образом, массив NUMBES используется в подпрограмме только в том случае, когда число проверяемых на точность компонент производной или число проверяемых на точность компонент решения меньше числа уравнений M (тип: целый);
DELTY, RO -	одномерные рабочие массивы длины M (тип: с двойной точностью);
MESTER -	целый указатель способа оценки абсолютной погрешности компоненты приближенного решения и его производной:

MESTER=1 -	для оценки абсолютной погрешности используется первая формула (асимптотическая);
MESTER=2 -	для оценки абсолютной погрешности используется вторая формула (завышенная оценка).

	(См. "Математическое описание");
K2 -	порядок частичной суммы смещенного ряда Чебышёва второй производной решения, с помощью которой вычисляется частичная сумма порядка (K2 + 2) для оценивающего решения более высокого порядка точности на элементарном сегменте [X, X + H] разбиения области интегрирования. Этот параметр относится только ко второму приближенному решению более высокого порядка точности O(H^K2+3), K2 > K (т.е. к оценивающему решению) и его второй производной. Если обращение к подпрограмме DE98D осуществляется со значением параметра JSTART = 1 (см. выше), то значение параметра K2 должно совпадать с его значением при предыдущем обращении к подпрограмме. Если же при очередном обращении к подпрограмме DE98D значение параметра K2 необходимо изменить, то это обращение должно выполняться при нулевом значении параметра JSTART (JSTART = 0) (см. "Математическое описание", "Замечания по использованию" и "Примеры"; тип: целый);
IMAX2 -	целая переменная, задающая количество итераций, которое предполагается выполнить в дополнительном итерационном процессе для вычисления второго, оценивающего, приближенного решения более высокого порядка точности. Для получения максимального порядка точности оценивающего решения необходимо выполнить не менее (K2 - K) итераций. Значение IMAX2 может изменяться от сегмента к сегменту. Если правая часть дифференциального уравнения не зависит от переменных Y и DY, т.е. дифференциальное уравнение имеет вид Y'' = F(X), то число итераций при вычислении второго приближенного решения можно положить равным 1. В этом случае параметр IMAX2 = 1 (см. "Математическое описание" и "Замечания по использованию");
YX2, DYX2 -	одномерные рабочие массивы длины M (тип: с двойной точностью);
NATTEM -	целая переменная, значение которой ограничивает число последовательных сокращений длины элементарного сегмента [X, X + H], если приближенное решение или его производная на этом элементарном сегменте не достигает требуемой точности;
AU -	двумерный массив размера M * (K + 3). На выходе из подпрограммы содержит коэффициенты Чебышёва a_i*[Y] для решения Y (X + αH), 0 ≤ α ≤ 1, на элементарном сегменте [X, X + H]. При этом переменная с индексом AU (N, I + 1) представляет I-й коэффициент Чебышёва N-й компоненты решения y_N(x)(I = 0, 1, ... , K + 2) (тип: с удвоенной точностью);
ADU -	двумерный массив размера M * (K + 2). На выходе из подпрограммы содержит коэффициенты Чебышёва a_i*[Y'] для первой производной решения Y' (X + αH), 0 ≤ α ≤ 1, на элементарном сегменте [X, X + H]. При этом переменная с индексом ADU (N, I + 1) представляет I-й коэффициент Чебышёва N-й компоненты первой производной решения y'_N(x) (I = 0, 1, ... , K + 1) (тип: с удвоенной точностью);
AJK -	двумерный массив размера M * (K + 1). На выходе из подпрограммы содержит коэффициенты Чебышёва a_i*[Y''] второй производной решения Y''(X + αH), 0 ≤ α ≤ 1, на элементарном сегменте [X, X + H]. При этом переменная с индексом AJK(N, I + 1) представляет I-й коэффициент Чебышёва N-й компоненты второй производной решения y''_N(x) (I = 0, 1, ... , K). Если при JSTART = 1 (т.е. при повторных обращениях к подпрограмме) параметр INIAPR = 2, то содержащиеся в массиве AJK при входе в подпрограмму коэффициенты Чебышёва второй производной решения, вычисленные на предыдущем элементарном сегменте [X - H', X] во время предыдущего обращения к подпрограмме DE98D, используются в подпрограмме при вычислении коэффициентов Чебышёва второй производной решения на текущем элементарном сегменте [X, X + H]. Заметим, что длина H текущего элементарного сегмента может быть больше или меньше длины H' предыдущего сегмента [X - H', X] или равна ей. Независимо от значения параметра INIAPR при повторных обращениях к подпрограмме значения массива AJK при входе в подпрограмму обязательно запоминаются в массиве AJKP и могут быть снова доступны на выходе из подпрограммы (тип: с удвоенной точностью);
AU2 -	двумерный рабочий массив размера M * (K2 + 3) (тип: с двойной точностью);
ADU2 -	двумерный рабочий массив размера M * (K2 + 2) (тип: с двойной точностью);
AJK2 -	двумерный рабочий массив размера M * (K2 + 1) (тип: с двойной точностью);
XJ0 -	одномерный рабочий массив длины K (тип: с удвоенной точностью);
XJ02 -	одномерный рабочий массив длины K2 (тип: с удвоенной точностью);
AJKP -	двумерный массив размера M * (K + 1). Если обращение к подпрограмме было выполнено со значением JSTART = 1 (т.е. при повторных обращениях к подпрограмме), то на выходе из подпрограммы содержит коэффициенты Чебышёва второй производной решения, относящиеся к предыдущему элементарному сегменту [X - H', X]; тем самым значение этого параметра определено только при повторных обращениях к подпрограмме (тип: с двойной точностью);
U, DU -	двумерные рабочие массивы размера M * K (тип: с двойной точностью);
U2, DU2 -	двумерные рабочие массивы размера M * K2 (тип: с двойной точностью);
ZFI -	двумерный рабочий массив размера M *(K + 4) (тип: с двойной точностью);
ZFI2 -	двумерный рабочий массив размера M *(K2 + 3) (тип: с двойной точностью);
HD4I -	одномерный рабочий массив длины K2 + 1 (тип: с двойной точностью);
RABC, RABC1 -	одномерные рабочие массивы длины K2 (тип: с двойной точностью);
RABC2, RABC3 -	одномерные рабочие массивы длины K2 + 2 (тип: с двойной точностью);
ALP, ALPN -	одномерные рабочие массивы длины K (тип: с двойной точностью);
CMAR -	двумерный рабочий массив размера K * K (тип: с двойной точностью);
P -	одномерный рабочий массив размера (K + 2) * (K - 1) / 2 (тип: с двойной точностью);
S -	одномерный рабочий массив размера (K + 3) * K / 2 (тип: с двойной точностью);
ALP2, ALPN2 -	одномерные рабочие массивы длины K2 (тип: с двойной точностью);
CMAR2 -	двумерный рабочий массив размера K2 * K2 (тип: с двойной точностью);
IERR -	целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом:

IERR=0 -

нулевое значение параметра IERR свидетельствует об успешном окончании работы подпрограммы DE98D;

IERR=65 -
IERR=66

когда какая-нибудь компонента решения или производной не может быть вычислена с требуемой точностью EPS; при этом IERR=65 указывает, что требуемая точность не может быть достигнута на сегменте [X, X + H] при минимальной длине частичного сегмента H, равной HMIN, а IERR=66 показывает, что требуемая точность не может быть достигнута на сегменте [X, X + H], т.к. исчерпано заданное число сокращений NATTEM длины H элементарного сегмента. В обоих случаях интегрирование системы уравнений прекращается.

   Версии:  нет

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

DE98D использует рабочие подпрограммы DE70DH, DE80DH, DE80DI, DE70DF, DE70DQ, DE71DE, DE70DP, DE71DT, DE71DP, DE71DI, DE71DF, DE71DS, DE70DA, DE70DC, DE75DK, DE75DE, DE75DN, DE98DK, DE98D1, DE98D0.

UTDE75 - подпрограмма выдачи диагностических сообщений.

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

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

Разбиение промежутка интегрирования на элементарные сегменты выполняется для того, чтобы на каждом таком сегменте ряды Чебышёва для решения и его первой и второй производных были быстросходящимися рядами; другими словами, чтобы убывание коэффициентов этих рядов Чебышёва на элементарном сегменте происходило достаточно быстро, вследствие чего можно было бы считать частичные суммы этих рядов близкими к многочленам наилучшего равномерного приближения на элементарном сегменте для решения и его производных.

Если при вычислении первого приближенного решения начальное приближение для коэффициентов Чебышёва правой части системы (т.е. функции Φ(α)) определяется первым способом (т.е. при INIAPR = 1), то для получения максимального порядка точности первого приближенного решения в конце элементарного сегмента необходимо выполнить в итерационном процессе не менее K итераций; тогда IMAX ≥ K. Если начальное приближение коэффициентов Чебышёва функции Φ(α) определяется вторым способом (т.е. при INIAPR = 2), то для получения максимального порядка точности первого приближенного решения необходимо выполнить в итерационном процессе не менее K + 1 итераций; в этом случае IMAX ≥ K + 1. Однако в некоторых случаях при втором способе определения начального приближения итерационный процесс может сойтись за значительно меньшее число итераций.

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

Для получения максимального порядка точности второго приближенного (оценивающего) решения необходимо выполнить в дополнительном итерационном процессе не менее (K2 - K) итераций. Указанное в данном разделе число итераций IMAX и IMAX2 носит асимптотический характер (когда все наши оценки справедливы при H --> 0). На практике же количество итераций в каждом из двух итерационных процессов зависит от длины сегмента, от заданной точности, от используемых порядков частичных сумм и от конкретной системы дифференциальных уравнений. Поэтому числа итераций в обоих итерационных процессах, т.е. значения параметров IMAX и IMAX2, могут быть как меньше, так и больше рекомендуемых здесь значений.

При работе подпрограммы значения параметров M, K, INIAPR, IMAX, IU, EPS, THRESH, MEXACT, NUMBES, MESTER, HMIN, HMAX, K2, IMAX2, NATTEM сохраняются. При многократном использовании подпрограммы DE98D со значением параметра JSTART = 1 для вычисления коэффициентов Чебышёва решения задачи Коши (1), (2) и его производных на последовательности элементарных сегментов, образующей промежуток интегрирования [XN, XK], значения параметров M, K, YX, DYX, X, AU, ADU, AJK, HD4I, RABC, RABC1, RABC2, RABC3, ALP, ALPN, CMAR, P, S, K2, ALP2, ALPN2, CMAR2 не должны изменяться в вызывающей программе между последовательными обращениями к подпрограмме.

Подпрограмма DE98D использует общие блоки с именами COM70D, COM80D, COM75D, COM98D. Поэтому, если вызывающая (главная) программа содержит хотя бы одно обращение к DE98D со значением параметра JSTART = 1 (и, следовательно, число обращений к DE98D более одного), то данные общие блоки должны быть описаны в вызывающей (главной) программе, например, таким способом:

COMMON / COM70D / WC1, WC2, WC3, LASN, HD2, HD4, HOLD

COMMON / COM80D / JLAS, JLAS1, JLAS2, JLAS3, H2D4, H2D8, H2D16, H2D32, H2D96

COMMON / COM75D / LASN0, LASN2, WC30, WC12, WC22, WC32

COMMON / COM98D / JLASA2, JLAS12, JLAS22, JLAS32, JLAS0, JLAS10, JLAS20, JLAS30

Глобальные (общие) переменные, составляющие содержание общих блоков, описываются как переменные с двойной точностью:

DOUBLE PRECISION WC1, WC2, WC3, HD2, HD4, HOLD, H2D4, H2D8, H2D16, H2D32, H2D96, WC30, WC12, WC22, WC32

Остальные общие переменные описываются как целые переменные (явно или неявно). Пользователю не рекомендуется применять для своих целей общие блоки с указанными именами. Пользователь не должен изменять содержимое указанных переменных из данных блоков. Если вызывающая программа содержит одно обращение к подпрограмме DE98D со значением параметра JSTART = 0 и при этом других обращений к DE98D нет или если вызывающая программа содержит несколько обращений к DE98D и при этом все обращения осуществляются с нулевым значением параметра JSTART, то указанные общие блоки описывать в вызывающей программе вовсе не обязательно (см. "Примеры использования").

При работе подпрограммы DE98D используется общий блок

COMMON/STAT76/NACCEP, NREJEC

для сбора полезной статистической информации. Такой же блок должен содержаться в вызывающей (главной) программе. Если подпрограмма DE98D заканчивает свою работу со значением IERR = 0, то на выходе из подпрограммы переменная NACCEP содержит значение, увеличенное на 1 по сравнению со значением этой переменной на входе в подпрограмму. Если же подпрограмма DE98D заканчивает свою работу со значением IERR = 65 или IERR = 66, то значение переменной NACCEP остается без изменения.

Переменная NREJEC учитывает, сколько раз сокращалась длина сегмента [X, X + H] при вычислении решения и его производной с требуемой точностью. Если заданная точность была достигнута при первоначальной длине сегмента (задаваемой параметром H при входе в подпрограмму DE98D), то значение переменной NREJEC на выходе из подпрограммы равно значению этой переменной на входе в подпрограмму; если же первоначально заданная длина H подвергалась сокращению, то значение переменной NREJEC (которое эта переменная имела на входе в подпрограмму) увеличивается на столько, сколько раз сокращалась длина сегмента.

Таким образом, переменная NACCEP подсчитывает число выбранных (или принятых) подпрограммой DE98D сегментов, а переменная NREJEC подсчитывает число отклоненных сегментов при решении дифференциальных уравнений на всем промежутке интегрирования.

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

Данные примеры не только иллюстрируют правила использования подпрогрвммы DE98D; они также показывают, как можно применять подпрограмму DE98D для решения обыкновенных дифференциальных уравнений второго порядка на заданном промежутке интегрирования. Вычисления во всех примерах проводились на Фортране с 15-16 значащими цифрами.

1) В первом примере решается задача Коши

       y'' = qy',   y(0) = e^q,   y'(0) = qe^q,   q = 4,   0 ≤ x ≤ x_f,   x_f = 7.

Точное решение задачи Коши y(x) = e^q(1+x). Используется шесть сегментов разбиения промежутка интегрирования [0, x_f] и на всех шести сегментах решение приближается частичными суммами ряда Чебышёва одного и того же порядка K + 2 при K = 18; погрешность этого решения вычисляется с помощью оценивающего решения, которое представляется частичной суммой тоже одного и того же порядка K2 + 2 при K2 = 25. Производная приближается соответственно частичной суммой ряда Чебышёва порядка K + 1, а ее погрешность вычисляется с помощью производной оценивающего решения, которая представляется частичной суммой порядка K2 + 1. Пример содержит шесть обращений к подпрограмме DE98D. Приводятся вызывающая программа, подпрограмма вычисления правой части дифференциального уравнения, а также результаты счета на произвольном элементарном сегменте после каждого обращения к подпрогрмме. Результаты включают границы элементарного сегмента [XP, X]; приближенное значение YX решения в конце X сегмента, точное значение YT решения в X, относительную погрешность (YT - YX) / YX приближенного решения; приближенное значение DYX производной в конце сегмента X, точное значение производной DYT в X, относительную погрешность (DYT - DYX) / DYX производной; рекомендуемую длину следующего элементарного сегмента H, значения параметров K, IMAX, K2, IMAX2, значения параметра JSTART на входе в подпрограмму и на выходе из нее, значения переменных NACCEP, NREJEC из общего блока и число обращений к правой части N; значение логической переменной BUL. Приводятся также вычисленные на каждом элементарном сегменте [XP, X] коэффициенты Чебышёва AU для решения Y, коэффициенты Чебышёва ADU для его первой производной Y' и коэффициенты Чебышёва AJK для второй производной Y''.

 DOUBLE PRECISION X,YX,DYX,EPS,THRESH,H,HMIN,HMAX,XP,YP,DYP,DELTY,
     1                 RO,YX2,DYX2,AU,ADU,AJK,AU2,ADU2,AJK2,XJ0,XJ02,
     2                 AJKP,U,DU,U2,DU2,ZFI,ZFI2,HD4I,RABC,RABC1,RABC2,
     3                 RABC3,ALP,ALPN,CMAR,P,S,ALP2,ALPN2,CMAR2
C
C
      LOGICAL BUL
C
      PARAMETER (K=18,K2=25,M=1)
      DIMENSION YP(M),DYP(M),DELTY(M),RO(M),YX2(M),DYX2(M),
     1          IU(2),EPS(2),THRESH(2),MEXACT(2),NUMBES(2),
     2          XJ0 (K),  AU (M,K+3),  ADU (M,K+2), AJK (M,K+1), 
     3          XJ02(K2), AU2(M,K2+3), ADU2(M,K2+2),AJK2(M,K2+1),
     4          AJKP(M,K+1), U(M,K), DU(M,K), U2(M,K2),DU2(M,K2),
     5          HD4I(K2+1),RABC(K2),RABC1(K2),RABC2(K2+2),RABC3(K2+2),
     6          ALP (K), ALPN(K),  CMAR (K,K),  P((K+2)*(K-1)/2),
     7          ALP2(K2),ALPN2(K2),CMAR2(K2,K2),S((K+3)*K/2),
     8          ZFI(M,K+4), ZFI2(M,K2+3)
C
      DOUBLE PRECISION WC1,WC2,WC3,HD2,HD4,HOLD,WC30,WC12,WC22,WC32,
     1                 H2D4,H2D8,H2D16,H2D32,H2D96
C
      COMMON/COM70D/WC1,WC2,WC3,LASN,HD2,HD4,HOLD
      COMMON/COM80D/JLAS,JLAS1,JLAS2,JLAS3,H2D4,H2D8,H2D16,H2D32,H2D96
      COMMON/COM75D/LASN0,LASN2,WC30,WC12,WC22,WC32
      COMMON/COM98D/JLASA2,JLAS12,JLAS22,JLAS32,
     1              JLAS0, JLAS10,JLAS20,JLAS30                         25.03.21   
C
      COMMON/STAT76/NACCEP,NREJEC
      COMMON/BLOCK/N,Q
      DOUBLE PRECISION Q,YT,DYT,XK      
C
      EXTERNAL F
C
C
      DATA IU/2,2/, EPS/2*0.5D-11/, THRESH/2*1.D0/, MEXACT/1,1/,
     1     NUMBES/1,1/ 
C
      NREJEC=0
      NACCEP=0
C
      M0=M
      Q=4.0D0
      N=0
      X=0.
      YX=EXP(Q)
      DYX=Q*YX
      XK=7.D0
      K0=K
      INIAPR=1
      IMAX=28
      MESTER=1
      H=1.D0
      HMIN=1.D-3
      HMAX=XK
      K20=K2
      IMAX2=3
      NATTEM=3
C
      JSTART=0
      BUL=.FALSE.
C
      KP1=K0+1
      KP2=K0+2
      KP3=K0+3
C
      DO 20 NX=1,6
      PRINT 54,JSTART
      CALL DE98D(F,M0,K0,INIAPR,IMAX,JSTART,YX,DYX,X,IU,EPS,THRESH,     
     1                 H,HMIN,HMAX,BUL,XP,YP,DYP,MEXACT,NUMBES,DELTY,RO,
     2                 MESTER, K20,IMAX2,YX2,DYX2,NATTEM,AU,ADU,AJK,AU2,
     3                 ADU2,AJK2,XJ0,XJ02,AJKP,U,DU,U2,DU2,ZFI,ZFI2,
     4                 HD4I,RABC,RABC1,RABC2,RABC3,
     5                 ALP,ALPN,CMAR,P,S,ALP2,ALPN2,CMAR2,IERR)         
C
C
C
      YT=EXP(Q*(X+1.D0))
      DYT=Q*YT
      PRINT 44,NX
      PRINT 50,IERR,N,NACCEP,NREJEC,BUL
      PRINT 51,XP,X
      PRINT 35, H,K0,INIAPR,IMAX,JSTART,K20,IMAX2
      PRINT 36, H
      PRINT 52, YX,DYX,YT,DYT,(YT-YX)/YX,(DYT-DYX)/DYX
C
      PRINT 110,NX
      DO 10 L=1,M
      PRINT 47
      PRINT 150,L
      PRINT 85,(J-1,AU(L,J  ),ADU(L,J  ),AJK(L,J),J=1,KP1)
      PRINT 86, KP1,AU(L,KP2),ADU(L,KP2),
     1          KP2,AU(L,KP3)
   10 CONTINUE
C
      PRINT 55
C
      IF (NX .NE. 5) GO TO 20
C
      H=XK-X
      BUL=.TRUE.
   20 CONTINUE
C
      STOP
C
   35 FORMAT(1X,'H =',D25.16,3X,'K =',I2,3X,'INIAPR=',I1,3X,'IMAX =',I2,
     1       3X,'JSTART (exit) =',I1/1X,31X,'K2=',I2,13X,' IMAX2=',I2)
   36 FORMAT(1X,'It is advisable to take the next segment H = ',D22.16/)
   44 FORMAT(/1X,'RESULTS AFTER ',I1,'-nd    CALL SUBROUTINE DE98D'//)
   47 FORMAT(1X,64(1H-)/)
   50 FORMAT(1X,'IERR=',I2,3X,'N=',I5,3X,'NACCEP=',I3,3X,'NREJEC=',I3,
     1           3X,'BUL=',L3)
C
   51 FORMAT(1X,'XP=',D25.16,5X,'X=',D25.16)
C
   52 FORMAT( 1X,'YX =',D25.16,8X,'DYX =',D25.16
     1       /1X,'YT =',D25.16,8X,'DYT =',D25.16
     2       /1X,'(YT-YX)/YX =',D23.16,4X,'(DYT-DYX)/DYX =',D23.16/)
C
   54 FORMAT(1X,'JSTART (input) =',I1)
   55 FORMAT(1X,64(1H*)/)
   85 FORMAT(5X,I2,4X,3D23.16)
   86 FORMAT(5X,I2,4X,2D23.16)
  110 FORMAT(/12X,'Coefficients AU, ADU and AJK  on ',I3,'  segment  ')
C      
  150 FORMAT('   Number of',15X,'Chebyshev coefficients'
     *       ' for ',I1,' component'/  
     1        14X,72(1H-)/
     2       ' coefficient',9X,'for Y',18X,'for Y''',18X,'for Y'''''/         
     3        1X,85(1H-)/)
C   
      END
C         
C
      SUBROUTINE F(X,Y,DY,Z,M)
      DOUBLE PRECISION X,Y,DY,Z,Q
      COMMON/BLOCK/N,Q
      Z=Q*DY
      N=N+1
      RETURN
      END
              
Результаты:


JSTART (input) =0

RESULTS AFTER 1-nd    CALL SUBROUTINE DE98D


IERR= 0   N=  666   NACCEP=  1   NREJEC=  0   BUL=  F
XP=    .0000000000000000D+00     X=    .1000000000000000D+01
H =    .1266777668108495D+01   K =18   INIAPR=1   IMAX =28   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1266777668108495D+01

YX =    .2980957987041730D+04        DYX =    .1192383194816691D+05
YT =    .2980957987041728D+04        DYT =    .1192383194816691D+05
(YT-YX)/YX = -.6102029654403152D-15    (DYT-DYX)/DYX =  .3051014827201579D-15


           Coefficients AU, ADU and AJK  on   1  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .1839300696370424D+04  .7357202785481693D+04  .2942881114192675D+05
     1      .1283417414302834D+04  .5133669657211334D+04  .2053467862884531D+05
     2      .5558832820675894D+03  .2223533128270356D+04  .8894132513081413D+04
     3      .1716508501676547D+03  .6866034006706182D+03  .2746413602682470D+04
     4      .4093073156462488D+02  .1637229262584994D+03  .6548917050339942D+03
     5      .7927923909155063D+01  .3171169563662009D+02  .1268467825464795D+03
     6      .1291112018849555D+01  .5164448075398153D+01  .2065779230159242D+02
     7      .1812517960576860D+00  .7250071842307584D+00  .2900028736923844D+01
     8      .2234944644573676D-01  .8939778578294545D-01  .3575911431311871D+00
     9      .2456224491795118D-02  .9824897967182089D-02  .3929959186958927D-01
    10      .2434260195714399D-03  .9737040783212100D-03  .3894816312631839D-02
    11      .2196429608241198D-04  .8785718432449386D-04  .3514287367408706D-03
    12      .1818762682150529D-05  .7275050695082988D-05  .2910020235410982D-04
    13      .1391439026149442D-06  .5565755812684584D-06  .2226303376887118D-05
    14      .9891939892834541D-08  .3956775910589151D-07  .1582721281499777D-06
    15      .6567361518380587D-09  .2626947269724252D-08  .1050886695719289D-07
    16      .4089617328584386D-10  .1635899956079975D-09  .6552919665225865D-09
    17      .2397878555651944D-11  .9592179430244616D-11  .3910723828104423D-10
    18      .1326804619373470D-12  .5342538236653549D-12  .3023765265952605D-11
    19      .6879334367125364D-14  .3918617075563577D-13
    20      .3843404585035932D-15
****************************************************************

JSTART (input) =1

RESULTS AFTER 2-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 1332   NACCEP=  2   NREJEC=  0   BUL=  F
XP=    .1000000000000000D+01     X=    .2266777668108495D+01
H =    .1558785559134970D+01   K =18   INIAPR=1   IMAX =28   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1558785559134970D+01

YX =    .4731229983354373D+06        DYX =    .1892491993341750D+07
YT =    .4731229983354323D+06        DYT =    .1892491993341729D+07
(YT-YX)/YX = -.1045743114380780D-13    (DYT-DYX)/DYX = -.1119560275395893D-13


           Coefficients AU, ADU and AJK  on   2  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .2535386626550867D+06  .1014154650620348D+07  .4056618602481394D+07
     1      .1948693511966433D+06  .7794774047865742D+06  .3117909619146296D+07
     2      .9970791863225074D+05  .3988316745290030D+06  .1595326698116010D+07
     3      .3744959058287953D+05  .1497983623315179D+06  .5991934493260708D+06
     4      .1101929190867297D+05  .4407716763469180D+05  .1763086705387671D+06
     5      .2654875816154660D+04  .1061950326461864D+05  .4247801305847469D+05
     6      .5404372406757509D+03  .2161748962703010D+04  .8646995850811998D+04
     7      .9513425637781927D+02  .3805370255112789D+03  .1522148102045129D+04
     8      .1474136569313394D+02  .5896546277254007D+02  .2358618510900778D+03
     9      .2039052291090314D+01  .8156209164361051D+01  .3262483665734912D+02
    10      .2546320840598162D+00  .1018528336240570D+01  .4074113344885602D+01
    11      .2897514440319330D-01  .1159005776135954D+00  .4636023103180876D+00
    12      .3027878792826832D-02  .1211151517034356D-01  .4844606057613965D-01
    13      .2924904316840163D-03  .1169961725251281D-02  .4679846801229903D-02
    14      .2626634882333199D-04  .1050653958704121D-03  .4202615391960540D-03
    15      .2203593575893328D-05  .8814373875164464D-05  .3525737021448094D-04
    16      .1734484746773895D-06  .6937939057049479D-06  .2775151628497952D-05
    17      .1285766942723091D-07  .5142939104825400D-07  .2055919245957938D-06
    18      .9006546218808489D-09  .3600556814288522D-08  .1444730024502405D-07
    19      .5979335629575432D-10  .2387721985797255D-09
    20      .3773600483019617D-11
****************************************************************

JSTART (input) =1

RESULTS AFTER 3-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 1998   NACCEP=  3   NREJEC=  0   BUL=  F
XP=    .2266777668108495D+01     X=    .3825563227243465D+01
H =    .1499768364507135D+01   K =18   INIAPR=1   IMAX =28   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1499768364507135D+01

YX =    .2414693019681911D+09        DYX =    .9658772078727559D+09
YT =    .2414693019681855D+09        DYT =    .9658772078727418D+09
(YT-YX)/YX = -.2344994252895975D-13    (DYT-DYX)/DYX = -.1456364851798566D-13


           Coefficients AU, ADU and AJK  on   3  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .1148186049675955D+09  .4592744198703790D+09  .1837097679481496D+10
     1      .9395775972889505D+08  .3758310389155773D+09  .1503324155662291D+10
     2      .5454234748094279D+08  .2181693899237690D+09  .8726775596950629D+09
     3      .2397719420298554D+08  .9590877681194101D+08  .3836351072477569D+09
     4      .8396434601896802D+07  .3358573840758678D+08  .1343429536303431D+09
     5      .2431101341934393D+07  .9724405367737409D+07  .3889762147094692D+08
     6      .5983724255845899D+06  .2393489702338263D+07  .9573958809351094D+07
     7      .1278759030870270D+06  .5115036123480425D+06  .2046014449390936D+07
     8      .2412325040878489D+05  .9649300163510424D+05  .3859720065396324D+06
     9      .4070545679564563D+04  .1628218271823626D+05  .6512873087243654D+05
    10      .6210368416923072D+03  .2484147366756549D+04  .9936589466749574D+04
    11      .8643870568716390D+02  .3457548227419132D+03  .1383019290784544D+04
    12      .1105828691940950D+02  .4423314767319312D+02  .1769325906064701D+03
    13      .1308687477660078D+01  .5234749906355312D+01  .2093899957918464D+02
    14      .1440613471752191D+00  .5762453882575345D+00  .2304981609672438D+01
    15      .1482197389530176D-01  .5928789714270521D-01  .2371515582614450D+00
    16      .1431331221441803D-02  .5725326337264877D-02  .2290130188694661D-01
    17      .1302167576457617D-03  .5208668507960205D-03  .2083377498792061D-02
    18      .1119761805315144D-04  .4478903178353178D-04  .1791605079972669D-03
    19      .9128448110679586D-06  .3651063735283907D-05
    20      .7072788871607701D-07
****************************************************************

JSTART (input) =1

RESULTS AFTER 4-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 2664   NACCEP=  4   NREJEC=  0   BUL=  F
XP=    .3825563227243465D+01     X=    .5325331591750600D+01
H =    .1516983021034061D+01   K =18   INIAPR=1   IMAX =28   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1516983021034061D+01

YX =    .9732545125420154D+11        DYX =    .3893018050168082D+12
YT =    .9732545125419646D+11        DYT =    .3893018050167858D+12
(YT-YX)/YX = -.5220809862510805D-13    (DYT-DYX)/DYX = -.5753865524148514D-13


           Coefficients AU, ADU and AJK  on   4  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .4730440256562344D+11  .1892176102624946D+12  .7568704410499822D+12
     1      .3831425001125868D+11  .1532570000450355D+12  .6130280001801449D+12
     2      .2175762419771881D+11  .8703049679087582D+11  .3481219871635050D+12
     3      .9299603866393669D+10  .3719841546557498D+11  .1487936618622998D+12
     4      .3155543863567262D+10  .1262217545426907D+11  .5048870181707564D+11
     5      .8835205878155690D+09  .3534082351262247D+10  .1413632940504857D+11
     6      .2100203789926457D+09  .8400815159705670D+09  .3360326063881996D+10
     7      .4330932334474330D+08  .1732372933789652D+09  .6929491735156933D+09
     8      .7878987970015530D+07  .3151595188005761D+08  .1260638075201073D+09
     9      .1281564090840197D+07  .5126256363357104D+07  .2050502545335722D+08
    10      .1884158186599733D+06  .7536632746376917D+06  .3014653098516062D+07
    11      .2526463074338977D+05  .1010585229721487D+06  .4042340918761138D+06
    12      .3113244772707851D+04  .1245297909017443D+05  .4981191636737222D+05
    13      .3548258980961588D+03  .1419303592370974D+04  .5677214378029341D+04
    14      .3761203851103527D+02  .1504481542963784D+03  .6017926243336697D+03
    15      .3725987040294371D+01  .1490394875657188D+02  .5961579415341839D+02
    16      .3464134783858164D+00  .1385653907440029D+01  .5542598958461895D+01
    17      .3033962573737014D-01  .1213575635791763D+00  .4854296418925514D+00
    18      .2511480467026771D-02  .1004511413203067D-01  .4020637659414206D-01
    19      .1970677022405468D-03  .7878823502006402D-03
    20      .1468212863085056D-04
****************************************************************

JSTART (input) =1

RESULTS AFTER 5-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 3330   NACCEP=  5   NREJEC=  0   BUL=  F
XP=    .5325331591750600D+01     X=    .6842314612784661D+01
H =    .1511638178190924D+01   K =18   INIAPR=1   IMAX =28   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1511638178190924D+01

YX =    .4202385088055835D+14        DYX =    .1680954035222321D+15
YT =    .4202385088055630D+14        DYT =    .1680954035222252D+15
(YT-YX)/YX = -.4889336548047213D-13    (DYT-DYX)/DYX = -.4127120584283230D-13


           Coefficients AU, ADU and AJK  on   5  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .2029308896667521D+14  .8117235586670027D+14  .3246894234667985D+15
     1      .1648743083921197D+14  .6594972335684739D+14  .2637988934273874D+15
     2      .9424522469488379D+13  .3769808987795319D+14  .1507923595118115D+15
     3      .4062080866521887D+13  .1624832346608741D+14  .6499329386434926D+14
     4      .1391312848421592D+13  .5565251393686348D+13  .2226100557474521D+14
     5      .3934495657622850D+12  .1573798263049138D+13  .6295193052196369D+13
     6      .9449686463395317D+11  .3779874585358072D+12  .1511949834143076D+13
     7      .1969377552461196D+11  .7877510209844315D+11  .3151004083936783D+12
     8      .3621471329739305D+10  .1448588531895462D+11  .5794354127575281D+11
     9      .5954927909220593D+09  .2381971163686481D+10  .9527884654706583D+10
    10      .8851476794086495D+08  .3540590717623690D+09  .1416236287025485D+10
    11      .1200064424951305D+08  .4800257699713802D+08  .1920103079779711D+09
    12      .1495279314730381D+07  .5981117258540477D+07  .2392446903668686D+08
    13      .1723300704272383D+06  .6893202817573298D+06  .2757281128505383D+07
    14      .1847246539507231D+05  .7388986181723700D+05  .2955594481356256D+06
    15      .1850564740838997D+04  .7402259089611206D+04  .2960902908749133D+05
    16      .1739935630804311D+03  .6959744033535871D+03  .2783893486421555D+04
    17      .1541107453105731D+02  .6164427537488879D+02  .2465634586624801D+03
    18      .1290174160178684D+01  .5160430104530478D+01  .2063861318305135D+02
    19      .1023926069091879D+00  .4092198437171664D+00
    20      .7717210803787389D-02
****************************************************************

JSTART (input) =1

RESULTS AFTER 6-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 3996   NACCEP=  6   NREJEC=  0   BUL=  T
XP=    .6842314612784661D+01     X=    .7000000000000000D+01
H =    .1419168484938049D+01   K =18   INIAPR=1   IMAX =28   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1419168484938049D+01

YX =    .7896296018268428D+14        DYX =    .3158518407307358D+15
YT =    .7896296018268069D+14        DYT =    .3158518407307227D+15
(YT-YX)/YX = -.4551184494205513D-13    (DYT-DYX)/DYX = -.4135641562125896D-13


           Coefficients AU, ADU and AJK  on   6  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .1180923676460555D+15  .4723694705842194D+15  .1889477882336878D+16
     1      .1839370571446434D+14  .7357482285785738D+14  .2942992914314297D+15
     2      .1444236584008648D+13  .5776946336034592D+13  .2310778534413834D+14
     3      .7575489462882669D+11  .3030195785153076D+12  .1212078314061357D+13
     4      .2982654904093234D+10  .1193061961637344D+11  .4772247846541236D+11
     5      .9398636673021084D+08  .3759454669208558D+09  .1503781867759209D+10
     6      .2468585310972509D+07  .9874341243930537D+07  .3949736489275169D+08
     7      .5558387583402714D+05  .2223355032788762D+06  .8893420827369690D+06
     8      .1095217502169747D+04  .4380870066756029D+04  .1752340745925903D+05
     9      .1918356807912148D+02  .7673416686587859D+02  .3069945030212402D+03
    10      .3024279289175070D+00  .1209636298621574D+01  .4789812088012695D+01
    11      .4334005925283078D-02  .1737420467796951D-01  .1464633941650391D+00
    12      .5732037065694029D-04  .2898646870441749D-03 -.5822753906250000D-01
    13      .6697103908737626D-06 -.7432267209996087D-04  .5822753906250000D-01
    14     -.4269282612692931D-06  .6901390980710654D-04 -.3371810913085937D-01
    15      .5441242545589035D-07  .7729557898395931D-04  .3371810913085937D-01
    16      .3584818618270419D-06  .4830973686497457D-04 -.6312942504882812D-01
    17      .3111821717248628D-06 -.6820198145643468D-04  .1411056518554687D-01
    18     -.1716415347198190D-06 -.8588397664884368D-04 -.3371810913085937D-01
    19     -.3586128858993471D-06  .1017047091894202D-04
    20      .7445910851858668D-07
****************************************************************

Из приведенных результатов наглядно видно, что на всех выбранных подпрограммой DE98D элементарных сегментах коэффициенты Чебышёва, как для решения Y, так и для его производных Y', Y'', быстро убывают. Отметим, однако, что в значениях последних коэффициентов Чебышёва (как для решения, так и для его производных), вычисленных на шестом элементарном сегменте, по мере уменьшения их абсолютных величин начинает заметнее проявлять себя вычислительная погрешность. Полезно сделать следующее замечание. При вычислении решения данной задачи Коши (с используемыми здесь вычислительными параметрами) на более длинной разрядной сетке ЭВМ уменьшается вычислительная погрешность, а вместе с ней уменьшается и полная погрешность приближенного решения на каждом элементарном сегменте. Поэтому длины элементарных сегментов будут несколько увеличины вычислительной подпрограммой DE98D и суммарная длина первых пяти сегментов может превысить длину промежутка интегрирования, на котором находится решение задачи Коши. Вследствие этого последний, шестой, (нестандартный) сегмент, который выводит в конец интервала интегрирования будет отрицательным. Такая ситуация складывается при решении задачи Коши из данного примера на языке Pascal с 19-20 значащими цифрами.

2) Последний, шестой, сегмент в первом примере был нестандартный, его длина к тому же в несколько раз меньше длины всех предыдущих сегментов. Поэтому естественно попытаться на шестом сегменте приближать решение (и, соответственно, его производные) частичной суммой мéньшего порядка. Поэтому во втором примере мы используем те же сегменты, что и в первом примере, но на шестом сегменте изменим значение параметра K и возьмем его равным K0 = 12, а параметр IMAX положим равным IMAX = 23. Приводится та же информация, что и в примере 1.

   DOUBLE PRECISION X,YX,DYX,EPS,THRESH,H,HMIN,HMAX,XP,YP,DYP,DELTY,
     1                 RO,YX2,DYX2,AU,ADU,AJK,AU2,ADU2,AJK2,XJ0,XJ02,
     2                 AJKP,U,DU,U2,DU2,ZFI,ZFI2,HD4I,RABC,RABC1,RABC2,
     3                 RABC3,ALP,ALPN,CMAR,P,S,ALP2,ALPN2,CMAR2
C
C
      LOGICAL BUL
C
      PARAMETER (K=18,K2=25,M=1)
      DIMENSION YP(M),DYP(M),DELTY(M),RO(M),YX2(M),DYX2(M),
     1          IU(2),EPS(2),THRESH(2),MEXACT(2),NUMBES(2),
     2          XJ0 (K),  AU (M,K+3),  ADU (M,K+2), AJK (M,K+1), 
     3          XJ02(K2), AU2(M,K2+3), ADU2(M,K2+2),AJK2(M,K2+1),
     4          AJKP(M,K+1), U(M,K), DU(M,K), U2(M,K2),DU2(M,K2),
     5          HD4I(K2+1),RABC(K2),RABC1(K2),RABC2(K2+2),RABC3(K2+2),
     6          ALP (K), ALPN(K),  CMAR (K,K),  P((K+2)*(K-1)/2),
     7          ALP2(K2),ALPN2(K2),CMAR2(K2,K2),S((K+3)*K/2),
     8          ZFI(M,K+4), ZFI2(M,K2+3)
C
      DOUBLE PRECISION WC1,WC2,WC3,HD2,HD4,HOLD,WC30,WC12,WC22,WC32,
     1                 H2D4,H2D8,H2D16,H2D32,H2D96
C
      COMMON/COM70D/WC1,WC2,WC3,LASN,HD2,HD4,HOLD
      COMMON/COM80D/JLAS,JLAS1,JLAS2,JLAS3,H2D4,H2D8,H2D16,H2D32,H2D96
      COMMON/COM75D/LASN0,LASN2,WC30,WC12,WC22,WC32
      COMMON/COM98D/JLASA2,JLAS12,JLAS22,JLAS32,
     1              JLAS0, JLAS10,JLAS20,JLAS30                         25.03.21   
C
      COMMON/STAT76/NACCEP,NREJEC
      COMMON/BLOCK/N,Q
      DOUBLE PRECISION Q,YT,DYT,XK      
C
      EXTERNAL F
C
C
      DATA IU/2,2/, EPS/2*0.5D-11/, THRESH/2*1.D0/, MEXACT/1,1/,
     1     NUMBES/1,1/ 
C
      NREJEC=0
      NACCEP=0
C
      M0=M
      Q=4.0D0
      N=0
      X=0.
      YX=EXP(Q)
      DYX=Q*YX
      XK=7.D0
      K0=K
      INIAPR=1
      IMAX=28
      MESTER=1
      H=1.D0
      HMIN=1.D-3
      HMAX=XK
      K20=K2
      IMAX2=3
      NATTEM=3
C
      JSTART=0
      BUL=.FALSE.
C
      KP1=K0+1
      KP2=K0+2
      KP3=K0+3
C
      DO 20 NX=1,6
      PRINT 54,JSTART
      CALL DE98D(F,M0,K0,INIAPR,IMAX,JSTART,YX,DYX,X,IU,EPS,THRESH,     
     1                 H,HMIN,HMAX,BUL,XP,YP,DYP,MEXACT,NUMBES,DELTY,RO,
     2                 MESTER, K20,IMAX2,YX2,DYX2,NATTEM,AU,ADU,AJK,AU2,
     3                 ADU2,AJK2,XJ0,XJ02,AJKP,U,DU,U2,DU2,ZFI,ZFI2,
     4                 HD4I,RABC,RABC1,RABC2,RABC3,
     5                 ALP,ALPN,CMAR,P,S,ALP2,ALPN2,CMAR2,IERR)         
C
C
C
      YT=EXP(Q*(X+1.D0))
      DYT=Q*YT
      PRINT 44,NX
      PRINT 50,IERR,N,NACCEP,NREJEC,BUL
      PRINT 51,XP,X
      PRINT 35, H,K0,INIAPR,IMAX,JSTART,K20,IMAX2
      PRINT 36, H
      PRINT 52, YX,DYX,YT,DYT,(YT-YX)/YX,(DYT-DYX)/DYX
C
      PRINT 110,NX
      DO 10 L=1,M
      PRINT 47
      PRINT 150,L
      PRINT 85,(J-1,AU(L,J  ),ADU(L,J  ),AJK(L,J),J=1,KP1)
      PRINT 86, KP1,AU(L,KP2),ADU(L,KP2),
     1          KP2,AU(L,KP3)
   10 CONTINUE
C
      PRINT 55
C
      IF (NX .NE. 5) GO TO 20
C
      H=XK-X
      BUL=.TRUE.
C
C  take into consideraitoin: parameter JSTART must be equal to zero    
C
      JSTART=0
C
      K0=  12
      IMAX=23
C
      KP1=K0+1
      KP2=K0+2
      KP3=K0+3
C
   20 CONTINUE
C
      STOP
C
   35 FORMAT(1X,'H =',D25.16,3X,'K =',I2,3X,'INIAPR=',I1,3X,'IMAX =',I2,
     1       3X,'JSTART (exit) =',I1/1X,31X,'K2=',I2,13X,' IMAX2=',I2)
   36 FORMAT(1X,'It is advisable to take the next segment H = ',D22.16/)
   44 FORMAT(/1X,'RESULTS AFTER ',I1,'-nd    CALL SUBROUTINE DE98D'//)
   47 FORMAT(1X,64(1H-)/)
   50 FORMAT(1X,'IERR=',I2,3X,'N=',I5,3X,'NACCEP=',I3,3X,'NREJEC=',I3,
     1           3X,'BUL=',L3)
C
   51 FORMAT(1X,'XP=',D25.16,5X,'X=',D25.16)
C
   52 FORMAT( 1X,'YX =',D25.16,8X,'DYX =',D25.16
     1       /1X,'YT =',D25.16,8X,'DYT =',D25.16
     2       /1X,'(YT-YX)/YX =',D23.16,4X,'(DYT-DYX)/DYX =',D23.16/)
C
   54 FORMAT(1X,'JSTART (input) =',I1)
   55 FORMAT(1X,64(1H*)/)
   85 FORMAT(5X,I2,4X,3D23.16)
   86 FORMAT(5X,I2,4X,2D23.16)
  110 FORMAT(/12X,'Coefficients AU, ADU and AJK  on ',I3,'  segment  ')
C      
  150 FORMAT('   Number of',15X,'Chebyshev coefficients'
     *       ' for ',I1,' component'/  
     1        14X,72(1H-)/
     2       ' coefficient',9X,'for Y',18X,'for Y''',18X,'for Y'''''/         
     3        1X,85(1H-)/)
C   
      END
C         
C
      SUBROUTINE F(X,Y,DY,Z,M)
      DOUBLE PRECISION X,Y,DY,Z,Q
      COMMON/BLOCK/N,Q
      Z=Q*DY
      N=N+1
      RETURN
      END

Результаты:


JSTART (input) =0

RESULTS AFTER 1-nd    CALL SUBROUTINE DE98D


IERR= 0   N=  666   NACCEP=  1   NREJEC=  0   BUL=  F
XP=    .0000000000000000D+00     X=    .1000000000000000D+01
H =    .1266777668108495D+01   K =18   INIAPR=1   IMAX =28   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1266777668108495D+01

YX =    .2980957987041730D+04        DYX =    .1192383194816691D+05
YT =    .2980957987041728D+04        DYT =    .1192383194816691D+05
(YT-YX)/YX = -.6102029654403152D-15    (DYT-DYX)/DYX =  .3051014827201579D-15


           Coefficients AU, ADU and AJK  on   1  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .1839300696370424D+04  .7357202785481693D+04  .2942881114192675D+05
     1      .1283417414302834D+04  .5133669657211334D+04  .2053467862884531D+05
     2      .5558832820675894D+03  .2223533128270356D+04  .8894132513081413D+04
     3      .1716508501676547D+03  .6866034006706182D+03  .2746413602682470D+04
     4      .4093073156462488D+02  .1637229262584994D+03  .6548917050339942D+03
     5      .7927923909155063D+01  .3171169563662009D+02  .1268467825464795D+03
     6      .1291112018849555D+01  .5164448075398153D+01  .2065779230159242D+02
     7      .1812517960576860D+00  .7250071842307584D+00  .2900028736923844D+01
     8      .2234944644573676D-01  .8939778578294545D-01  .3575911431311871D+00
     9      .2456224491795118D-02  .9824897967182089D-02  .3929959186958927D-01
    10      .2434260195714399D-03  .9737040783212100D-03  .3894816312631839D-02
    11      .2196429608241198D-04  .8785718432449386D-04  .3514287367408706D-03
    12      .1818762682150529D-05  .7275050695082988D-05  .2910020235410982D-04
    13      .1391439026149442D-06  .5565755812684584D-06  .2226303376887118D-05
    14      .9891939892834541D-08  .3956775910589151D-07  .1582721281499777D-06
    15      .6567361518380587D-09  .2626947269724252D-08  .1050886695719289D-07
    16      .4089617328584386D-10  .1635899956079975D-09  .6552919665225865D-09
    17      .2397878555651944D-11  .9592179430244616D-11  .3910723828104423D-10
    18      .1326804619373470D-12  .5342538236653549D-12  .3023765265952605D-11
    19      .6879334367125364D-14  .3918617075563577D-13
    20      .3843404585035932D-15
****************************************************************

JSTART (input) =1

RESULTS AFTER 2-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 1332   NACCEP=  2   NREJEC=  0   BUL=  F
XP=    .1000000000000000D+01     X=    .2266777668108495D+01
H =    .1558785559134970D+01   K =18   INIAPR=1   IMAX =28   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1558785559134970D+01

YX =    .4731229983354373D+06        DYX =    .1892491993341750D+07
YT =    .4731229983354323D+06        DYT =    .1892491993341729D+07
(YT-YX)/YX = -.1045743114380780D-13    (DYT-DYX)/DYX = -.1119560275395893D-13


           Coefficients AU, ADU and AJK  on   2  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .2535386626550867D+06  .1014154650620348D+07  .4056618602481394D+07
     1      .1948693511966433D+06  .7794774047865742D+06  .3117909619146296D+07
     2      .9970791863225074D+05  .3988316745290030D+06  .1595326698116010D+07
     3      .3744959058287953D+05  .1497983623315179D+06  .5991934493260708D+06
     4      .1101929190867297D+05  .4407716763469180D+05  .1763086705387671D+06
     5      .2654875816154660D+04  .1061950326461864D+05  .4247801305847469D+05
     6      .5404372406757509D+03  .2161748962703010D+04  .8646995850811998D+04
     7      .9513425637781927D+02  .3805370255112789D+03  .1522148102045129D+04
     8      .1474136569313394D+02  .5896546277254007D+02  .2358618510900778D+03
     9      .2039052291090314D+01  .8156209164361051D+01  .3262483665734912D+02
    10      .2546320840598162D+00  .1018528336240570D+01  .4074113344885602D+01
    11      .2897514440319330D-01  .1159005776135954D+00  .4636023103180876D+00
    12      .3027878792826832D-02  .1211151517034356D-01  .4844606057613965D-01
    13      .2924904316840163D-03  .1169961725251281D-02  .4679846801229903D-02
    14      .2626634882333199D-04  .1050653958704121D-03  .4202615391960540D-03
    15      .2203593575893328D-05  .8814373875164464D-05  .3525737021448094D-04
    16      .1734484746773895D-06  .6937939057049479D-06  .2775151628497952D-05
    17      .1285766942723091D-07  .5142939104825400D-07  .2055919245957938D-06
    18      .9006546218808489D-09  .3600556814288522D-08  .1444730024502405D-07
    19      .5979335629575432D-10  .2387721985797255D-09
    20      .3773600483019617D-11
****************************************************************

JSTART (input) =1

RESULTS AFTER 3-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 1998   NACCEP=  3   NREJEC=  0   BUL=  F
XP=    .2266777668108495D+01     X=    .3825563227243465D+01
H =    .1499768364507135D+01   K =18   INIAPR=1   IMAX =28   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1499768364507135D+01

YX =    .2414693019681911D+09        DYX =    .9658772078727559D+09
YT =    .2414693019681855D+09        DYT =    .9658772078727418D+09
(YT-YX)/YX = -.2344994252895975D-13    (DYT-DYX)/DYX = -.1456364851798566D-13


           Coefficients AU, ADU and AJK  on   3  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .1148186049675955D+09  .4592744198703790D+09  .1837097679481496D+10
     1      .9395775972889505D+08  .3758310389155773D+09  .1503324155662291D+10
     2      .5454234748094279D+08  .2181693899237690D+09  .8726775596950629D+09
     3      .2397719420298554D+08  .9590877681194101D+08  .3836351072477569D+09
     4      .8396434601896802D+07  .3358573840758678D+08  .1343429536303431D+09
     5      .2431101341934393D+07  .9724405367737409D+07  .3889762147094692D+08
     6      .5983724255845899D+06  .2393489702338263D+07  .9573958809351094D+07
     7      .1278759030870270D+06  .5115036123480425D+06  .2046014449390936D+07
     8      .2412325040878489D+05  .9649300163510424D+05  .3859720065396324D+06
     9      .4070545679564563D+04  .1628218271823626D+05  .6512873087243654D+05
    10      .6210368416923072D+03  .2484147366756549D+04  .9936589466749574D+04
    11      .8643870568716390D+02  .3457548227419132D+03  .1383019290784544D+04
    12      .1105828691940950D+02  .4423314767319312D+02  .1769325906064701D+03
    13      .1308687477660078D+01  .5234749906355312D+01  .2093899957918464D+02
    14      .1440613471752191D+00  .5762453882575345D+00  .2304981609672438D+01
    15      .1482197389530176D-01  .5928789714270521D-01  .2371515582614450D+00
    16      .1431331221441803D-02  .5725326337264877D-02  .2290130188694661D-01
    17      .1302167576457617D-03  .5208668507960205D-03  .2083377498792061D-02
    18      .1119761805315144D-04  .4478903178353178D-04  .1791605079972669D-03
    19      .9128448110679586D-06  .3651063735283907D-05
    20      .7072788871607701D-07
****************************************************************

JSTART (input) =1

RESULTS AFTER 4-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 2664   NACCEP=  4   NREJEC=  0   BUL=  F
XP=    .3825563227243465D+01     X=    .5325331591750600D+01
H =    .1516983021034061D+01   K =18   INIAPR=1   IMAX =28   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1516983021034061D+01

YX =    .9732545125420154D+11        DYX =    .3893018050168082D+12
YT =    .9732545125419646D+11        DYT =    .3893018050167858D+12
(YT-YX)/YX = -.5220809862510805D-13    (DYT-DYX)/DYX = -.5753865524148514D-13


           Coefficients AU, ADU and AJK  on   4  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .4730440256562344D+11  .1892176102624946D+12  .7568704410499822D+12
     1      .3831425001125868D+11  .1532570000450355D+12  .6130280001801449D+12
     2      .2175762419771881D+11  .8703049679087582D+11  .3481219871635050D+12
     3      .9299603866393669D+10  .3719841546557498D+11  .1487936618622998D+12
     4      .3155543863567262D+10  .1262217545426907D+11  .5048870181707564D+11
     5      .8835205878155690D+09  .3534082351262247D+10  .1413632940504857D+11
     6      .2100203789926457D+09  .8400815159705670D+09  .3360326063881996D+10
     7      .4330932334474330D+08  .1732372933789652D+09  .6929491735156933D+09
     8      .7878987970015530D+07  .3151595188005761D+08  .1260638075201073D+09
     9      .1281564090840197D+07  .5126256363357104D+07  .2050502545335722D+08
    10      .1884158186599733D+06  .7536632746376917D+06  .3014653098516062D+07
    11      .2526463074338977D+05  .1010585229721487D+06  .4042340918761138D+06
    12      .3113244772707851D+04  .1245297909017443D+05  .4981191636737222D+05
    13      .3548258980961588D+03  .1419303592370974D+04  .5677214378029341D+04
    14      .3761203851103527D+02  .1504481542963784D+03  .6017926243336697D+03
    15      .3725987040294371D+01  .1490394875657188D+02  .5961579415341839D+02
    16      .3464134783858164D+00  .1385653907440029D+01  .5542598958461895D+01
    17      .3033962573737014D-01  .1213575635791763D+00  .4854296418925514D+00
    18      .2511480467026771D-02  .1004511413203067D-01  .4020637659414206D-01
    19      .1970677022405468D-03  .7878823502006402D-03
    20      .1468212863085056D-04
****************************************************************

JSTART (input) =1

RESULTS AFTER 5-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 3330   NACCEP=  5   NREJEC=  0   BUL=  F
XP=    .5325331591750600D+01     X=    .6842314612784661D+01
H =    .1511638178190924D+01   K =18   INIAPR=1   IMAX =28   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1511638178190924D+01

YX =    .4202385088055835D+14        DYX =    .1680954035222321D+15
YT =    .4202385088055630D+14        DYT =    .1680954035222252D+15
(YT-YX)/YX = -.4889336548047213D-13    (DYT-DYX)/DYX = -.4127120584283230D-13


           Coefficients AU, ADU and AJK  on   5  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .2029308896667521D+14  .8117235586670027D+14  .3246894234667985D+15
     1      .1648743083921197D+14  .6594972335684739D+14  .2637988934273874D+15
     2      .9424522469488379D+13  .3769808987795319D+14  .1507923595118115D+15
     3      .4062080866521887D+13  .1624832346608741D+14  .6499329386434926D+14
     4      .1391312848421592D+13  .5565251393686348D+13  .2226100557474521D+14
     5      .3934495657622850D+12  .1573798263049138D+13  .6295193052196369D+13
     6      .9449686463395317D+11  .3779874585358072D+12  .1511949834143076D+13
     7      .1969377552461196D+11  .7877510209844315D+11  .3151004083936783D+12
     8      .3621471329739305D+10  .1448588531895462D+11  .5794354127575281D+11
     9      .5954927909220593D+09  .2381971163686481D+10  .9527884654706583D+10
    10      .8851476794086495D+08  .3540590717623690D+09  .1416236287025485D+10
    11      .1200064424951305D+08  .4800257699713802D+08  .1920103079779711D+09
    12      .1495279314730381D+07  .5981117258540477D+07  .2392446903668686D+08
    13      .1723300704272383D+06  .6893202817573298D+06  .2757281128505383D+07
    14      .1847246539507231D+05  .7388986181723700D+05  .2955594481356256D+06
    15      .1850564740838997D+04  .7402259089611206D+04  .2960902908749133D+05
    16      .1739935630804311D+03  .6959744033535871D+03  .2783893486421555D+04
    17      .1541107453105731D+02  .6164427537488879D+02  .2465634586624801D+03
    18      .1290174160178684D+01  .5160430104530478D+01  .2063861318305135D+02
    19      .1023926069091879D+00  .4092198437171664D+00
    20      .7717210803787389D-02
****************************************************************

JSTART (input) =0

RESULTS AFTER 6-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 3756   NACCEP=  6   NREJEC=  0   BUL=  T
XP=    .6842314612784661D+01     X=    .7000000000000000D+01
H =    .2789580198714488D+00   K =12   INIAPR=1   IMAX =23   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .2789580198714488D+00

YX =    .7896296018268428D+14        DYX =    .3158518407307358D+15
YT =    .7896296018268069D+14        DYT =    .3158518407307227D+15
(YT-YX)/YX = -.4551184494205513D-13    (DYT-DYX)/DYX = -.4135641562125896D-13


           Coefficients AU, ADU and AJK  on   6  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .1180923676460555D+15  .4723694705842194D+15  .1889477882336878D+16
     1      .1839370571446434D+14  .7357482285785738D+14  .2942992914314297D+15
     2      .1444236584008648D+13  .5776946336034592D+13  .2310778534413834D+14
     3      .7575489462882669D+11  .3030195785153076D+12  .1212078314061357D+13
     4      .2982654904093234D+10  .1193061961637344D+11  .4772247846541236D+11
     5      .9398636673021084D+08  .3759454669208558D+09  .1503781867759209D+10
     6      .2468585310972509D+07  .9874341243930537D+07  .3949736489275169D+08
     7      .5558387583402714D+05  .2223355032788762D+06  .8893420827369690D+06
     8      .1095217502169747D+04  .4380870066756029D+04  .1752340745925903D+05
     9      .1918356807912148D+02  .7673416686587859D+02  .3069945030212402D+03
    10      .3024279289175070D+00  .1209636298621574D+01  .4789812088012695D+01
    11      .4334005925283078D-02  .1737420467796951D-01  .1464633941650391D+00
    12      .5732037065694029D-04  .2898646870441749D-03 -.5822753906250000D-01
    13      .6697103908737626D-06 -.7432267209996087D-04
    14     -.4269282612692931D-06
****************************************************************

Результаты второго примера показывают, что в конце x_f интервала интегрирования получены те же значения для решения YX и его производной DYX, какие были вычислены в первом примере, а также те же значения коэффициентов Чебышёва (как в первом примере) для решения Y (коэффициента с номерами i = 0, ... , 14), для первой производной Y' (коэффициенты с номерами i = 0, ... , 13) и для второй производной Y'' (коэффициенты с номерами i = 0, ... , 12). Как следствие, при этом используется меньшее число обращений к правой части дифференциального уравнения.

3) Третий пример только для Фортрана. Данный пример отличается от первого примера только тем, что шестое обращение к подпрограмме DE98D выполняется с тем значением длины частичного сегмента H, которое было выбрано и рекомендовано подпрограммой DE98D после вычисления решения на пятом элементарном сегменте. В результате точка XK = 7, которая является концом промежутка интегрирования, оказалась промежуточной точкой шестого элементарного сегмента. Чтобы найти искомое решение задачи Коши и его производную в этой точке, достаточно вычислить значения соответствующих частичных сумм, коэффициенты которых получены на шестом элементарном сегменте, в точке XK = 7. Такой способ иллюстрируется в третьем примере на языке Фортран. Приводятся вызывающая программа на Фортране и результаты, которые (помимо уже выше указанных) включают приближенные значения решения Y и его производной DY в точке XK, точные значения решения YT и его производной DYT в точке XK и их относительные погрешности.

      DOUBLE PRECISION X,YX,DYX,EPS,THRESH,H,HMIN,HMAX,XP,YP,DYP,DELTY,
     1                 RO,YX2,DYX2,AU,ADU,AJK,AU2,ADU2,AJK2,XJ0,XJ02,
     2                 AJKP,U,DU,U2,DU2,ZFI,ZFI2,HD4I,RABC,RABC1,RABC2,
     3                 RABC3,ALP,ALPN,CMAR,P,S,ALP2,ALPN2,CMAR2
C
C
      DOUBLE PRECISION T,DY,Y
      LOGICAL BUL
C
      PARAMETER (K=18,K2=25,M=1)
      DIMENSION YP(M),DYP(M),DELTY(M),RO(M),YX2(M),DYX2(M),
     1          IU(2),EPS(2),THRESH(2),MEXACT(2),NUMBES(2),
     2          XJ0 (K),  AU (M,K+3),  ADU (M,K+2), AJK (M,K+1), 
     3          XJ02(K2), AU2(M,K2+3), ADU2(M,K2+2),AJK2(M,K2+1),
     4          AJKP(M,K+1), U(M,K), DU(M,K), U2(M,K2),DU2(M,K2),
     5          HD4I(K2+1),RABC(K2),RABC1(K2),RABC2(K2+2),RABC3(K2+2),
     6          ALP (K), ALPN(K),  CMAR (K,K),  P((K+2)*(K-1)/2),
     7          ALP2(K2),ALPN2(K2),CMAR2(K2,K2),S((K+3)*K/2),
     8          ZFI(M,K+4), ZFI2(M,K2+3)
C
      DOUBLE PRECISION WC1,WC2,WC3,HD2,HD4,HOLD,WC30,WC12,WC22,WC32,
     1                 H2D4,H2D8,H2D16,H2D32,H2D96
C
      COMMON/COM70D/WC1,WC2,WC3,LASN,HD2,HD4,HOLD
      COMMON/COM80D/JLAS,JLAS1,JLAS2,JLAS3,H2D4,H2D8,H2D16,H2D32,H2D96
      COMMON/COM75D/LASN0,LASN2,WC30,WC12,WC22,WC32
      COMMON/COM98D/JLASA2,JLAS12,JLAS22,JLAS32,
     1              JLAS0, JLAS10,JLAS20,JLAS30                         25.03.21   
C
      COMMON/STAT76/NACCEP,NREJEC
      COMMON/BLOCK/N,Q
      DOUBLE PRECISION Q,YT,DYT,XK      
C
      EXTERNAL F
C
C
      DATA IU/2,2/, EPS/2*0.5D-11/, THRESH/2*1.D0/, MEXACT/1,1/,
     1     NUMBES/1,1/ 
C
      NREJEC=0
      NACCEP=0
C
      M0=M
      Q=4.0D0
      N=0
      X=0.
      YX=EXP(Q)
      DYX=Q*YX
      XK=7.D0
      K0=K
      INIAPR=1
      IMAX=28
      MESTER=1
      H=1.D0
      HMIN=1.D-3
      HMAX=XK
      K20=K2
      IMAX2=3
      NATTEM=3
C
      JSTART=0
      BUL=.FALSE.
C
      KP1=K0+1
      KP2=K0+2
      KP3=K0+3
C
      DO 20 NX=1,6
      PRINT 54,JSTART
      CALL DE98D(F,M0,K0,INIAPR,IMAX,JSTART,YX,DYX,X,IU,EPS,THRESH,     
     1                 H,HMIN,HMAX,BUL,XP,YP,DYP,MEXACT,NUMBES,DELTY,RO,
     2                 MESTER, K20,IMAX2,YX2,DYX2,NATTEM,AU,ADU,AJK,AU2,
     3                 ADU2,AJK2,XJ0,XJ02,AJKP,U,DU,U2,DU2,ZFI,ZFI2,
     4                 HD4I,RABC,RABC1,RABC2,RABC3,
     5                 ALP,ALPN,CMAR,P,S,ALP2,ALPN2,CMAR2,IERR)         
C
C
C
      YT=EXP(Q*(X+1.D0))
      DYT=Q*YT
      PRINT 44,NX
      PRINT 50,IERR,N,NACCEP,NREJEC,BUL
      PRINT 51,XP,X
      PRINT 35, H,K0,INIAPR,IMAX,JSTART,K20,IMAX2
      PRINT 36, H
      PRINT 52, YX,DYX,YT,DYT,(YT-YX)/YX,(DYT-DYX)/DYX
C
      PRINT 110,NX
      DO 10 L=1,M
      PRINT 47
      PRINT 150,L
      PRINT 85,(J-1,AU(L,J  ),ADU(L,J  ),AJK(L,J),J=1,KP1)
      PRINT 86, KP1,AU(L,KP2),ADU(L,KP2),
     1          KP2,AU(L,KP3)
   10 CONTINUE
C
      PRINT 55
C
   20 CONTINUE
C
C        Calculation value of partial sum for derivative on sixth segment
C                          at point XK
C      For this map point XK onto interval  [-1, 1]
C
      T=2.D0*(XK-XP)/(X-XP)-1.D0
      CALL DE70DC(M,KP1,T,ADU,DY,ZFI(1,KP2))
C
C        Calculation value of partial sum for solution on sixth segment
C                          at point XK
C
      CALL DE70DC(M,KP2,T, AU, Y,ZFI(1,KP2)) 
C
      YT=EXP(Q*(XK+1.D0))
      DYT=Q*YT
      PRINT 56, XK
      PRINT 53, Y,DY,YT,DYT,(YT-Y)/Y,(DYT-DY)/DY
C      
      STOP
C
   35 FORMAT(1X,'H =',D25.16,3X,'K =',I2,3X,'INIAPR=',I1,3X,'IMAX =',I2,
     1       3X,'JSTART (exit) =',I1/1X,31X,'K2=',I2,13X,' IMAX2=',I2)
   36 FORMAT(1X,'It is advisable to take the next segment H = ',D22.16/)
   44 FORMAT(/1X,'RESULTS AFTER ',I1,'-nd    CALL SUBROUTINE DE98D'//)
   47 FORMAT(1X,64(1H-)/)
   50 FORMAT(1X,'IERR=',I2,3X,'N=',I5,3X,'NACCEP=',I3,3X,'NREJEC=',I3,
     1           3X,'BUL=',L3)
C
   51 FORMAT(1X,'XP=',D25.16,5X,'X=',D25.16)
C
   52 FORMAT( 1X,'YX =',D25.16,8X,'DYX =',D25.16
     1       /1X,'YT =',D25.16,8X,'DYT =',D25.16
     2       /1X,'(YT-YX)/YX =',D23.16,4X,'(DYT-DYX)/DYX =',D23.16/)
C
   53 FORMAT( 1X,'Y  =',D25.16,8X,'DY  =',D25.16
     1       /1X,'YT =',D25.16,8X,'DYT =',D25.16
     2       /1X,'(YT-Y)/Y  =',D23.16,4X,'(DYT-DY)/DY  =',D23.16/)
   54 FORMAT(1X,'JSTART (input) =',I1)
   55 FORMAT(1X,64(1H*)/)
   56 FORMAT(1X,'Results at point XK=',D25.16/)
   85 FORMAT(5X,I2,4X,3D23.16)
   86 FORMAT(5X,I2,4X,2D23.16)
  110 FORMAT(/12X,'Coefficients AU, ADU and AJK  on ',I3,'  segment  ')
C      
  150 FORMAT('   Number of',15X,'Chebyshev coefficients'
     *       ' for ',I1,' component'/  
     1        14X,72(1H-)/
     2       ' coefficient',9X,'for Y',18X,'for Y''',18X,'for Y'''''/         
     3        1X,85(1H-)/)
C   
      END
C         
C
      SUBROUTINE F(X,Y,DY,Z,M)
      DOUBLE PRECISION X,Y,DY,Z,Q
      COMMON/BLOCK/N,Q
      Z=Q*DY
      N=N+1
      RETURN
      END

Результаты:


JSTART (input) =0

RESULTS AFTER 1-nd    CALL SUBROUTINE DE98D


IERR= 0   N=  666   NACCEP=  1   NREJEC=  0   BUL=  F
XP=    .0000000000000000D+00     X=    .1000000000000000D+01
H =    .1266777668108495D+01   K =18   INIAPR=1   IMAX =28   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1266777668108495D+01

YX =    .2980957987041730D+04        DYX =    .1192383194816691D+05
YT =    .2980957987041728D+04        DYT =    .1192383194816691D+05
(YT-YX)/YX = -.6102029654403152D-15    (DYT-DYX)/DYX =  .3051014827201579D-15


           Coefficients AU, ADU and AJK  on   1  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .1839300696370424D+04  .7357202785481693D+04  .2942881114192675D+05
     1      .1283417414302834D+04  .5133669657211334D+04  .2053467862884531D+05
     2      .5558832820675894D+03  .2223533128270356D+04  .8894132513081413D+04
     3      .1716508501676547D+03  .6866034006706182D+03  .2746413602682470D+04
     4      .4093073156462488D+02  .1637229262584994D+03  .6548917050339942D+03
     5      .7927923909155063D+01  .3171169563662009D+02  .1268467825464795D+03
     6      .1291112018849555D+01  .5164448075398153D+01  .2065779230159242D+02
     7      .1812517960576860D+00  .7250071842307584D+00  .2900028736923844D+01
     8      .2234944644573676D-01  .8939778578294545D-01  .3575911431311871D+00
     9      .2456224491795118D-02  .9824897967182089D-02  .3929959186958927D-01
    10      .2434260195714399D-03  .9737040783212100D-03  .3894816312631839D-02
    11      .2196429608241198D-04  .8785718432449386D-04  .3514287367408706D-03
    12      .1818762682150529D-05  .7275050695082988D-05  .2910020235410982D-04
    13      .1391439026149442D-06  .5565755812684584D-06  .2226303376887118D-05
    14      .9891939892834541D-08  .3956775910589151D-07  .1582721281499777D-06
    15      .6567361518380587D-09  .2626947269724252D-08  .1050886695719289D-07
    16      .4089617328584386D-10  .1635899956079975D-09  .6552919665225865D-09
    17      .2397878555651944D-11  .9592179430244616D-11  .3910723828104423D-10
    18      .1326804619373470D-12  .5342538236653549D-12  .3023765265952605D-11
    19      .6879334367125364D-14  .3918617075563577D-13
    20      .3843404585035932D-15
****************************************************************

JSTART (input) =1

RESULTS AFTER 2-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 1332   NACCEP=  2   NREJEC=  0   BUL=  F
XP=    .1000000000000000D+01     X=    .2266777668108495D+01
H =    .1558785559134970D+01   K =18   INIAPR=1   IMAX =28   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1558785559134970D+01

YX =    .4731229983354373D+06        DYX =    .1892491993341750D+07
YT =    .4731229983354323D+06        DYT =    .1892491993341729D+07
(YT-YX)/YX = -.1045743114380780D-13    (DYT-DYX)/DYX = -.1119560275395893D-13


           Coefficients AU, ADU and AJK  on   2  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .2535386626550867D+06  .1014154650620348D+07  .4056618602481394D+07
     1      .1948693511966433D+06  .7794774047865742D+06  .3117909619146296D+07
     2      .9970791863225074D+05  .3988316745290030D+06  .1595326698116010D+07
     3      .3744959058287953D+05  .1497983623315179D+06  .5991934493260708D+06
     4      .1101929190867297D+05  .4407716763469180D+05  .1763086705387671D+06
     5      .2654875816154660D+04  .1061950326461864D+05  .4247801305847469D+05
     6      .5404372406757509D+03  .2161748962703010D+04  .8646995850811998D+04
     7      .9513425637781927D+02  .3805370255112789D+03  .1522148102045129D+04
     8      .1474136569313394D+02  .5896546277254007D+02  .2358618510900778D+03
     9      .2039052291090314D+01  .8156209164361051D+01  .3262483665734912D+02
    10      .2546320840598162D+00  .1018528336240570D+01  .4074113344885602D+01
    11      .2897514440319330D-01  .1159005776135954D+00  .4636023103180876D+00
    12      .3027878792826832D-02  .1211151517034356D-01  .4844606057613965D-01
    13      .2924904316840163D-03  .1169961725251281D-02  .4679846801229903D-02
    14      .2626634882333199D-04  .1050653958704121D-03  .4202615391960540D-03
    15      .2203593575893328D-05  .8814373875164464D-05  .3525737021448094D-04
    16      .1734484746773895D-06  .6937939057049479D-06  .2775151628497952D-05
    17      .1285766942723091D-07  .5142939104825400D-07  .2055919245957938D-06
    18      .9006546218808489D-09  .3600556814288522D-08  .1444730024502405D-07
    19      .5979335629575432D-10  .2387721985797255D-09
    20      .3773600483019617D-11
****************************************************************

JSTART (input) =1

RESULTS AFTER 3-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 1998   NACCEP=  3   NREJEC=  0   BUL=  F
XP=    .2266777668108495D+01     X=    .3825563227243465D+01
H =    .1499768364507135D+01   K =18   INIAPR=1   IMAX =28   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1499768364507135D+01

YX =    .2414693019681911D+09        DYX =    .9658772078727559D+09
YT =    .2414693019681855D+09        DYT =    .9658772078727418D+09
(YT-YX)/YX = -.2344994252895975D-13    (DYT-DYX)/DYX = -.1456364851798566D-13


           Coefficients AU, ADU and AJK  on   3  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .1148186049675955D+09  .4592744198703790D+09  .1837097679481496D+10
     1      .9395775972889505D+08  .3758310389155773D+09  .1503324155662291D+10
     2      .5454234748094279D+08  .2181693899237690D+09  .8726775596950629D+09
     3      .2397719420298554D+08  .9590877681194101D+08  .3836351072477569D+09
     4      .8396434601896802D+07  .3358573840758678D+08  .1343429536303431D+09
     5      .2431101341934393D+07  .9724405367737409D+07  .3889762147094692D+08
     6      .5983724255845899D+06  .2393489702338263D+07  .9573958809351094D+07
     7      .1278759030870270D+06  .5115036123480425D+06  .2046014449390936D+07
     8      .2412325040878489D+05  .9649300163510424D+05  .3859720065396324D+06
     9      .4070545679564563D+04  .1628218271823626D+05  .6512873087243654D+05
    10      .6210368416923072D+03  .2484147366756549D+04  .9936589466749574D+04
    11      .8643870568716390D+02  .3457548227419132D+03  .1383019290784544D+04
    12      .1105828691940950D+02  .4423314767319312D+02  .1769325906064701D+03
    13      .1308687477660078D+01  .5234749906355312D+01  .2093899957918464D+02
    14      .1440613471752191D+00  .5762453882575345D+00  .2304981609672438D+01
    15      .1482197389530176D-01  .5928789714270521D-01  .2371515582614450D+00
    16      .1431331221441803D-02  .5725326337264877D-02  .2290130188694661D-01
    17      .1302167576457617D-03  .5208668507960205D-03  .2083377498792061D-02
    18      .1119761805315144D-04  .4478903178353178D-04  .1791605079972669D-03
    19      .9128448110679586D-06  .3651063735283907D-05
    20      .7072788871607701D-07
****************************************************************

JSTART (input) =1

RESULTS AFTER 4-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 2664   NACCEP=  4   NREJEC=  0   BUL=  F
XP=    .3825563227243465D+01     X=    .5325331591750600D+01
H =    .1516983021034061D+01   K =18   INIAPR=1   IMAX =28   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1516983021034061D+01

YX =    .9732545125420154D+11        DYX =    .3893018050168082D+12
YT =    .9732545125419646D+11        DYT =    .3893018050167858D+12
(YT-YX)/YX = -.5220809862510805D-13    (DYT-DYX)/DYX = -.5753865524148514D-13


           Coefficients AU, ADU and AJK  on   4  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .4730440256562344D+11  .1892176102624946D+12  .7568704410499822D+12
     1      .3831425001125868D+11  .1532570000450355D+12  .6130280001801449D+12
     2      .2175762419771881D+11  .8703049679087582D+11  .3481219871635050D+12
     3      .9299603866393669D+10  .3719841546557498D+11  .1487936618622998D+12
     4      .3155543863567262D+10  .1262217545426907D+11  .5048870181707564D+11
     5      .8835205878155690D+09  .3534082351262247D+10  .1413632940504857D+11
     6      .2100203789926457D+09  .8400815159705670D+09  .3360326063881996D+10
     7      .4330932334474330D+08  .1732372933789652D+09  .6929491735156933D+09
     8      .7878987970015530D+07  .3151595188005761D+08  .1260638075201073D+09
     9      .1281564090840197D+07  .5126256363357104D+07  .2050502545335722D+08
    10      .1884158186599733D+06  .7536632746376917D+06  .3014653098516062D+07
    11      .2526463074338977D+05  .1010585229721487D+06  .4042340918761138D+06
    12      .3113244772707851D+04  .1245297909017443D+05  .4981191636737222D+05
    13      .3548258980961588D+03  .1419303592370974D+04  .5677214378029341D+04
    14      .3761203851103527D+02  .1504481542963784D+03  .6017926243336697D+03
    15      .3725987040294371D+01  .1490394875657188D+02  .5961579415341839D+02
    16      .3464134783858164D+00  .1385653907440029D+01  .5542598958461895D+01
    17      .3033962573737014D-01  .1213575635791763D+00  .4854296418925514D+00
    18      .2511480467026771D-02  .1004511413203067D-01  .4020637659414206D-01
    19      .1970677022405468D-03  .7878823502006402D-03
    20      .1468212863085056D-04
****************************************************************

JSTART (input) =1

RESULTS AFTER 5-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 3330   NACCEP=  5   NREJEC=  0   BUL=  F
XP=    .5325331591750600D+01     X=    .6842314612784661D+01
H =    .1511638178190924D+01   K =18   INIAPR=1   IMAX =28   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1511638178190924D+01

YX =    .4202385088055835D+14        DYX =    .1680954035222321D+15
YT =    .4202385088055630D+14        DYT =    .1680954035222252D+15
(YT-YX)/YX = -.4889336548047213D-13    (DYT-DYX)/DYX = -.4127120584283230D-13


           Coefficients AU, ADU and AJK  on   5  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .2029308896667521D+14  .8117235586670027D+14  .3246894234667985D+15
     1      .1648743083921197D+14  .6594972335684739D+14  .2637988934273874D+15
     2      .9424522469488379D+13  .3769808987795319D+14  .1507923595118115D+15
     3      .4062080866521887D+13  .1624832346608741D+14  .6499329386434926D+14
     4      .1391312848421592D+13  .5565251393686348D+13  .2226100557474521D+14
     5      .3934495657622850D+12  .1573798263049138D+13  .6295193052196369D+13
     6      .9449686463395317D+11  .3779874585358072D+12  .1511949834143076D+13
     7      .1969377552461196D+11  .7877510209844315D+11  .3151004083936783D+12
     8      .3621471329739305D+10  .1448588531895462D+11  .5794354127575281D+11
     9      .5954927909220593D+09  .2381971163686481D+10  .9527884654706583D+10
    10      .8851476794086495D+08  .3540590717623690D+09  .1416236287025485D+10
    11      .1200064424951305D+08  .4800257699713802D+08  .1920103079779711D+09
    12      .1495279314730381D+07  .5981117258540477D+07  .2392446903668686D+08
    13      .1723300704272383D+06  .6893202817573298D+06  .2757281128505383D+07
    14      .1847246539507231D+05  .7388986181723700D+05  .2955594481356256D+06
    15      .1850564740838997D+04  .7402259089611206D+04  .2960902908749133D+05
    16      .1739935630804311D+03  .6959744033535871D+03  .2783893486421555D+04
    17      .1541107453105731D+02  .6164427537488879D+02  .2465634586624801D+03
    18      .1290174160178684D+01  .5160430104530478D+01  .2063861318305135D+02
    19      .1023926069091879D+00  .4092198437171664D+00
    20      .7717210803787389D-02
****************************************************************

JSTART (input) =1

RESULTS AFTER 6-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 3996   NACCEP=  6   NREJEC=  0   BUL=  F
XP=    .6842314612784661D+01     X=    .8353952790975585D+01
H =    .1513287547049641D+01   K =18   INIAPR=1   IMAX =28   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1513287547049641D+01

YX =    .1776152798388182D+17        DYX =    .7104611193552725D+17
YT =    .1776152798388110D+17        DYT =    .7104611193552441D+17
(YT-YX)/YX = -.4019924415556695D-13    (DYT-DYX)/DYX = -.3997403830595594D-13


           Coefficients AU, ADU and AJK  on   6  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .8594192070744370D+16  .3437676828297744D+17  .1375070731319097D+18
     1      .6975832266300016D+16  .2790332906520004D+17  .1116133162608002D+18
     2      .3979442081664001D+16  .1591776832665602D+17  .6367107330662413D+17
     3      .1710760056458815D+16  .6843040225835279D+16  .2737216090334108D+17
     4      .5842644237944612D+15  .2337057695177849D+16  .9348230780711298D+16
     5      .1647170093225317D+15  .6588680372901237D+15  .2635472149160414D+16
     6      .3943494112136971D+14  .1577397644854761D+15  .6309590579418454D+15
     7      .8191691199523721D+13  .3276676479809305D+14  .1310670591923342D+15
     8      .1501367317821462D+13  .6005469271284706D+13  .2402187708511370D+14
     9      .2460473835238289D+12  .9841895340944611D+12  .3936758136363990D+13
    10      .3644900359573869D+11  .1457960143824269D+12  .5831840575249320D+12
    11      .4924844270546087D+10  .1969937708189944D+11  .7879750832771019D+11
    12      .6115341816248140D+09  .2446136726281303D+10  .9784546908642752D+10
    13      .7023667633845988D+08  .2809467051802657D+09  .1123786827754307D+10
    14      .7502869318538503D+07  .3001147733392174D+08  .1200459188241253D+09
    15      .7490356239539328D+06  .2996142737479480D+07  .1198457576901245D+08
    16      .7018168587332315D+05  .2807268900416248D+06  .1122907456478119D+07
    17      .6194600413695525D+04  .2477831926393555D+05  .9911189156150818D+05
    18      .5167894934183059D+03  .2067066408821073D+04  .8271868242263794D+04
    19      .4087083657321656D+02  .1634054155381356D+03
    20      .3069451319833283D+01
****************************************************************

Results at point XK=    .7000000000000000D+01

Y  =    .7896296018268516D+14        DY  =    .3158518407307306D+15
YT =    .7896296018268069D+14        DYT =    .3158518407307227D+15
(YT-Y)/Y  = -.5659298979751140D-13    (DYT-DY)/DY  = -.2473469833807395D-13

4) Решается задача Коши из примера 1; используется разбиение отрезка интегрирования на семь сегментов, при этом порядок частичных сумм ряда Чебышёва, представляющих приближенное решение, или порядок частичных сумм для оценивающего решения изменяется от сегмента к сегменту. Поэтому каждое обращение к подпрограмме DE98D осуществляется с нулевым значением параметра JSTART. Приводится та же информация, что и в примере 1.

      DOUBLE PRECISION X,YX,DYX,EPS,THRESH,H,HMIN,HMAX,XP,YP,DYP,DELTY,
     1                 RO,YX2,DYX2,AU,ADU,AJK,AU2,ADU2,AJK2,XJ0,XJ02,
     2                 AJKP,U,DU,U2,DU2,ZFI,ZFI2,HD4I,RABC,RABC1,RABC2,
     3                 RABC3,ALP,ALPN,CMAR,P,S,ALP2,ALPN2,CMAR2
C
C
C       take into consideraitoin: all arraies dependent on K
C              must be described for maximal value K
C
C
      LOGICAL BUL
C
      PARAMETER (K=18,K2=27,M=1)
      DIMENSION YP(M),DYP(M),DELTY(M),RO(M),YX2(M),DYX2(M),
     1          IU(2),EPS(2),THRESH(2),MEXACT(2),NUMBES(2),
     2          XJ0 (K),  AU (M,K+3),  ADU (M,K+2), AJK (M,K+1), 
     3          XJ02(K2), AU2(M,K2+3), ADU2(M,K2+2),AJK2(M,K2+1),
     4          AJKP(M,K+1), U(M,K), DU(M,K), U2(M,K2),DU2(M,K2),
     5          HD4I(K2+1),RABC(K2),RABC1(K2),RABC2(K2+2),RABC3(K2+2),
     6          ALP (K), ALPN(K),  CMAR (K,K),  P((K+2)*(K-1)/2),
     7          ALP2(K2),ALPN2(K2),CMAR2(K2,K2),S((K+3)*K/2),
     8          ZFI(M,K+4), ZFI2(M,K2+3)
C
      DOUBLE PRECISION WC1,WC2,WC3,HD2,HD4,HOLD,WC30,WC12,WC22,WC32,
     1                 H2D4,H2D8,H2D16,H2D32,H2D96
C
      COMMON/COM70D/WC1,WC2,WC3,LASN,HD2,HD4,HOLD
      COMMON/COM80D/JLAS,JLAS1,JLAS2,JLAS3,H2D4,H2D8,H2D16,H2D32,H2D96
      COMMON/COM75D/LASN0,LASN2,WC30,WC12,WC22,WC32
      COMMON/COM98D/JLASA2,JLAS12,JLAS22,JLAS32,
     1              JLAS0, JLAS10,JLAS20,JLAS30                         25.03.21   
C
      COMMON/STAT76/NACCEP,NREJEC
      COMMON/BLOCK/N,Q
      DOUBLE PRECISION Q,YT,DYT,XK      
C
      EXTERNAL F
C
C
      DATA IU/2,2/, EPS/2*0.5D-11/, THRESH/2*1.D0/, MEXACT/1,1/,
     1     NUMBES/1,1/ 
C
      M0=M
      Q=4.0D0
C
      N=0
      NREJEC=0
      NACCEP=0
C
      X=0.
      YX=EXP(Q)
      DYX=Q*YX
      XK=7.D0
      K0=12
      INIAPR=1
      IMAX=23
      MESTER=1
      H=1.D0
      HMIN=1.D-3
      HMAX=XK
      K20=25
      IMAX2=3
      NATTEM=3
C
      JSTART=0
      BUL=.FALSE.
C
      DO 20 NX=1,7
      PRINT 54,JSTART
      CALL DE98D(F,M0,K0,INIAPR,IMAX,JSTART,YX,DYX,X,IU,EPS,THRESH,     
     1                 H,HMIN,HMAX,BUL,XP,YP,DYP,MEXACT,NUMBES,DELTY,RO,
     2                 MESTER, K20,IMAX2,YX2,DYX2,NATTEM,AU,ADU,AJK,AU2,
     3                 ADU2,AJK2,XJ0,XJ02,AJKP,U,DU,U2,DU2,ZFI,ZFI2,
     4                 HD4I,RABC,RABC1,RABC2,RABC3,
     5                 ALP,ALPN,CMAR,P,S,ALP2,ALPN2,CMAR2,IERR)         
C
C
C
      YT=EXP(Q*(X+1.D0))
      DYT=Q*YT
      PRINT 44,NX
      PRINT 50,IERR,N,NACCEP,NREJEC,BUL
      PRINT 51,XP,X
      PRINT 35, H,K0,INIAPR,IMAX,JSTART,K20,IMAX2
      PRINT 36, H
      PRINT 52, YX,DYX,YT,DYT,(YT-YX)/YX,(DYT-DYX)/DYX
C
      PRINT 110,NX
      KP1=K0+1
      KP2=K0+2
      KP3=K0+3
C
      DO 10 L=1,M
      PRINT 47
      PRINT 150,L
      PRINT 85,(J-1,AU(L,J  ),ADU(L,J  ),AJK(L,J),J=1,KP1)
      PRINT 86, KP1,AU(L,KP2),ADU(L,KP2),
     1          KP2,AU(L,KP3)
   10 CONTINUE
C
      PRINT 55
C
      GO TO (1,2,3,4,5,6,20),NX
C
C  take into consideraitoin: parameter JSTART must be equal to zero    
C
    1 JSTART=0
C
      K0=  16
      IMAX=25
      GO TO 20     
C
C  take into consideraitoin: parameter JSTART must be equal to zero    
C
    2 JSTART=0
C
      K0=  17
      IMAX=24
      GO TO 20
C
C  take into consideraitoin: parameter JSTART must be equal to zero    
C
    3 JSTART=0
C
      K0=  18
      IMAX=25
      GO TO 20
C
C  take into consideraitoin: parameter JSTART must be equal to zero    
C
    4 JSTART=0
C
      K20=26
      GO TO 20
C
C  take into consideraitoin: parameter JSTART must be equal to zero    
C
    5 JSTART=0
C
      K20=27
      GO TO 20
C
C  take into consideraitoin: parameter JSTART must be equal to zero    
C
    6 JSTART=0
C
      K0=  12
      IMAX=23
      K20= 25
C
      H=XK-X
      BUL=.TRUE.
C
   20 CONTINUE
C2
      STOP
   35 FORMAT(1X,'H =',D25.16,3X,'K =',I2,3X,'INIAPR=',I1,3X,'IMAX =',I2,
     1       3X,'JSTART (exit) =',I1/1X,31X,'K2=',I2,13X,' IMAX2=',I2)
   36 FORMAT(1X,'It is advisable to take the next segment H = ',D22.16/)
   44 FORMAT(/1X,'RESULTS AFTER ',I1,'-nd    CALL SUBROUTINE DE98D'//)
   47 FORMAT(1X,64(1H-)/)
   50 FORMAT(1X,'IERR=',I2,3X,'N=',I5,3X,'NACCEP=',I3,3X,'NREJEC=',I3,
     1           3X,'BUL=',L3)
C
   51 FORMAT(1X,'XP=',D25.16,5X,'X=',D25.16)
C
   52 FORMAT( 1X,'YX =',D25.16,8X,'DYX =',D25.16
     1       /1X,'YT =',D25.16,8X,'DYT =',D25.16
     2       /1X,'(YT-YX)/YX =',D23.16,4X,'(DYT-DYX)/DYX =',D23.16/)
C
   54 FORMAT(1X,'JSTART (input) =',I1)
   55 FORMAT(1X,64(1H*)/)
   85 FORMAT(5X,I2,4X,3D23.16)
   86 FORMAT(5X,I2,4X,2D23.16)
  110 FORMAT(/12X,'Coefficients AU, ADU and AJK  on ',I3,'  segment  ')
C      
  150 FORMAT('   Number of',15X,'Chebyshev coefficients'
     *       ' for ',I1,' component'/  
     1        14X,72(1H-)/
     2       ' coefficient',9X,'for Y',18X,'for Y''',18X,'for Y'''''/         
     3        1X,85(1H-)/)
C   
      END
C         
C
      SUBROUTINE F(X,Y,DY,Z,M)
      DOUBLE PRECISION X,Y,DY,Z,Q
      COMMON/BLOCK/N,Q
      Z=Q*DY
      N=N+1
      RETURN
      END

Результаты:


JSTART (input) =0

RESULTS AFTER 1-nd    CALL SUBROUTINE DE98D


IERR= 0   N=  426   NACCEP=  1   NREJEC=  0   BUL=  F
XP=    .0000000000000000D+00     X=    .1000000000000000D+01
H =    .9294500310870136D+00   K =12   INIAPR=1   IMAX =23   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .9294500310870136D+00

YX =    .2980957987048024D+04        DYX =    .1192383194820469D+05
YT =    .2980957987041728D+04        DYT =    .1192383194816691D+05
(YT-YX)/YX = -.2112065014125831D-11    (DYT-DYX)/DYX = -.3168478898038800D-11


           Coefficients AU, ADU and AJK  on   1  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .1839300696373307D+04  .7357202785500214D+04  .2942881114203202D+05
     1      .1283417414305285D+04  .5133669657226690D+04  .2053467862892990D+05
     2      .5558832820690890D+03  .2223533128279073D+04  .8894132513125262D+04
     3      .1716508501683074D+03  .6866034006739789D+03  .2746413602697313D+04
     4      .4093073156482424D+02  .1637229262593850D+03  .6548917050375144D+03
     5      .7927923909198030D+01  .3171169563679106D+02  .1268467825471533D+03
     6      .1291112018856193D+01  .5164448075424373D+01  .2065779230169317D+02
     7      .1812517960582762D+00  .7250071842424272D+00  .2900028736968347D+01
     8      .2234944644679166D-01  .8939778579263880D-01  .3575911429052044D+00
     9      .2456224492633306D-02  .9824897945094148D-02  .3929959160390536D-01
    10      .2434260183902198D-03  .9737040578397459D-03  .3894816881815023D-02
    11      .2196429511203351D-04  .8785720948535675D-04  .3514292903155219D-03
    12      .1818763512157910D-05  .7275072910271300D-05  .2909966445932555D-04
    13      .1391445827456549D-06  .5565609017770620D-06
    14      .9891601049636025D-08
****************************************************************

JSTART (input) =0

RESULTS AFTER 2-nd    CALL SUBROUTINE DE98D


IERR= 0   N=  984   NACCEP=  2   NREJEC=  0   BUL=  F
XP=    .1000000000000000D+01     X=    .1929450031087014D+01
H =    .1200586501114462D+01   K =16   INIAPR=1   IMAX =25   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1200586501114462D+01

YX =    .1227371215853634D+06        DYX =    .4909484863414662D+06
YT =    .1227371215849774D+06        DYT =    .4909484863399095D+06
(YT-YX)/YX = -.3145321921411109D-11    (DYT-DYX)/DYX = -.3170812675419977D-11


           Coefficients AU, ADU and AJK  on   2  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .7916468317759733D+05  .3166587327104146D+06  .1266634930841657D+07
     1      .5329262534466617D+05  .2131705013786647D+06  .8526820055146578D+06
     2      .2182687741913895D+05  .8730750967655571D+05  .3492300387062221D+06
     3      .6325329225235419D+04  .2530131690094162D+05  .1012052676037663D+06
     4      .1410516086068238D+04  .5642064344272942D+04  .2256825737709174D+05
     5      .2550035965682835D+03  .1020014386273133D+04  .4080057545092543D+04
     6      .3871777502804620D+02  .1548711001121835D+03  .6194844004487265D+03
     7      .5063694047015595D+01  .2025477618806269D+02  .8101910475229565D+02
     8      .5813963701693226D+00  .2325585480678708D+01  .9302341922697829D+01
     9      .5947563205167522D-01  .2379025282063283D+00  .9516101128213716D+00
    10      .5485164013769804D-02  .2194065605460805D-01  .8776262421586090D-01
    11      .4604743745398079D-03  .1841897497444412D-02  .7367590006113200D-02
    12      .3547016204965109D-04  .1418806475471880D-03  .5675226214568729D-03
    13      .2524048537367260D-05  .1009619405206186D-04  .4038482618518024D-04
    14      .1668857038717749D-06  .6675435037173291D-06  .2670210732635958D-05
    15      .1030382708316387D-07  .4121626494026600D-07  .1648735703518867D-06
    16      .5966850719436963D-09  .2387116270875637D-08  .9523429711499887D-08
    17      .3253890747458102D-10  .1297693491743727D-09
    18      .1669565440731301D-11
****************************************************************

JSTART (input) =0

RESULTS AFTER 3-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 1552   NACCEP=  3   NREJEC=  0   BUL=  F
XP=    .1929450031087014D+01     X=    .3130036532201475D+01
H =    .1126362819601387D+01   K =17   INIAPR=1   IMAX =24   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1126362819601387D+01

YX =    .1494886789443820D+08        DYX =    .5979547157775281D+08
YT =    .1494886789439088D+08        DYT =    .5979547157756354D+08
(YT-YX)/YX = -.3164992143198924D-11    (DYT-DYX)/DYX = -.3165365946451929D-11


           Coefficients AU, ADU and AJK  on   3  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .8267992423691862D+07  .3307196969476748D+08  .1322878787790697D+09
     1      .6232455676145561D+07  .2492982270458226D+08  .9971929081832884D+08
     2      .3076816552265530D+07  .1230726620906213D+08  .4922906483624831D+08
     3      .1106933192910194D+07  .4427732771640779D+07  .1771093108656285D+08
     4      .3108354458246381D+06  .1243341783298546D+07  .4973367133193947D+07
     5      .7132119648644638D+05  .2852847859457793D+06  .1141139143782936D+07
     6      .1380896576576041D+05  .5523586306303672D+05  .2209434522520147D+06
     7      .2310096896659963D+04  .9240387586636421D+04  .3696155034646319D+05
     8      .3399835128282131D+03  .1359934051310942D+04  .5439736205191602D+04
     9      .4464738502408613D+02  .1785895400951748D+03  .7143581603491168D+03
    10      .5291706078670853D+01  .2116682431389036D+02  .8466729723846940D+02
    11      .5713765583789090D+00  .2285506233099871D+01  .9142024927239369D+01
    12      .5664626705525970D-01  .2265850680426766D+00  .9063402703159582D+00
    13      .5190611666154220D-02  .2076244651848334D-01  .8304978804481777D-01
    14      .4421104589365038D-03  .1768441844172883D-02  .7073771851302979D-02
    15      .3517596829609333D-04  .1407039972555332D-03  .5628175668448421D-03
    16      .2625640988050752D-05  .1050262379480544D-04  .4200876926319097D-04
    17      .1845663535725190D-06  .7382196299630176D-06  .2951266296236099D-05
    18      .1225832173222870D-07  .4897299056843021D-07
    19      .7709125735525809D-09
****************************************************************

JSTART (input) =0

RESULTS AFTER 4-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 2164   NACCEP=  4   NREJEC=  0   BUL=  F
XP=    .3130036532201475D+01     X=    .4256399351802862D+01
H =    .1242070842037200D+01   K =18   INIAPR=1   IMAX =25   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .1242070842037200D+01

YX =    .1353009770155917D+10        DYX =    .5412039080623665D+10
YT =    .1353009770151628D+10        DYT =    .5412039080606512D+10
(YT-YX)/YX = -.3169728331218697D-11    (DYT-DYX)/DYX = -.3169375904230570D-11


           Coefficients AU, ADU and AJK  on   4  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .7772348948943951D+09  .3108939579577579D+10  .1243575831831031D+11
     1      .5714259557758487D+09  .2285703823103394D+10  .9142815292413572D+10
     2      .2699152765337030D+09  .1079661106134811D+10  .4318644424539247D+10
     3      .9215715919177034D+08  .3686286367670816D+09  .1474514547068326D+10
     4      .2446019513004591D+08  .9784078052018365D+08  .3913631220807337D+09
     5      .5292803570720291D+07  .2117121428288113D+08  .8468485713152386D+08
     6      .9650855675988009D+06  .3860342270395191D+07  .1544136908158039D+08
     7      .1519081995094055D+06  .6076327980376146D+06  .2430531192150084D+07
     8      .2102262619049431D+05  .8409050476196800D+05  .3363620190476874D+06
     9      .2594846300999915D+04  .1037938520400314D+05  .4151754081590103D+05
    10      .2889724311631665D+03  .1155889724643028D+04  .4623558898276469D+04
    11      .2931034638622335D+02  .1172413855332744D+03  .4689655423638587D+03
    12      .2729131368774104D+01  .1091652546531651D+02  .4366610201950607D+02
    13      .2348337211146330D+00  .9393348783029950D+00  .3757340160744207D+01
    14      .1878055933698562D-01  .7512224651846877D-01  .3004894284922557D+00
    15      .1402866673095950D-02  .5611477131613684D-02  .2244610012257908D-01
    16      .9830209700086963D-04  .3932170077602066D-03  .1572754364133289D-02
    17      .6486548498112528D-05  .2594634298012396D-04  .1034871929732617D-03
    18      .4043982999041439D-06  .1615571530494036D-05  .6339624633255880D-05
    19      .2390490460359785D-07  .9616656591795409D-07
    20      .1369867211456028D-08
****************************************************************

JSTART (input) =0

RESULTS AFTER 5-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 2781   NACCEP=  5   NREJEC=  0   BUL=  F
XP=    .4256399351802862D+01     X=    .5498470193840062D+01
H =    .1221473107377168D+01   K =18   INIAPR=1   IMAX =25   JSTART (exit) =1
                               K2=26              IMAX2= 3
It is advisable to take the next segment H =  .1221473107377168D+01

YX =    .1945355530562408D+12        DYX =    .7781422122249636D+12
YT =    .1945355530556247D+12        DYT =    .7781422122224987D+12
(YT-YX)/YX = -.3167443740958422D-11    (DYT-DYX)/DYX = -.3167757489038353D-11


           Coefficients AU, ADU and AJK  on   5  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .1054610070206654D+12  .4218440280826618D+12  .1687376112330647D+13
     1      .8049423111592047D+11  .3219769244636820D+12  .1287907697854728D+13
     2      .4065453351558197D+11  .1626181340623280D+12  .6504725362493110D+12
     3      .1503172746532054D+11  .6012690986128216D+11  .2405076394451275D+12
     4      .4348083939808908D+10  .1739233575923560D+11  .6956934303694147D+11
     5      .1029035210900213D+10  .4116140843600835D+10  .1646456337440255D+11
     6      .2056663897260654D+09  .8226655589042438D+09  .3290662235616311D+10
     7      .3553444057987130D+08  .1421377623194702D+09  .5685510492774353D+09
     8      .5403187627926272D+07  .2161275051169515D+08  .8645100204645504D+08
     9      .7332846711510379D+06  .2933138684595097D+07  .1173255473819118D+08
    10      .8983365728125292D+05  .3593346291196261D+06  .1437338516415389D+07
    11      .1002755697403386D+05  .4011022789576746D+05  .1604409115673014D+06
    12      .1027831572261896D+04  .4111326289467624D+04  .1644530510779115D+05
    13      .9738354958813493D+02  .3895341971278362D+03  .1558136756510416D+04
    14      .8577193625598261D+01  .3430877388266841D+02  .1372350967242382D+03
    15      .7057211218063326D+00  .2822886647594625D+01  .1129154253035085D+02
    16      .5447733932432916D-01  .2179105636889022D+00  .8715390886063688D+00
    17      .3960446959575960D-02  .1584086666432404D-01  .6329724506940693D-01
    18      .2720693548075422D-03  .1086865606229109D-02  .4294728045351803D-02
    19      .1768013624871478D-04  .6962972423624721D-04
    20      .1069628725404883D-05
****************************************************************

JSTART (input) =0

RESULTS AFTER 6-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 3403   NACCEP=  6   NREJEC=  0   BUL=  F
XP=    .5498470193840062D+01     X=    .6719943301217230D+01
H =    .1224981073594249D+01   K =18   INIAPR=1   IMAX =25   JSTART (exit) =1
                               K2=27              IMAX2= 3
It is advisable to take the next segment H =  .1224981073594249D+01

YX =    .2575817593992597D+14        DYX =    .1030327037597044D+15
YT =    .2575817593984434D+14        DYT =    .1030327037593774D+15
(YT-YX)/YX = -.3168896749147471D-11    (DYT-DYX)/DYX = -.3174052879003456D-11


           Coefficients AU, ADU and AJK  on   6  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .1410212592895493D+14  .5640850371582002D+14  .2256340148632810D+15
     1      .1069851155926320D+14  .4279404623705303D+14  .1711761849482128D+15
     2      .5343429978427024D+13  .2137371991370822D+14  .8549487965483295D+14
     3      .1949354584533278D+13  .7797418338133129D+13  .3118967335253240D+14
     4      .5557160956734667D+12  .2222864382693863D+13  .8891457530775346D+13
     5      .1295319709458621D+12  .5181278837834449D+12  .2072511535133706D+13
     6      .2548759467927679D+11  .1019503787171055D+12  .4078015148683719D+12
     7      .4334316448143456D+10  .1733726579257266D+11  .6934906317025052D+11
     8      .6485581472707449D+09  .2594232589082024D+10  .1037693035630500D+11
     9      .8660510122133690D+08  .3464204048850234D+09  .1385681619525028D+10
    10      .1043856339718969D+08  .4175425358853891D+08  .1670170143406147D+09
    11      .1146294685569840D+07  .4585178742252185D+07  .1834071496114313D+08
    12      .1155841476326329D+06  .4623365905668981D+06  .1849346349705592D+07
    13      .1077253353473588D+05  .4309013394912242D+05  .1723605265017897D+06
    14      .9332927920246583D+03  .3733171119452675D+04  .1493267999678850D+05
    15      .7553246136857689D+02  .3021298271105998D+03  .1208512243978679D+04
    16      .5734998938211107D+01  .2293987920381993D+02  .9175593783706427D+02
    17      .4100767932413790D+00  .1640254455460054D+01  .6560060031712055D+01
    18      .2770748550440179D-01  .1107053382643813D+00  .4421771839261055D+00
    19      .1770404098845527D-02  .7030650311694545D-02
    20      .1051944517966898D-03
****************************************************************

JSTART (input) =0

RESULTS AFTER 7-nd    CALL SUBROUTINE DE98D


IERR= 0   N= 3829   NACCEP=  7   NREJEC=  0   BUL=  T
XP=    .6719943301217230D+01     X=    .7000000000000000D+01
H =    .2520510289044934D+01   K =12   INIAPR=1   IMAX =23   JSTART (exit) =1
                               K2=25              IMAX2= 3
It is advisable to take the next segment H =  .2520510289044934D+01

YX =    .7896296018293119D+14        DYX =    .3158518407317252D+15
YT =    .7896296018268069D+14        DYT =    .3158518407307227D+15
(YT-YX)/YX = -.3172373470038028D-11    (DYT-DYX)/DYX = -.3173956490731654D-11


           Coefficients AU, ADU and AJK  on   7  segment  
----------------------------------------------------------------

  Number of               Chebyshev coefficients for 1 component
             ------------------------------------------------------------------------
coefficient         for Y                  for Y'                  for Y''
-------------------------------------------------------------------------------------

     0      .9741284365080634D+14  .3896513746032264D+15  .1558605498412905D+16
     1      .2626434532537098D+14  .1050573813014839D+15  .4202295252059357D+15
     2      .3630600842384231D+13  .1452240336953693D+14  .5808961347814763D+14
     3      .3367323943883587D+12  .1346929577553435D+13  .5387718310213760D+13
     4      .2348418509388529D+11  .9393674037554123D+11  .3757469615020809D+12
     5      .1311957084681376D+10  .5247828338724727D+10  .2099131335491418D+11
     6      .6112307852159707D+08  .2444923140861736D+09  .9779692563162088D+09
     7      .2442001758403434D+07  .9768007033852480D+07  .3907202816908455D+08
     8      .8539442709565666D+05  .3415777085626237D+06  .1366310781895638D+07
     9      .2654941447434037D+04  .1061976571001539D+05  .4247907594776154D+05
    10      .7430046812609403D+02  .2972017115393317D+03  .1188764739036560D+04
    11      .1890544134337072D+01  .7562045681527610D+01  .3028428077697754D+02
    12      .4409945485899004D-01  .1763510591887546D+00  .6838560104370117D+00
    13      .9490487368791738D-03  .3669672555850556D-02
    14      .1715062974786023D-04
****************************************************************

Из приведенных результатов четвертого примера наглядно видно, что на всех выбранных подпрограммой DE98D элементарных сегментах коэффициенты Чебышёва, как для решения Y, так и для его производных Y', Y'', быстро убывают.