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

Текст подпрограммы и версий (Си)
de79d_c.zip

Тексты тестовых примеров (Си)
tde79d1_c.zip , tde79d2_c.zip , tde79d3_c.zip , tde79d4_c.zip

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

Назначение

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

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

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

          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.

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

int de79d_c (U_fp f, int *m, double *xn, double *n, double *dyn, double *xk, int *k, int *iniapr,
                     int *imax, double *h, double *y, double *dy, double *rab)

   Параметры

f -	имя функции вычисления значений правой части дифференциального уравнения. Первый оператор функции должен иметь вид: int f (x, y, dy, d2y, m). Здесь: x, y, dy - значения независимой, зависимой переменных и производной решения соответственно. Вычисленное значение правой части должно быть помещено в d2y. В случае системы уравнений, т.е. когда 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 * m + 5 (тип: с удвоенной точностью).

   Версии:  нет

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

de78d_c -

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

Кроме того, используются рабочие функции de70dk_c, de70dh_c, de80dk_c, de80dh_c, de80d0_c, de80di_c, de70df_c, de70dq_c, de71de_c, de70dp_c, de71dt_c, de71dp_c, de71di_c, de71df_c, de71ds_c, de70da_c, de70dc_c.

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

Разбиение промежутка интегрирования на элементарные сегменты (шаги) выполняется для того, чтобы на каждом таком сегменте ряды Чебышёва для решения и его первой и второй производных были быстросходящимися рядами. Другими словами, длина элементарных сегментов, задаваемая параметром 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_c используются внешние структуры, элементы которых нельзя портить: struct { double wc1, wc2, wc3; double hd2, hd4, hold; int lasn; }com70d_;
struct { int jlas, jlas1, jlas2, jlas3; double h2d4, h2d8, h2d16, h2d32, h2d96; }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 при двух способах определения начального приближения коэффициентов Чебышёва второй производной решения.

 int main(void)
{
    /* Builtin functions */
    double sqrt(double), exp(double);
    /* Local variables */
    extern int f_c();
    static double h__;
    static int m;
    static double y[2], h0;
    static int k0;
    static double dy[2], xk, xn, yn[2], yt[2], rab[469], dyn[2], dyt[2];
    int de79d_c( );
    static double dely[2];
    static int imax, ncheb;
    static double deldy[2], reldy[2];
    static int iniapr;
    static double relerr[2];

    m = 2;
    xn = 0.f;
    yn[0] = 1.;
    yn[1] = .5;
    dyn[0] = 0.;
    dyn[1] = 0.;
    h0 = .1;
    xk = sqrt(2.) * 3.;
    yt[0] = exp(xk * xk);
    yt[1] = exp(-xk * xk) * .5;
    dyt[0] = xk / yt[1];
    dyt[1] = -xk / yt[0];
    k0 = 10;
    imax = 15;
    for (iniapr = 1; iniapr <= 2; ++iniapr) {
	h__ = h0;
	ncheb = 0;
	de79d_c((U_fp)f_c, &m, &xn, yn, dyn, &xk, &k0, &iniapr, &imax, &h__, y,
		dy, rab);
	dely[0] = yt[0] - y[0];
	dely[1] = yt[1] - y[1];
	relerr[0] = dely[0] / yt[0];
	relerr[1] = dely[1] / yt[1];
	deldy[0] = dyt[0] - dy[0];
	deldy[1] = dyt[1] - dy[1];
	reldy[0] = deldy[0] / dyt[0];
	reldy[1] = deldy[1] / dyt[1];
    /*  Операторы вывода на печать:  h0, k0, iniapr, imax, y, yt, dely, relerr, dy, dyt, deldy, reldy  */
    . . . . . . . . . . . . . . . . . . . . . . . . . . 
    return 0;
} /* main */
/*..........................................................................*/
int f_c(double *x, double *y, double *dy,
	double *z__, int *m)
{
    /* Parameter adjustments */
    --z__;
    --dy;
    --y;

    /* Function Body */
    z__[1] = 1. / y[2] + *x * *x / (y[1] * y[2] * y[2]);
    z__[2] = -1. / y[1] + *x * *x / (y[1] * y[1] * y[2]);
    return 0;
} /* f_c */
 
Результаты:
---------------------------------------------------------
 h0=1.00000000000000e-001 k0= 10 iniapr= 1 imax= 15 

 y=       6.56599691373181e+007    7.61498987235807e-009
 yt=      6.56599691373307e+007    7.61498987235629e-009
 dely=    1.26585364341736e-005   -1.78671012311454e-021 
 relerr=  1.92789253490170e-013   -2.34630663082114e-013

 dy=      5.57143313153942e+008   -6.46153317289326e-008
 dyt=     5.57143313154072e+008   -6.46153317289202e-008
 deldy=   1.29699707031250e-004    1.24010917433950e-020 
 reldy=   2.32794155415777e-013   -1.91921815017853e-013
------------------------------------------------------
 h0=1.00000000000000e-001 k0= 10 iniapr= 2 imax= 15 

 y=       6.56599691372846e+007    7.61498987236259e-009
 yt=      6.56599691373307e+007    7.61498987235629e-009
 dely=    4.61041927337646e-005   -6.30807935132939e-021 
 relerr=  7.02165921481560e-013   -8.28376591048242e-013

 dy=      5.57143313153613e+008   -6.46153317289653e-008
 dyt=     5.57143313154072e+008   -6.46153317289202e-008
 deldy=   4.58955764770508e-004    4.51309742208932e-020 
 reldy=   8.23766083042961e-013   -6.98456125091651e-013
------------------------------------------------------

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

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

 int main(void)
{
    /* Builtin functions */
    double sqrt(double), exp(double);
    /* Local variables */
    extern /* Subroutine */ int f_c();
    static double h__;
    static int m;
    static double y[2], h0;
    static int k0;
    static double dy[2], xk, xn, yn[2], yt[2], rab[839], dyn[2], dyt[2];
    int de79d_c( );
    static double dely[2];
    static int imax, ncheb;
    static double deldy[2], reldy[2];
    static int iniapr;
    static double relerr[2];

    m = 2;
    xn = 0.f;
    yn[0] = 1.;
    yn[1] = .5;
    dyn[0] = 0.;
    dyn[1] = 0.;
    h0 = .5;
    xk = sqrt(2.) * 3.;
    yt[0] = exp(xk * xk);
    yt[1] = exp(-xk * xk) * .5;
    dyt[0] = xk / yt[1];
    dyt[1] = -xk / yt[0];
    k0 = 15;
    imax = 28;
    for (iniapr = 1; iniapr <= 2; ++iniapr) {
	h__ = h0;
	ncheb = 0;

	de79d_c((U_fp)f_c, &m, &xn, yn, dyn, &xk, &k0, &iniapr, &imax, &h__, y,
		dy, rab);
	dely[0] = yt[0] - y[0];
	dely[1] = yt[1] - y[1];
	relerr[0] = dely[0] / yt[0];
	relerr[1] = dely[1] / yt[1];
	deldy[0] = dyt[0] - dy[0];
	deldy[1] = dyt[1] - dy[1];
	reldy[0] = deldy[0] / dyt[0];
	reldy[1] = deldy[1] / dyt[1];

    /*  Операторы вывода на печать:  h0, k0, iniapr, imax, y, yt, dely, relerr, dy, dyt, deldy, reldy  */
    . . . . . . . . . . . . . . . . . . . . . . . . . . 
    return 0;
} /* main */
/*..........................................................................*/
int f_c(double *x, double *y, double *dy,
	double *z__, int *m)
{
    /* Parameter adjustments */
    --z__;
    --dy;
    --y;

    /* Function Body */
    z__[1] = 1. / y[2] + *x * *x / (y[1] * y[2] * y[2]);
    z__[2] = -1. / y[1] + *x * *x / (y[1] * y[1] * y[2]);
    return 0;
} /* f_c */
 
Результаты:
---------------------------------------------------------
 h0=5.00000000000000e-001 k0= 15 iniapr= 1 imax= 28 

 y= 6.56599691373779e+007 7.61498987236258e-009
 yt=6.56599691373307e+007 7.61498987235629e-009
 dely=  -4.71696257591248e-005  -6.28988137785322e-021 
 relerr= -7.18392444877142e-013 -8.25986834294628e-013

 dy= 5.57143313153783e+008 -6.46153317288589e-008
 dyt= 5.57143313154072e+008 -6.46153317289202e-008
 deldy=   2.89320945739746e-004  -6.12378351085257e-020 
 reldy= 5.19293580141628e-013 9.47729176187410e-013

------------------------------------------------------
 h0=5.00000000000000e-001 k0= 15 iniapr= 2 imax= 28 

 y= 6.56599691371035e+007 7.61498987238218e-009
 yt=6.56599691373307e+007 7.61498987235629e-009
 dely=   2.27294862270355e-004  -2.58940618953600e-020 
 relerr= 3.46169614845474e-012 -3.40040660977894e-012

 dy= 5.57143313152141e+008 -6.46153317291376e-008
 dyt= 5.57143313154072e+008 -6.46153317289202e-008
 deldy=   1.93142890930176e-003   2.17488944097342e-019 
 reldy= 3.46666443570443e-012 -3.36590308024373e-012
------------------------------------------------------

Здесь при втором способе выбора начального приближения (т.е. при 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.

 int main(void)
{
    /* Builtin functions */
    double sqrt(double), exp(double);
    /* Local variables */
    extern int f_c();
    static double h__;
    static int m;
    static double y[2], h0;
    static int k0;
    static double dy[2], xk, xn, yn[2], yt[2], rab[469], dyn[2], dyt[2];
    int de79d_c( );
    static double dely[2];
    static int imax;
    static double deldy[2];
    static int iniapr;
    static double relerr[2];

    m = 2;
    xn = sqrt(2.) * -3.;
    yn[0] = exp(xn * xn);
    yn[1] = exp(-xn * xn) * .5;
    dyn[0] = xn / yn[1];
    dyn[1] = -xn / yn[0];
    yt[0] = 1.;
    yt[1] = .5;
    dyt[0] = 0.;
    dyt[1] = 0.;
    h0 = .1;
    xk = 0.;
    k0 = 10;
    imax = 14;
    for (iniapr = 1; iniapr <= 2; ++iniapr) {
	h__ = h0;
	de79d_c((U_fp)f_c, &m, &xn, yn, dyn, &xk, &k0, &iniapr, &imax, &h__, y,
		dy, rab);
	dely[0] = yt[0] - y[0];
	dely[1] = yt[1] - y[1];
	relerr[0] = dely[0] / yt[0];
	relerr[1] = dely[1] / yt[1];
	deldy[0] = dyt[0] - dy[0];
	deldy[1] = dyt[1] - dy[1];
    /*  Операторы вывода на печать:  h0, k0, iniapr, imax, y, yt, dely, relerr, dy, dyt, deldy  */
    . . . . . . . . . . . . . . . . . . . . . . . . . . 
    return 0;
} /* main */
/*..........................................................................*/
int f_c(double *x, double *y, double *dy,
	double *z__, int *m)
{
    /* Parameter adjustments */
    --z__;
    --dy;
    --y;

    /* Function Body */
    z__[1] = 1. / y[2] + *x * *x / (y[1] * y[2] * y[2]);
    z__[2] = -1. / y[1] + *x * *x / (y[1] * y[1] * y[2]);
    return 0;
} /* f_c */
 
Результаты:
---------------------------------------------------------
 h0=1.00000000000000e-001 k0= 10 iniapr= 1 imax= 14 

 y=      9.99999999998606e-001    5.00000000000424e-001
 yt=     1.00000000000000e+000    5.00000000000000e-001
 dely=   1.39410705202181e-012   -4.24105195406810e-013 
 relerr= 1.39410705202181e-012   -8.48210390813620e-013

 dy=    -1.38562772367123e-013   -6.17492168508704e-014
 dyt=    0.00000000000000e+000    0.00000000000000e+000
 deldy=  1.38562772367123e-013    6.17492168508704e-014 
------------------------------------------------------
 h0=1.00000000000000e-001 k0= 10 iniapr= 2 imax= 14 

 y=      9.99999999998513e-001    5.00000000000457e-001
 yt=     1.00000000000000e+000    5.00000000000000e-001
 dely=   1.48692169688047e-012   -4.57078819238177e-013 
 relerr= 1.48692169688047e-012   -9.14157638476354e-013

 dy=    -1.47125367444545e-013   -6.51527443107369e-014
 dyt=    0.00000000000000e+000    0.00000000000000e+000
 deldy=  1.47125367444545e-013    6.51527443107369e-014 
------------------------------------------------------

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.

 int main(void)
{
    /* Builtin functions */
    double sqrt(double), exp(double);

    /* Local variables */
    extern int f_c();
    static double h__;
    static int m;
    static double y[2], h0;
    static int k0;
    static double dy[2], xk, xn, yn[2], yt[2], rab[469], dyn[2], dyt[2];
    int de79d_c( );
    static double dely[2];
    static int imax;
    static double deldy[2];
    static int iniapr;
    static double relerr[2];

    m = 2;
    xn = sqrt(2.) * 3.;
    yn[0] = exp(xn * xn);
    yn[1] = exp(-xn * xn) * .5;
    dyn[0] = xn / yn[1];
    dyn[1] = -xn / yn[0];
    xk = 0.;
    yt[0] = 1.;
    yt[1] = .5;
    dyt[0] = 0.;
    dyt[1] = 0.;
    h0 = -.1;
    k0 = 10;
    imax = 14;
    for (iniapr = 1; iniapr <= 2; ++iniapr) {
	h__ = h0;

	de79d_c((U_fp)f_c, &m, &xn, yn, dyn, &xk, &k0, &iniapr, &imax, &h__, y,
		dy, rab);
	dely[0] = yt[0] - y[0];
	dely[1] = yt[1] - y[1];
	relerr[0] = dely[0] / yt[0];
	relerr[1] = dely[1] / yt[1];
	deldy[0] = dyt[0] - dy[0];
	deldy[1] = dyt[1] - dy[1];
    /*  Операторы вывода на печать:  h0, k0, iniapr, imax, y, yt, dely, relerr, dy, dyt, deldy  */
    . . . . . . . . . . . . . . . . . . . . . . . . . . 
    return 0;
} /* main */
/*..........................................................................*/
int f_c(double *x, double *y, double *dy,
	double *z__, int *m)
{
    /* Parameter adjustments */
    --z__;
    --dy;
    --y;

    /* Function Body */
    z__[1] = 1. / y[2] + *x * *x / (y[1] * y[2] * y[2]);
    z__[2] = -1. / y[1] + *x * *x / (y[1] * y[1] * y[2]);
    return 0;
} /* f_c */
 
Результаты:
---------------------------------------------------------
 h0=-1.00000000000000e-001 k0= 10 iniapr= 1 imax= 14 

 y=       9.99999999998606e-001    5.00000000000424e-001
 yt=      1.00000000000000e+000    5.00000000000000e-001
 dely=    1.39410705202181e-012   -4.24105195406810e-013 
 relerr=  1.39410705202181e-012   -8.48210390813620e-013

 dy=      1.38562772367123e-013    6.17492168508704e-014
 dyt=     0.00000000000000e+000    0.00000000000000e+000
 deldy=  -1.38562772367123e-013   -6.17492168508704e-014 
------------------------------------------------------
 h0=-1.00000000000000e-001 k0= 10 iniapr= 2 imax= 14 

 y=       9.99999999998513e-001    5.00000000000457e-001
 yt=      1.00000000000000e+000    5.00000000000000e-001
 dely=    1.48692169688047e-012   -4.57078819238177e-013 
 relerr=  1.48692169688047e-012   -9.14157638476354e-013

 dy=      1.47125367444545e-013    6.51527443107369e-014
 dyt=     0.00000000000000e+000    0.00000000000000e+000
 deldy=  -1.47125367444545e-013   -6.51527443107369e-014 
------------------------------------------------------

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