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

Текст подпрограммы и версий ( Фортран ) de79d.zip	Тексты тестовых примеров ( Фортран ) tde79d1.zip , tde79d2.zip , tde79d3.zip , tde79d4.zip
Текст подпрограммы и версий ( Си ) de79d_c.zip	Тексты тестовых примеров ( Си ) tde79d1_c.zip , tde79d2_c.zip , tde79d3_c.zip , tde79d4_c.zip
Текст подпрограммы и версий ( Паскаль ) de79e_p.zip	Тексты тестовых примеров ( Паскаль ) tde79e1_p.zip , tde79e2_p.zip , tde79e3_p.zip , tde79e4_p.zip

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

Назначение

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

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

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

          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 :

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

методом рядов Чебышёва. Решение Y и его производная Y' вычисляются в одной точке XK, которая является концом интервала интегрирования. Предполагается, что правая часть системы имеет непрерывные ограниченные частные производные по переменным X, Y, Y'. Тогда решение системы и его производные Y', Y'' разлагаются на промежутке интегрирования в равномерно сходящиеся ряды по смещенным многочленам Чебышёва первого рода.

Интервал интегрирования разбивается на элементарные сегменты длины H: [x_s, x_s + H], x₀ = XN, s = 0,1,...,[(XK - XN) / H] - 1 (квадратные скобки означают целую часть числа). Если длина интервала интегрирования не является целым кратным H, то последний элементарный сегмент считается нестандартным; в этом случае s принимает значения s = 0,1,...,[(XK - XN)/H]. На каждом элементарном сегменте решение исходной задачи Коши приближенно представляется в виде (K + 2) - й частичной суммы смещенного ряда Чебышёва

                                   _K+2
           Y(x_s + αH)   ≈   ∑   a_i^*[Y] T_i^*(α) ,     0 ≤ α ≤ 1 ,
                                    ⁱ⁼⁰

а его производные Y'и Y'' приближенно представляются соответственно в виде (K + 1) - й и К - й частичных сумм рядов Чебышёва

                                    _K+1
           Y'(x_s + αH)   ≈   ∑   a_i^*[Y'] T_i^*(α) ,     0 ≤ α ≤ 1 ,
                                     ⁱ⁼⁰

                                       _K
           Y''(x_s + αH)   ≈   ∑   a_i^*[Y''] T_i^*(α) ,     0 ≤ α ≤ 1 .
                                      ⁱ⁼⁰

Здесь a_i^*[Y] , a_i^*[Y'] , a_i^*[Y''] - коэффициенты рядов, T_i^*(α) - смещенные многочлены Чебышёва первого рода на [0,1]. Значения H и K задаются пользователем при обращении к подпрограмме. Они выбираются таким образом, чтобы на каждом элементарном сегменте [x_s, x_s + H] ряды Чебышёва для решения задачи Коши и его первой и второй производных являлись быстросходящимися рядами, а их суммы хорошо приближались частичными суммами (K + 2) - го, (K + 1)- го и К - го порядков соответственно. Решение и его производная Y' в конце каждого сегмента вычисляются по формулам

                                                 _K+2                                                              _K+1
        Y(x_s + H)  =  Y(x_s+1)   ≈   ∑   a_i^*[Y] ,       Y'(x_s + H)  =  Y'(x_s+1)   ≈   ∑   a_i^*[Y'] .
                                                 ⁱ⁼⁰                                                                ⁱ⁼⁰

В качестве решения и его первой производной в конце интервала интегрирования [XN, XK] принимаются значения решения и его производной в конце последнего элементарного сегмента.

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

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

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

    SUBROUTINE  DE79D(F, M, XN, YN, DYN, XK, K, INIAPR, IMAX, H, Y, DY, RAB)

   Параметры

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

INIAPR=1 -

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

INIAPR=2 -

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

IMAX -	целая переменная, задающая количество итераций, которое предполагается выполнить в итерационном процессе вычисления коэффициентов Чебышёва для второй производной решения задачи Коши на каждом элементарном сегменте, исходя из некоторого начального приближения, способ определения которого задается параметром INIAPR; IMAX≥1. Для получения максимального порядка точности приближенного решения необходимо выполнить не менее K итераций (см. "Математическое описание" и "Замечания по использованию");
H -	переменная с двойной точностью, содержащая значение длины элементарных сегментов, на которые разбивается интервал интегрирования (диаметр разбиения промежутка интегрирования, или аналог шага интегрирования для разностных методов). Предполагается, что на каждом элементарном сегменте ряды Чебышёва для решения и его производных являются быстросходящимися рядами. Может задаваться с учетом направления интегрирования, т.е. положительным, если XK > XN, отрицательным, если XK < XN , или без такого учета в виде абсолютной величины;
Y, DY -	искомое решение задачи Коши и его первая производная, вычисленные подпрограммой при значении аргумента XK; для системы уравнений (когда M ≠ 1) задаются одномерными массивами длины M; в случае совпадения значений параметров XN и XK значения Y и DY полагаются равными начальным значениям YN и DYN (тип: с двойной точностью);
RAB -	одномерный рабочий массив длины 2 * K² + 10 * K + 7 * M * K + 12 * М + 5 (тип: с удвоенной точностью).

   Версии:  нет

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

DE78D -

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

Кроме того, используются рабочие подпрограммы DE70DK, DE70DH, DE80DK, DE80DH, DE80D0, DE80DI, DE70DF, DE70DQ, DE71DE, DE70DP, DE71DT, DE71DP, DE71DI, DE71DF, DE71DS, DE70DA, DE70DC.

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

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

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

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

Если правая часть дифференциального уравнения не зависит от переменных Y, Y', т.е. дифференциальное уравнение имеет вид Y'' = F(X), то число итераций можно положить равным 1 при любых H и K, удовлетворяющих описанным выше условиям. В этом случае параметр IMAX = 1.

Как следует из вышеописанного, управлять точностью приближенного решения задачи Коши можно с помощью четырех параметров H, K, IMAX, INIAPR, подбирая для каждой конкретной задачи наиболее подходящий набор их значений.

При работе подпрограммы значения параметров M, XN, YN, DYN, XK, K, INIAPR, IMAX сохраняются. Значение параметра H сохраняется, если он задан с учетом направления интегрирования, иначе его знак меняется на противоположный. Если после работы подпрограммы нет необходимости иметь начальные значения решения YN и/или производной DYN, то параметры YN и Y, а в случае производной, параметры DYN и DY при обращении к ней можно совместить.

Так как при интегрировании системы уравнений с помощью подпрограммы DE79D используются общие блоки с именами COM70D и COM80D, то пользователю не рекомендуется использовать для своих целей общие блоки с указанными именами (пользователь не должен портить содержимое их элементов).

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

Вычисления во всех четырех примерах на Фортране проводились с 15-16 значащими цифрами.

 1)         y''₁  =   1/y₂ + x²/(y₁y₂²) ,     y₁(0)  =  1 ,       y'₁(0)  =  0 ,  
             y'₂  =  -1/y₁ + x²/(y₁²y₂) ,     y₂(0)  =  1/2 ,     y'₂(0)  =  0 ,    0 ≤ x ≤ 3√2 .
     
        Точное решение системы:
                 y₁(x)  =  e^x*x  ,     y₂(x)  =  (1/2)*e^-x*x .

Решение и его первая производная вычисляются в точке XK = 3√2. Приводятся подпрограмма вычисления значений правой части и фрагмент вызывающей программы, а также результаты счета, включая абсолютную DELY и относительную RELERR погрешности приближенного решения,точные значения решения YT, абсолютную DELDY и относительную RELDY погрешности первой производной, точные значения первой производной DYT, значения параметров H, K, IMAX, INIAPR при двух способах определения начального приближения коэффициентов Чебышёва второй производной решения.

      SUBROUTINE F(X,Y,DY,Z,M)
      DOUBLE PRECISION X,Y,Z,DY
      DIMENSION Y(2),Z(2),DY(2)
      Z(1)= 1.D0/Y(2)+X*X/(Y(1)*Y(2)*Y(2))
      Z(2)=-1.D0/Y(1)+X*X/(Y(1)*Y(1)*Y(2))
      RETURN
      END

      PARAMETER (K=10)
      DIMENSION YN(2), Y(2), RAB(2*K*K+10*K+7*2*K+12*2+5), DYN(2), DY(2),
     1          YT(2), DYT(2), DELY(2), RELERR(2), DELDY(2), RELDY(2)
      DOUBLE PRECISION YN, Y, RAB, YT, DELY, XN, H, XK, H0, RELERR,
     1          DYN, DY, DYT, DELDY, RELDY
      EXTERNAL F
      M=2
      XN=0.
      YN(1)=1.D0
      YN(2)=0.5D0
      DYN(1)=0.D0
      DYN(2)=0.D0
      H0=0.1D0
      XK=3.D0*DSQRT(2.D0)
      YT(1)=DEXP(XK*XK)
      YT(2)=0.5D0*DEXP(-XK*XK)
      DYT(1)= XK/YT(2)
      DYT(2)=-XK/YT(1)
      K0=K
      IMAX=15
      DO 42 INIAPR=1,2
      H =H0
      CALL DE79D(F, M, XN, YN, DYN, XK, K0, INIAPR, IMAX, H, Y, DY, RAB)
      DELY(1) =YT(1)-Y(1)
      DELY(2) =YT(2)-Y(2)
      RELERR(1)=DELY(1)/YT(1)
      RELERR(2)=DELY(2)/YT(2)
      DELDY(1)=DYT(1)-DY(1)
      DELDY(2)=DYT(2)-DY(2)
      RELDY(1)=DELDY(1)/DYT(1)
      RELDY(2)=DELDY(2)/DYT(2)
      PRINT 35, H0,K0,INIAPR,IMAX
      PRINT 60, Y,YT , DELY, RELERR
      PRINT 61,DY,DYT,DELDY, RELDY
 42 CONTINUE
      
Результаты:

-----------------------------------------------------------------------
H0=    .1000000000000000D+00   K=10   INIAPR=1   IMAX= 15
Y=      .6565996913732499D+08    .7614989872357003D-08
YT=    .6565996913733028D+08    .7614989872356341D-08
DELTA Y  =    .5289912223815918D-05   -.6617444900424221D-21
RELERR   =    .8056525602002446D-13   -.8690024558596761D-13

DY=       .5571433131540191D+09   -.6461533172892591D-07
DYT=     .5571433131540680D+09   -.6461533172892060D-07
DELTA DY =    .4887580871582031D-04    .5307190810140226D-20
RELDY    =       .8772573871366663D-13   -.8213516309728744D-13
-----------------------------------------------------------------------
H0=    .1000000000000000D+00   K=10   INIAPR=2   IMAX= 15
Y=      .6565996913732578D+08    .7614989872356884D-08
YT=    .6565996913733028D+08    .7614989872356341D-08
DELTA Y  =    .4500150680541992D-05   -.5426304818347862D-21
RELERR   =    .6853720371280954D-13   -.7125820138049343D-13

DY=      .5571433131540281D+09   -.6461533172892506D-07
DYT=    .5571433131540680D+09   -.6461533172892060D-07
DELTA DY =    .3993511199951172D-04    .4460157862885925D-20
RELDY    =       .7167834748555688D-13   -.6902630913662311D-13
----------------------------------------------------------------------------

Видно, что при втором способе выбора начального приближения коэффициентов Чебышёва второй производной (т.е. при INIAPR = 2) приближенное решение и его первая производная вычислены с лучшей точностью, чем при INIAPR = 1, т.е. с меньшей абсолютной и относительной погрешностью.

2) Решается задача Коши из примера 1 с другими значениями параметров H, K, IMAX: H = 0,5, K=15, IMAX = 28. Так как значение H здесь увеличено по сравнению с примером 1, то порядок K частичной суммы и число итераций также возрастают. Приводятся результаты счета, включающие те же данные, что и в примере 1.

H0=   .5000000000000000D+00   K=15   INIAPR=1   IMAX= 28
Y=     .6565996913738872D+08    .7614989872357818D-08
YT=   .6565996913733028D+08    .7614989872356341D-08
DELTA Y  =   -.5844235420227051D-04   -.1477344574019707D-20
RELERR   =   -.8900758707339040D-12   -.1940047982706727D-12

DY=      .5571433131541841D+09   -.6461533172882461D-07
DYT=    .5571433131540680D+09   -.6461533172892060D-07
DELTA DY =   -.1161098480224609D-03   -.9599265572555376D-19
RELDY    =      -.2084021207490520D-12    .1485601840260912D-11
-----------------------------------------------------------------------
H0=    .5000000000000000D+00   K=15   INIAPR=2   IMAX= 28
Y=      .6565996913741745D+08    .7614989872327381D-08
YT=    .6565996913733028D+08    .7614989872356341D-08
DELTA Y  =   -.8717179298400879D-04    .2895959324548150D-19
RELERR   =   -.1327624641456741D-11    .3802971997455907D-11

DY=      .5571433131560063D+09   -.6461533172884166D-07
DYT=    .5571433131540680D+09   -.6461533172892060D-07
DELTA DY =   -.1938343048095703D-02   -.7894611766206096D-19
RELDY    =      -.3479074418254194D-11    .1221786154302543D-11
----------------------------------------------------------------------------

Здесь при втором способе выбора начального приближения (т.е. при INIAPR = 2) приближенное решение и его первая производная вычислены с худшей точностью, чем при INIAPR = 1, т.е. с большей абсолютной и относительной погрешностью.

3) Решается система обыкновенных дифференциальных уравнений из примера 1 на отрезке -3√2 ≤ x ≤ 0. Начальные условия задаются в точке t = -3√2:

     y₁(t)  =  e^t*t ,    y₂(t)  =  (1/2)*e^-t*t  ,    y'₁(t)  =  2te^t*t  ,  y'₂(t)  =  -te^-t*t  .

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

      PARAMETER (K=10)
      DIMENSION YN(2), Y(2), RAB(2*K*K+10*K+7*2*K+12*2+5), DYN(2), DY(2),
     1          YT(2), DYT(2), DELY(2), RELERR(2), DELDY(2)
      DOUBLE PRECISION YN, Y, RAB, YT, DELY, XN, H, XK, H0, RELERR,
     1          DYN, DY, DYT, DELDY
      EXTERNAL F
      M=2
      XN=-3.D0*DSQRT(2.D0)
      YN(1)=DEXP(XN*XN)
      YN(2)=0.5D0*DEXP(-XN*XN)
      DYN(1)= XN/YN(2)
      DYN(2)=-XN/YN(1)
      XK=0.D0
      YT(1)=1.D0
      YT(2)=0.5D0
      DYT(1)=0.D0
      DYT(2)=0.D0
      H0=0.1D0
      K0=K
      IMAX=14
      DO 42 INIAPR=1,2
      H =H0
      CALL DE79D(F, M, XN, YN, DYN, XK, K0, INIAPR, IMAX, H, Y, DY, RAB)
      DELY(1) =YT(1)-Y(1)
      DELY(2) =YT(2)-Y(2)
      RELERR(1)=DELY(1)/YT(1)
      RELERR(2)=DELY(2)/YT(2)
      DELDY(1)=DYT(1)-DY(1)
      DELDY(2)=DYT(2)-DY(2)
      PRINT 35, H0,K0,INIAPR,IMAX
      PRINT 60, Y,YT, DELY, RELERR
      PRINT 61,DY,DYT,DELDY
 42 CONTINUE
  
Результаты:

-----------------------------------------------------------------------
H0=    .1000000000000000D+00   K=10   INIAPR=1   IMAX= 14
Y=      .9999999999997565D+00    .5000000000000473D+00
YT=    .1000000000000000D+01    .5000000000000000D+00
DELTA Y  =    .2434719093002968D-12   -.4729550084903167D-13
RELERR   =    .2434719093002968D-12   -.9459100169806334D-13

DY=      -.4957839694341715D-13   -.2043851199395874D-13
DYT=    .0000000000000000D+00    .0000000000000000D+00
DELTA DY =    .4957839694341715D-13    .2043851199395874D-13

-----------------------------------------------------------------------
H0=    .1000000000000000D+00   K=10   INIAPR=2   IMAX= 14
Y=      .9999999999996718D+00    .5000000000000471D+00
YT=    .1000000000000000D+01    .5000000000000000D+00
DELTA Y  =    .3281819260791963D-12   -.4707345624410664D-13
RELERR   =    .3281819260791963D-12   -.9414691248821327D-13

DY=      -.7678579994063739D-13   -.3276545701424993D-13
DYT=    .0000000000000000D+00    .0000000000000000D+00
DELTA DY =    .7678579994063739D-13    .3276545701424993D-13

-----------------------------------------------------------------------

4) Решается система обыкновенных дифференциальных уравнений из примера 1, но начальное условие задается в точке t = 3√2:

     y₁(t)  =  e^t*t ,    y₂(t)  =  (1/2)*e^-t*t  ,  y'₁(t)  =  2te^t*t  ,   y'₂(t)  =  -te^-t*t  ,

а решение задачи Коши и его первая производная вычисляются в точке XK = 0. Интегрирование системы выполняется справа налево с отрицательным шагом H = -0.1. Приводятся фрагмент вызывающей программы и результаты счета, включающие такие же данные, как и в примере 1.

      PARAMETER (K=10)
      DIMENSION YN(2), Y(2), RAB(2*K*K+10*K+7*2*K+12*2+5), DYN(2), DY(2),
     1          YT(2), DYT(2), DELY(2), RELERR(2), DELDY(2)
      DOUBLE PRECISION YN, Y, RAB, YT, DELY, XN, H, XK, H0, RELERR,
     1          DYN, DY, DYT, DELDY
      EXTERNAL F
      M=2
      XN=3.D0*DSQRT(2.D0)
      YN(1)=DEXP(XN*XN)
      YN(2)=0.5D0*DEXP(-XN*XN)
      DYN(1)= XN/YN(2)
      DYN(2)=-XN/YN(1)
      XK=0.D0
      YT(1)=1.D0
      YT(2)=0.5D0
      DYT(1)=0.D0
      DYT(2)=0.D0
      H0=-0.1D0
      K0=K
      IMAX=14
      DO 42 INIAPR=1,2
      H =H0
      CALL DE79D(F, M, XN, YN, DYN, XK, K0, INIAPR, IMAX, H, Y, DY, RAB)
      DELY(1) =YT(1)-Y(1)
      DELY(2) =YT(2)-Y(2)
      RELERR(1)=DELY(1)/YT(1)
      RELERR(2)=DELY(2)/YT(2)
      DELDY(1)=DYT(1)-DY(1)
      DELDY(2)=DYT(2)-DY(2)
      PRINT 35, H0,K0,INIAPR,IMAX
      PRINT 60, Y,YT , DELY, RELERR
      PRINT 61,DY,DYT,DELDY
 42 CONTINUE
    
Результаты:
   
-----------------------------------------------------------------------
H0=    -.1000000000000000D+00   K=10   INIAPR=1   IMAX= 14
Y=      .9999999999997565D+00    .5000000000000473D+00
YT=    .1000000000000000D+01    .5000000000000000D+00
DELTA Y  =    .2434719093002968D-12   -.4729550084903167D-13
RELERR   =    .2434719093002968D-12   -.9459100169806334D-13

DY=      .4957839694341715D-13    .2043851199395874D-13
DYT=    .0000000000000000D+00    .0000000000000000D+00
DELTA DY =   -.4957839694341715D-13   -.2043851199395874D-13

-----------------------------------------------------------------------
H0=    -.1000000000000000D+00   K=10   INIAPR=2   IMAX= 14
Y=      .9999999999996718D+00    .5000000000000471D+00
YT=    .1000000000000000D+01    .5000000000000000D+00
DELTA Y  =    .3281819260791963D-12   -.4707345624410664D-13
RELERR   =    .3281819260791963D-12   -.9414691248821327D-13

DY=      .7678579994063739D-13    .3276545701424993D-13
DYT=    .0000000000000000D+00    .0000000000000000D+00
DELTA DY =   -.7678579994063739D-13   -.3276545701424993D-13
-----------------------------------------------------------------------

Примеры 3 и 4 относятся к тому случаю, когда решение задачи Коши и его первая производная вычисляются при первом способе выбора начального приближения для коэффициентов Чебышёва второй производной (т.е. при INIAPR = 1) с лучшей точностью, чем при INIAPR = 2, т.е. с меньшей абсолютной и относительной погрешностью.