Текст подпрограммы и версий (Си)
de78d_c.zip
Тексты тестовых примеров (Си)
tde78d1_c.zip , tde78d2_c.zip , tde78d3_c.zip , tde78d4_c.zip

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

Назначение

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

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

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

(1)     Y ' = F (X, Y, Y'),
          Y  = ( y1, ... , yM),  Y'  = ( y'1, ... , y'M),   
          F = ( f1 (X, y1, ... , yM, y'1, ... , y'M), ... , fM (X, y1, ... , yM, y'1, ... , y'M)),    XN ≤ X ≤ XK,
с начальными условиями, заданными в точке XN :
(2)       Y(XN) = YN ,     YN = (y10, ... , yM0),     Y'(XN)=DYN,     DYN=(y'10, ... , y'M0),

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

                                                                                                            1
(3)     Y(XN + αΔ) = ∑ ' ai*[Y] Ti* (α),       0 ≤ α ≤ 1,       ai*[Y] = 2/π ∫ Y(XN + αΔ) Ti* (α) / √ (α (1 - α)) dα ,
                                  i=0                                                                       0
 
                                                                                                              1
(4)     Y'(XN + αΔ) = ∑ ' ai*[Y'] Ti* (α),       0 ≤ α ≤ 1,       ai*[Y'] = 2/π ∫ Y'(XN + αΔ) Ti* (α) / √ (α (1 - α)) dα ,
                                  i=0                                                                         0

                                                                                                   1
(5)     Φ(α) = ∑ ' ai*[Φ] Ti* (α),         0 ≤ α ≤ 1,         ai*[Φ] = 2/π ∫ Φ(α) Ti* (α) / √ (α (1 - α)) dα .
                     i=0                                                                           0

Здесь: штрих у знака суммы означает, что слагаемое с индексом 0 берется с дополнительным множителем 1/2; Ti* (α) - смещенный многочлен Чебышёва первого рода на [0, 1]: Ti* (α) = Ti (2α-1) ;  Ti (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: [xs, xs+H], x0 = XN, s = 0, 1, ... , на каждом из которых ряды Чебышёва для решения Y(X) и его производных Y'(X), Y''(X) будут сходиться значительно быстрее. На каждом подобном сегменте решение исходной задачи Коши приближенно представляется в виде (K + 2) - й частичной суммы смещенного ряда Чебышёва

                                   K+2                                                               1
(6)     Y (xs + αH)   ≈   ∑ ' ai*[Y] Ti*(α),   0 ≤ α ≤ 1,   ai*[Y] = 2/π ∫ Y(xs + αH) Ti* (α) / √ (α (1 - α)) dα ,
                                    i=0                                                               0
а его производные - в виде частичных сумм (K + 1)-го и К-го порядков
                                    K+1                                                                  1
(7)     Y' (xs + αH)   ≈   ∑ ' ai*[Y'] Ti*(α),   0 ≤ α ≤ 1,   ai*[Y ' ] = 2/π ∫ Y ' (xs + αH) Ti* (α) / √ (α (1 - α)) dα ,
                                     i=0                                                                  0

                                      K                                                                    1         
(8)     Y'' (xs + αH)   ≈   ∑ ' ai*[Y''] Ti*(α),   0 ≤ α ≤ 1,   ai*[Y''] = 2/π ∫ Y''(xs + αH) Ti* (α) / √ (α (1 - α)) dα .
                                     i=0                                                                   0

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

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

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

Начальное приближение коэффициентов ai*[Φ] ряда Чебышёва для второй производной на сегменте [xs, xs + H] может вычисляться двумя способами. В первом способе начальное приближение определяется только с использованием значения решения Y(X) и его первой производной Y'(X) в узле xs. При этом погрешность начального приближения для всех коэффициентов a0*[Φ], a1*[Φ], ..., aK*[Φ] является величиной O(H2) при H --> 0. Во втором способе начальное приближение определяется через коэффициенты ряда Чебышёва производной Φ(α) на предыдущем элементарном сегменте [xs - 1, xs] = [X - H', X], где H' - длина предыдущего элементарного сегмента. В этом случае погрешности начального приближения для коэффициентов a0*[Φ], a1*[Φ], ... , aK*[Φ] имеют, соответственно, порядки O(H), O(H2), ..., O(HK + 1). Заметим, что длина H текущего элементарного сегмента [X, X + H] может быть больше или меньше длины H' предыдущего элементарного сегмента [X - H', X] или равна ей. Второй способ определения начального приближения в некоторых случаях может привести к более быстрой сходимости итерационного процесса и, тем самым, к меньшему числу выполняемых итераций. Второй способ может быть применен только начиная со второго элементарного сегмента [x0 + H0, x1 + H], H0 - длина начального (первого) элементарного сегмента, x1 = x0 + H0. На начальном элементарном сегменте [x0, x0 + 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 de78d_c (S_fp f, int *m, int *k, int *iniapr, int *imax, int *jstart, double *yx, double *dyx, 
	double *x, double *xj0, double *au, double *adu, 	double *ajk, double *h, double *hd4i, double *rabc, 
	double *rabc1, double *rabc2, double *rabc3, double *alp, double *alpn, double *cmar, double *p,
              double *s, double *xp, double *yp, double *dyp, double *ajkp, double *u, double *du, double *zfi)

Параметры

f - имя функции вычисления значений правой части дифференциального уравнения. Первый оператор функции должен иметь вид:
int  f (x, y, dy, d2y, m).
Здесь: x, y, dy - значения независимой, зависимой переменных и производной решения соответственно. Вычисленное значение правой части должно быть помещено в d2y. В случае системы уравнений, т.е. когда m ≠ 1 , параметры y, dy, d2y представляют одномерные массивы длины m (тип параметров x, y, dy и d2y: с двойной точностью);
m - количество уравнений в системе (тип: целый);
k - порядок частичной суммы смещенного ряда Чебышёва, с помощью которой аппроксимируется вторая производная решения задачи Коши на элементарном сегменте [x, x + h] разбиения интервала интегрирования; при этом само решение задачи Коши приближается на элементарном сегменте [x, x + h] частичной суммой (k + 2) - го порядка, а его первая производная - частичной суммой (k + 1) - го порядка; k≥2. Если обращение к функции de78d_c осуществляется со значением параметра jstart = 1 (см. ниже), то значение параметра k должно совпадать с его значением при предыдущем обращении к функции. Если же при очередном обращении к функции de78d_c значение параметра k необходимо изменить, то это обращение должно выполняться при нулевом значении параметра jstart (jstart = 0) (см. "Математическое описание", "Замечания по использованию" и "Примеры использования"; тип: целый);
iniapr - целый указатель способа выбора начального приближения коэффициентов Чебышёва для второй производной решения на элементарном сегменте [x, x + h]:
iniapr=1 - для первого способа, когда начальное приближение определяется только с использованием значения решения и его первой производной в начале x элементарного сегмента;
iniapr=2 - для второго способа, когда начальное приближение коэффициентов Чебышёва на текущем элементарном сегменте [x, x + h] (начиная со второго) определяется через коэффициенты Чебышёва, вычисленные на предыдущем элементарном сегменте, т.е. путем экстраполяции коэффициентов с предыдущего сегмента на следующий (см. "Математическое описание").
Значение этого параметра может меняться от сегмента к сегменту.
imax - целая переменная, задающая количество итераций, которое предполагается выполнить в итерационном процессе вычисления коэффициентов Чебышёва для второй производной решения задачи Коши на элементарном сегменте [x, x + h], исходя из некоторого начального приближения, способ определения которого задается параметром iniapr; imax ≥ 1. Для получения максимального порядка точности приближенного решения необходимо выполнить не менее k итераций. Значение imax может изменяться от сегмента к сегменту (см. "Математическое описание", "Замечания по использованию" и "Примеры использования");
jstart - целый указатель режима использования функции, имеющий следующие значения:
0 - первое обращение к функции должно быть исполнено с нулевым значением jstart; выполнить первый (начальный) шаг интегрирования для значений независимой переменной x, зависимой переменной yx, производной dyx и шага интегрирования (длины элементарного сегмента) h. При данном значении параметра jstart будет применен исключительно первый способ определения начального приближения для коэффициентов Чебышёва второй производной решения независимо от значения параметра iniapr. Нулевое значение параметра jstart также может означать, что необходимо выполнить очередной шаг интегрирования с измененным значением параметра k;
1 - выполнить следующий (очередной) шаг интегрирования системы дифференциальных уравнений для значений независимой перемннной x, зависимой переменной yx, производной dyx и шага интегрирования (длины элементарного сегмента) h. При данном значении параметра jstart способ определения начального приближения для коэффициентов Чебышёва второй производной решения определяется параметром iniapr. При обращении к функции de78d_c со значением jstart = 1 значение параметра k (см. выше) должно совпадать с его значением при предыдущем обращении к функции.
На выходе из функции de78d_c параметр jstart всегда принимает значение, равное 1 (см. "Примеры использования");
x, yx, dyx - начальное значение аргумента, решения и первой производной (начало элементарного сегмента x, решение yx и производная dyx в нем); в результате работы функции в x получается новое значение аргумента, равное x + h (конец элементарного сегмента), а в yx и dyx - соответствующие значения решения и производной. В случае системы уравнений, т.е. когда m ≠ 1, yx и dyx задаются одномерными массивами длины m (тип параметров x, yx, dyx: с двойной точностью);
xj0 - одномерный рабочий массив длины k (тип: с удвоенной точностью);
au - двумерный массив размера m * (k + 3). На выходе из функции содержит коэффициенты Чебышёва ai*[Y] для решения Y(x + αh), 0 ≤ α ≤ 1, на элементарном сегменте [x, x + h]. При этом переменная с индексом au(n, i + 1) представляет i-й коэффициент Чебышёва n-й компоненты решения yN(x) (i = 0, 1, ... , k + 2) (тип: с удвоенной точностью);
adu - двумерный массив размера m * (k + 2). На выходе из функции содержит коэффициенты Чебышёва ai*[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). На выходе из функции содержит коэффициенты Чебышёва ai*[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] во время предыдущего обращения к функции de78d_c, используются в функции при вычислении коэффициентов Чебышёва второй производной решения на текущем элементарном сегменте [x, x + h]. Заметим, что длина h текущего элементарного сегмента может быть больше или меньше длины h' предыдущего сегмента [x - h', x] или равна ей. Независимо от значения параметра iniapr при повторных обращениях к функции значения массива ajk при входе в функцию обязательно запоминаются в массиве ajkp и могут быть снова доступны на выходе из функции (тип: с удвоенной точностью);
h - длина текущего элементарного сегмента [x, x + h]. Значение h может быть переменной величиной и изменяться от сегмента к сегменту (тип: с двойной точностью);
hd4i - одномерный рабочий массив длины k + 1 (тип: с двойной точностью);
rabc, rabc1, -
alp,   
alpn   
одномерные рабочие массивы длины k (тип: с двойной точностью);
rabc2, rabc3 - одномерные рабочие массивы длины k + 2 (тип: с двойной точностью);
cmar - двумерный рабочий массив размера k * k (тип: с двойной точностью);
p - одномерный рабочий массив размера (k + 2) * (k - 1) / 2 (тип: с двойной точностью);
s - одномерный рабочий массив размера (k + 3) * k / 2 (тип: с двойной точностью);
xp, yp, dyp - переменная xp и одномерные массивы yp, dyp длины m. На выходе из функции содержат начальные значения аргумента, решения и производной, т.е. те значения, которые имели параметры x, yx, dyx на входе в функцию. Таким образом, на выходе из функции de78d_c границы элементарного сегмента, к которому относятся вычисленные функцией коэффициенты Чебышёва решения и коэффициенты Чебышёва его производных и содержащиеся в массивах au, adu и ajk, показываются параметрами xp и x, а именно: параметр xp содержит начало сегмента, а параметр x содержит конец этого сегмента (тип: с двойной точностью);
ajkp - двумерный массив размера m * (k + 1). Если обращение к функции было выполнено со значением jstart = 1 (т.е. при повторных обращениях к функции), то на выходе из функции содержит коэффициенты Чебышёва второй производной решения, относящиеся к предыдущему элементарному сегменту [x - h', x]; тем самым значение этого параметра определено только при повторных обращениях к функции (тип: с двойной точностью);
u, du - двумерные рабочие массивы размера m * k (тип: с двойной точностью);
zfi - двумерный рабочий массив размера m * (k + 3) (тип: с двойной точностью).

Версии: нет

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

 

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, подбирается таким образом, чтобы убывание коэффициентов этих рядов Чебышёва на элементарном сегменте [x, x + h] происходило достаточно быстро, вследствие чего можно было бы считать частичные суммы этих рядов близкими к многочленам наилучшего равномерного приближения на элементарном сегменте [x, x + h] для решения и его производных. Порядок этих частичных сумм задается параметром k.

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

Если длина h элементарного сегмента [x, x + h] подобрана достаточно малой (или, вернее сказать, выбрана довольно удачной), то хорошая точность приближенного решения может быть получена и с существенно меньшим числом итераций при любом способе выбора начального приближения. Вообще, число итераций зависит от k и h. c увеличением h или k число итераций может также возрастать.

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

Как следует из вышеописанного, управлять точностью приближенного решения задачи Коши можно с помощью четырех параметров h, k, imax, iniapr, подбирая для любой конкретной задачи и для каждого элементарного сегмента наиболее подходящий набор их значений.

При работе функции значения параметров m, k, iniapr, imax, h сохраняются. При многократном использовании функции de78d_c для вычисления коэффициентов Чебышёва решения задачи Коши (1), (2) и его производных на последовательности элементарных сегментов, образующей промежуток интегрирования [xn, xk] системы дифференциальных уравнений, значения параметров m, k, yx, dyx, x, au, adu, ajk, hd4i, rabc, rabc1, rabc2, rabc3, alp, alpn, cmar, p, s не должны изменяться в вызывающей функции между последовательными обращениями к функции.

При интегрировании системы уравнений с помощью функции de78d_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_;

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

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

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

(9)    y''1 = -2qy'2 - ((1 -e3-y1 + y'2/(2q)) / (x+1))2 , 

       y''2 = 2qy'1 - (y'2 - 2q (y1 - 3))2 ,       q = 1/2 ,      0 ≤ x ≤ xf ,     xf = 1 .

Частное решение системы (9) имеет вид

      y1(x) = 3 + cos q(2x - 1) ,                  y2(x) = 2 + sin q(2x - 1) .

Функции y1(x), y2(x) разлагаются на отрезке [0, 1] в смещенные ряды Чебышёва, коэффициенты которых выражаются через цилиндрические функции:

 
                                                
(10)       y1(x) = 3 + J0 (q) + 2 ∑ (-1)i J2i (q) T2i* (x) ,
                                               i=1

                                   
(11)       y2(x) = 2 + 2 ∑ (-1)i J2i+1 (q)Tt2i+1* (x) .
                                  i=0

Задачи Коши отличаются начальными условиями. Вычисления во всех четырех примерах на Си проводились с 15-16 значащими цифрами.

1) Начальное условие задачи Коши для системы (9)

 
xn = 0,         y1(0) = 3 + cos q,         y2(0) = 2 - sin q,         y'1(0) = 2q sin q,       y'2(0) = 2q cos q.

Рассматривается один элементарный сегмент [0, 1]. Выполняется одно обращение к функции de78d_c из начальной точки x = 0 с шагом h =  Вычисляются решение yx и производная dyx в конце данного сегмента, т.е. при x = 1, коэффициенты Чебышёва au решения на данном сегменте, коэффициенты Чебышёва первой производной решения adu на этом же сегменте и коэффициенты Чебышёва ajk его второй производной. Приводятся функция вычисления правой части системы (9), фрагмент вызывающей функции, результаты счета, включая точные значения решения yt и пpoизвoднoй dyt в тoчке x = 1, абсолютную погрешность dely приближенного решения yx и абсолютную погрешность deldy пpиближeннoгo знaчeния пpoизвoднoй dyx, вычиcлeнныx в тoчке x = 1. Кроме вышеперечисленного приводятся: значения параметров xp, x на выходе из функции, которые представляют элементарный сегмент [0, 1]; значения параметров h, k, iniapr, imax, при которых были вычислены приближенное решение yx и коэффициенты Чебышёва. Даются также точные значения коэффициентов Чебышёва на данном сегменте для компонент решения y1(x), y2(x) в (10), (11) и абсолютные погрешности их приближенных значений, вычисленных функцией de78d_c.

 int main(void)
{
    /* System generated locals */
    int i__1, i__2, i__3;
    /* Builtin functions */
    double cos(double), sin(double);
    /* Local variables */
    extern int f_c();
    static double h__;
    static int j, l, m, k;
    static double p[65], q, s[77], u[22] /* was [2][11] */, x, x1,
            au[28] /* was [2][14] */, du[22] /* was [2][11] */, xk;
    static int nx;
    static double xp, yp[2], yt[2], yx[2];
    static int kp1, kp2, kp3;
    static double xj0[11], ajk[24] /* was [2][12] */, adu[26] /* was [2][13] */,
                  alp[11], zfi[28] /* was [2][14] */, dyp[2],
           	  dyt[2], dyx[2], hd4i[12], rabc[11];
    extern int de78d_c( );
    static double cmar[121] /* was [11][11] */, ajkp[24] /* was [2][12] */,
                  alpn[11], dely[2];
    static int imax;
    static double rabc1[11], rabc2[13], rabc3[13], deldy[2];
    static int iniapr, jstart;

    m = 2;
    q = .5;
    x = 0.f;
    yx[0] = cos(q) + 3.;
    yx[1] = -sin(q) + 2.;
    dyx[0] = q * 2. * sin(q);
    dyx[1] = q * 2. * cos(q);
    h__ = 1.;
    xk = x + h__;
    x1 = xk * 2. - 1.;
    yt[0] = cos(q * x1) + 3.;
    yt[1] = sin(q * x1) + 2.;
    dyt[0] = q * -2. * sin(q * x1);
    dyt[1] = q * 2. * cos(q * x1);
    imax = 16;
    jstart = 0;
    iniapr = 1;
    k = 11;

    de78d_c((U_fp)f_c, &m, &k, &iniapr, &imax, &jstart, yx, dyx, &x, xj0,
	    au, adu, ajk, &h__, hd4i, rabc, rabc1, rabc2, rabc3, alp, alpn,
	    cmar, p, s, &xp, yp, dyp, ajkp, u, du, zfi);

    dely[0] = yt[0] - yx[0];
    dely[1] = yt[1] - yx[1];
    deldy[0] = dyt[0] - dyx[0];
    deldy[1] = dyt[1] - dyx[1];

    /*  Операторы вывода на печать:  h0, k0, iniapr, imax, jstart, y, yt, dely, dy, dyt, deldy  */
    . . . . . . . . . . . . . . . . . . . . . . . . . . 
    return 0;
} /* main */
/*..........................................................................*/
 int f_c(double *x, double *y, double *dy,
	double *z__, int *m)
{
    /* System generated locals */
    double d__1;
    /* Builtin functions */
    double exp(double);
    /* Local variables */
    static real q;
    /* Parameter adjustments */
    --z__;
    --dy;
    --y;

    /* Function Body */
    q = .5f;
/* Computing 2nd power */
    d__1 = (1. - exp(3. - y[1] + dy[2] * .5 / q)) / (*x + 1.);
    z__[1] = q * -2. * dy[2] - d__1 * d__1;
/* Computing 2nd power */
    d__1 = dy[2] - q * 2. * (y[1] - 3.);
    z__[2] = q * 2. * dy[1] - d__1 * d__1;
    return 0;
} /* f_c */
 
Результаты:
----------------------------------------------------------------
 xp= 0.00000000000000e+000  x= 1.00000000000000e+000 
 h=1.00000000000000e+000 k = 11 iniapr= 1 imax= 16 jstart=  1

 yx= 3.87758256189037e+000 2.47942553860420e+000
 yt= 3.87758256189037e+000 2.47942553860420e+000
 dely=   0.00000000000000e+000 0.00000000000000e+000

 dyx= -4.79425538604203e-001 8.77582561890373e-001
 dyt= -4.79425538604203e-001 8.77582561890373e-001
 deldy=   1.11022302462516e-016 0.00000000000000e+000 

 ***********************************************************
  COEFFICIENTS AU, ADU AND AJK ON   1 SEGMENT 

                  Chebyshev coefficients  for 1 component 

              for Y                   for  Y'                for Y'' 
   0   7.87693961448163e+000  -2.91650384398601e-017  -1.87693961448163e+000
   1  -5.48877349820787e-018  -4.84536915349748e-001  -2.16305109674436e-016
   2  -6.12080469173653e-002  -7.24043763312976e-018   6.12080469173654e-002
   3  -3.51977482623631e-019   5.12745998917449e-003  -1.58381608609398e-016
   4   3.21472952728575e-004  -3.01670784164619e-018  -3.21472952728538e-004
   5   9.55375519816019e-019  -1.61072544827030e-005  -1.10114283143059e-016
   6  -6.72136925723497e-007  -2.21242182379666e-017   6.72136925521814e-007
   7  -1.86768243938183e-019   2.40317346609178e-008   4.20866954568139e-016
   8   7.51644631713979e-010  -1.68947074076975e-017  -7.51644983884234e-010
   9  -1.04586334371041e-018  -2.08935539295400e-011   9.61497591614457e-016
  10  -5.22635983226684e-013   2.07563729658772e-017   5.22957579204928e-013
  11   4.09594332094098e-019   1.18853995273847e-014   1.31242672979370e-016
  12   2.47612490153849e-016   2.73422235373688e-018
  13   5.25811991103246e-020 
 ------------------------------------------------------------

                  Chebyshev coefficients  for 2 component 

              for Y                   for  Y'                for Y'' 
   0   4.00000000000000e+000   1.87693961448163e+000  -1.14088561731576e-016
   1   4.84536915349748e-001  -2.17213128993893e-017  -4.84536915349748e-001
   2  -2.81571398187064e-018  -6.12080469173653e-002  -2.72033101340191e-017
   3  -5.12745998917449e-003   8.04398955575834e-019   5.12745998917451e-003
   4   6.68746629073058e-019   3.21472952728573e-004  -3.68560976009291e-017
   5   1.61072544827145e-005  -9.89554710959309e-018  -1.61072544826608e-005
   6  -3.76401274308686e-019  -6.72136925715806e-007   1.61054844590933e-016
   7  -2.40317346552878e-008  -8.61916526184626e-019   2.40317345184844e-008
   8  -2.95103169674683e-019   7.51644632253964e-010   1.85188507324102e-016
   9   2.08935352538740e-011   8.58138490340523e-018  -2.08937136424961e-011
  10   2.84842207356998e-019  -5.22636885500724e-013  -1.23741349198486e-016
  11  -1.18836800575604e-014  -2.81230339087469e-018   1.17617775328288e-014
  12  -5.85896539765560e-020   2.45037031933934e-016
  13   4.71225061411411e-018 

 ***********************************************************
           Chebyshev    coefficients   for    1  component 

            approximate                 exact               absolute error 
   0   7.87693961448163e+000   7.87693961448163e+000   0.00000000000000e+000
   1  -5.48877349820787e-018   0.00000000000000e+000   5.48877349820787e-018
   2  -6.12080469173653e-002  -6.12080469173653e-002   6.93889390390723e-018
   3  -3.51977482623631e-019   0.00000000000000e+000   3.51977482623631e-019
   4   3.21472952728575e-004   3.21472952728575e-004   4.87890977618477e-019
   5   9.55375519816019e-019   0.00000000000000e+000  -9.55375519816019e-019
   6  -6.72136925723497e-007  -6.72136925723770e-007  -2.72532850779071e-019
   7  -1.86768243938183e-019   0.00000000000000e+000   1.86768243938183e-019
   8   7.51644631713979e-010   7.51644630959522e-010  -7.54456754253744e-019
   9  -1.04586334371041e-018   0.00000000000000e+000   1.04586334371041e-018
  10  -5.22635983226684e-013  -5.22635472164561e-013   5.11062123000897e-019
  11   4.09594332094098e-019   0.00000000000000e+000  -4.09594332094098e-019
  12   2.47612490153849e-016   2.47676511895987e-016   6.40217421380354e-020
  13   5.25811991103246e-020   0.00000000000000e+000  -5.25811991103246e-020
 ------------------------------------------------------------

           Chebyshev    coefficients   for    2  component 

            approximate                 exact               absolute error 
   0   4.00000000000000e+000   4.00000000000000e+000   0.00000000000000e+000
   1   4.84536915349748e-001   4.84536915349748e-001   0.00000000000000e+000
   2  -2.81571398187064e-018   0.00000000000000e+000   2.81571398187064e-018
   3  -5.12745998917449e-003  -5.12745998917449e-003   0.00000000000000e+000
   4   6.68746629073058e-019   0.00000000000000e+000  -6.68746629073058e-019
   5   1.61072544827145e-005   1.61072544827149e-005   4.84502845829460e-019
   6  -3.76401274308686e-019   0.00000000000000e+000   3.76401274308686e-019
   7  -2.40317346552878e-008  -2.40317346555260e-008  -2.38204855358120e-019
   8  -2.95103169674683e-019   0.00000000000000e+000   2.95103169674683e-019
   9   2.08935352538740e-011   2.08935351786580e-011  -7.52160565496782e-020
  10   2.84842207356998e-019   0.00000000000000e+000  -2.84842207356998e-019
  11  -1.18836800575604e-014  -1.18837079244649e-014  -2.78669045142565e-020
  12  -5.85896539765560e-020   0.00000000000000e+000   5.85896539765560e-020
  13   4.71225061411411e-018   4.76464654243101e-018   5.23959283169000e-020
 ------------------------------------------------------------

2) Начальное условие задачи Коши для системы (9)

   xn = - 0,5,      y1(xn) = 3 + cos q (2 * xn - 1),    y2(xn) = 2 + sin q (2 * xn - 1),

                         y'1(xn) = - 2q sin q (2 * xn - 1),      y'2(xn) = 2q cos q (2*xn-1).

Рассматриваются два элементарных сегмента: [- 0,5, 0] и [0, 1]. Выполняются два обращения к функции de78d_c. Первое обращение осуществляется из начальной точки x = - 0,5 с шагом h = 0,5. Второе обращение осуществляется из точки x = 0 с шагом h = 1. Приводятся фрагмент вызывающей функции и результаты счета. В качестве выходных данных, полученных после каждого обращения к функции de78d_c, приводятся: значения параметров xp, x (которые последовательно представляют границы каждого из двух элементарных сегментов); значения параметров h, k, iniapr, imax, при которых было вычислено приближенное значение решения yx в конце каждого элементарного сегмента; приближенные значения решения yx и производной dyx в конце каждого элементарного сегмента; точные значения решения yt и производной dyt в конце каждого элементарного сегмента; абсолютные погрешности dely и deldy приближенных значений yx и dyx. После этого для каждого элементарного сегмента приводятся приближенные значения коэффициентов Чебышёва au для компонент решения y1(x), y2(x), приближенные значения коэффициентов Чебышёва adu их первых производных y'1(x), y'2(x) и приближенные значения коэффициентов Чебышёва ajk их вторых производных y''1(x), y''2(x). Даются также точные значения коэффициентов Чебышёва решения в (10), (11) и абсолютные погрешности приближенных значений au, вычисленных за два обращения к функции de78d_c.

/* Common Block Declarations */
struct {
    double wc1, wc2, wc3, hd2, hd4, hold;
    int lasn;
} com70d_;
struct {
    int jlas, jlas1, jlas2, jlas3;
    double h2d4, h2d8, h2d16, h2d32, h2d96;
} com80d_;
/*..........................................................................*/
 int main(void)
{
    /* System generated locals */
    int i__1, i__2, i__3;

    /* Builtin functions */
    double cos(double), sin(double);

    /* Local variables */
    extern int f_c();
    static double h__;
    static int j, l, m;
    static double p[65], q, s[77], u[22]  /* was [2][11] */, x, x1,
            au[28] /* was [2][14] */, du[22] /* was [2][11] */, xk;
    static int nx;
    static double xp, yp[2], yt[2], yx[2];
    static int kp3, kp2, kp1;
    static double xj0[11], ajk[24] /* was [2][12] */, adu[26] /* was [2][13] */,
                   alp[11], zfi[28] /* was [2][14] */, dyp[2],
	           dyt[2], dyx[2], hd4i[12], rabc[11];
    extern int de78d_c( );
    static double cmar[121] /* was [11][11] */, ajkp[24] /* was [2][12]*/,
                  alpn[11], dely[2];
    static int imax, k;
    static double rabc1[11], rabc2[13], rabc3[13], deldy[2];
    static int iniapr, jstart;
/* __________________________________________________ */
    m = 2;
    q = .5;
    x = -.5;
    x1 = x * 2. - 1.;
    yx[0] = cos(q * x1) + 3.;
    yx[1] = sin(q * x1) + 2.;
    dyx[0] = q * -2. * sin(q * x1);
    dyx[1] = q * 2. * cos(q * x1);
    h__ = .5;
    xk = x + h__;
    x1 = xk * 2. - 1.;
    yt[0] = cos(q * x1) + 3.;
    yt[1] = sin(q * x1) + 2.;
    dyt[0] = q * -2. * sin(q * x1);
    dyt[1] = q * 2. * cos(q * x1);
    imax = 14;
    jstart = 0;
    iniapr = 1;
    k = 11;

    de78d_c((U_fp)f_c, &m, &k, &iniapr, &imax, &jstart, yx, dyx, &x, xj0,
	    au, adu, ajk, &h__, hd4i, rabc, rabc1, rabc2, rabc3, alp, alpn,
	    cmar, p, s, &xp, yp, dyp, ajkp, u, du, zfi);

    dely[0] = yt[0] - yx[0];
    dely[1] = yt[1] - yx[1];
    deldy[0] = dyt[0] - dyx[0];
    deldy[1] = dyt[1] - dyx[1];

    /*  Операторы вывода на печать:  h, k, iniapr, imax, jstart, yx, yt, dely, dyx, dyt, deldy, au, adu, ajk  */
    . . . . . . . . . . . . . . . . . . . . . . . . . . 
L1:
    h__ = 1.;
    xk = x + h__;
    x1 = xk * 2. - 1.;
    yt[0] = cos(q * x1) + 3.;
    yt[1] = sin(q * x1) + 2.;
    dyt[0] = q * -2. * sin(q * x1);
    dyt[1] = q * 2. * cos(q * x1);
    imax = 16;
    iniapr = 2;

    de78d_c((U_fp)f_c, &m, &k, &iniapr, &imax, &jstart, yx, dyx, &x, xj0,
	    au, adu, ajk, &h__, hd4i, rabc, rabc1, rabc2, rabc3, alp, alpn,
	    cmar, p, s, &xp, yp, dyp, ajkp, u, du, zfi);

    dely[0] = yt[0] - yx[0];
    dely[1] = yt[1] - yx[1];
    deldy[0] = dyt[0] - dyx[0];
    deldy[1] = dyt[1] - dyx[1];

    /*  Операторы вывода на печать:  h, k, iniapr, imax, jstart, yx, yt, dely, dyx, dyt, deldy, au, adu, ajk  */
    . . . . . . . . . . . . . . . . . . . . . . . . . . 
    return 0;
} /* main */
/*..........................................................................*/
 int f_c(double *x, double *y, double *dy,
	double *z__, int *m)
{
    /* System generated locals */
    double d__1;
    /* Builtin functions */
    double exp(double);
    /* Local variables */
    static real q;
    /* Parameter adjustments */
    --z__;
    --dy;
    --y;

    /* Function Body */
    q = .5f;
/* Computing 2nd power */
    d__1 = (1. - exp(3. - y[1] + dy[2] * .5 / q)) / (*x + 1.);
    z__[1] = q * -2. * dy[2] - d__1 * d__1;
/* Computing 2nd power */
    d__1 = dy[2] - q * 2. * (y[1] - 3.);
    z__[2] = q * 2. * dy[1] - d__1 * d__1;
    return 0;
} /* f_c */
 
Результаты:
----------------------------------------------------------------
 RESULTS AFTER THE FIRST CALL FUNCTION DE78D_C 

 xp= -5.00000000000000e-001  x= 0.00000000000000e+000 
 h=5.00000000000000e-001 k = 11 iniapr= 1 imax= 14 jstart=  1

 yx= 3.87758256189037e+000 1.52057446139580e+000
 yt= 3.87758256189037e+000 1.52057446139580e+000
 dely=   0.00000000000000e+000 0.00000000000000e+000

 dyx= 4.79425538604203e-001 8.77582561890373e-001
 dyt= 4.79425538604203e-001 8.77582561890373e-001
 deldy=   0.00000000000000e+000 0.00000000000000e+000 

 ***********************************************************
  COEFFICIENTS AU, ADU AND AJK ON   1 SEGMENT 

                  Chebyshev coefficients  for 1 component 

              for Y                   for  Y'                for Y'' 
   0   7.44060162317046e+000   1.34205937233529e+000  -1.44060162317046e+000
   1   1.69081826785891e-001  -1.81496854116472e-001  -1.69081826785892e-001
   2  -1.13732097613172e-002  -1.05952419518445e-002   1.13732097613173e-002
   3  -4.42044443620038e-004   4.74502064602712e-004   4.42044443619909e-004
   4   1.48397891478288e-005   1.38246950364322e-005  -1.48397891478114e-005
   5   3.45797545845705e-007  -3.71188127810890e-007  -3.45797545922951e-007
   6  -7.73596476256090e-009  -7.20679739597085e-009   7.73596462422640e-009
   7  -1.28728735957424e-010   1.38180792032965e-010   1.28729083650152e-010
   8   2.15954368872478e-012   2.01181764486368e-012  -2.15972961967203e-012
   9   2.79466264948910e-014  -3.00040454205374e-014  -2.72456211238004e-014
  10  -3.75130347672437e-016  -3.39462768468969e-016   5.61650606665381e-016
  11  -3.84704381235065e-018   6.38239325756115e-018  -8.85996462827998e-017
  12   6.64832630995953e-020  -9.22912982112498e-019
  13  -8.87416328954325e-021 

                  Chebyshev coefficients  for 2 component 

              for Y                   for  Y'                for Y'' 
   0   2.65794062766471e+000   1.44060162317046e+000   1.34205937233529e+000
   1   1.81496854116472e-001   1.69081826785892e-001  -1.81496854116472e-001
   2   1.05952419518445e-002  -1.13732097613172e-002  -1.05952419518446e-002
   3  -4.74502064602710e-004  -4.42044443620041e-004   4.74502064602861e-004
   4  -1.38246950364337e-005   1.48397891478312e-005   1.38246950363745e-005
   5   3.71188127814586e-007   3.45797545838455e-007  -3.71188127737259e-007
   6   7.20679738700201e-009  -7.73596475224088e-009  -7.20679716365068e-009
   7  -1.38180791100930e-010  -1.28728737641912e-010   1.38180370303550e-010
   8  -2.01182331352307e-012   2.15954941118904e-012   2.01214429639647e-012
   9   2.99990951620585e-014   2.79544235644900e-014  -3.07920125486826e-014
  10   3.49511857091183e-016  -3.85440479175585e-016  -5.74200246811901e-016
  11  -4.38512213378557e-018  -6.52500280468070e-018   4.32257853642815e-017
  12  -6.79687792154239e-020   4.50268597544598e-019
  13   4.32950574562114e-021 

 ------------------------------------------------------------
 RESULTS AFTER THE SECOND CALL FUNCTION DE78D_C 

 xp= 0.00000000000000e+000  x= 1.00000000000000e+000 
 h=1.00000000000000e+000 k = 11 iniapr= 2 imax= 16 jstart=  1

 yx= 3.87758256189037e+000 2.47942553860420e+000
 yt= 3.87758256189037e+000 2.47942553860420e+000
 dely=   0.00000000000000e+000 0.00000000000000e+000

 dyx= -4.79425538604203e-001 8.77582561890373e-001
 dyt= -4.79425538604203e-001 8.77582561890373e-001
 deldy=   1.11022302462516e-016 0.00000000000000e+000 

************************************************************
  COEFFICIENTS AU, ADU AND AJK ON   2 SEGMENT 

                  Chebyshev coefficients  for 1 component 

              for Y                   for  Y'                for Y'' 
   0   7.87693961448163e+000  -4.51570204840213e-017  -1.87693961448163e+000
   1  -9.17506088465858e-018  -4.84536915349748e-001  -2.35617460871834e-016
   2  -6.12080469173653e-002  -8.44745958296714e-018   6.12080469173654e-002
   3  -4.52527352070610e-019   5.12745998917449e-003  -1.68037784208097e-016
   4   3.21472952728575e-004  -3.01713135811982e-018  -3.21472952728509e-004
   5   1.05591127204721e-018  -1.61072544827016e-005  -1.19763682478180e-016
   6  -6.72136925723422e-007  -2.41353567990640e-017   6.72136925521814e-007
   7  -2.80148584510188e-019   2.40317346605729e-008   4.59484880699357e-016
   8   7.51644631669679e-010  -1.62911964327788e-017  -7.51644974228058e-010
   9  -1.03580012172317e-018  -2.08935528568198e-011   9.80803166548277e-016
  10  -5.22635939953139e-013   2.09976079492552e-017   5.22928617454396e-013
  11   4.10504892512396e-019   1.18847413057817e-014   1.40898848578069e-016
  12   2.47598777203786e-016   2.93539267870978e-018
  13   5.64498592059573e-020 

                  Chebyshev coefficients  for 2 component 

              for Y                   for  Y'                for Y'' 
   0   4.00000000000000e+000   1.87693961448163e+000  -2.72033101340191e-017
   1   4.84536915349748e-001   3.62021881656488e-018  -4.84536915349748e-001
   2   4.52527352070610e-019  -6.12080469173653e-002  -4.16841854002786e-017
   3  -5.12745998917449e-003   0.00000000000000e+000   5.12745998917450e-003
   4   6.33559468722535e-019   3.21472952728574e-004  -4.16841854002786e-017
   5   1.61072544827145e-005  -1.01369514995606e-017  -1.61072544826802e-005
   6  -4.00824057621185e-019  -6.72136925716409e-007   1.61054844590933e-016
   7  -2.40317346552771e-008  -5.17174116652126e-019   2.40317345136563e-008
   8  -2.71762706239231e-019   7.51644631348803e-010   1.75535719857192e-016
   9   2.08935352220250e-011   8.17923248300327e-018  -2.08936895054453e-011
  10   2.72047090217653e-019  -5.22636644096334e-013  -1.18916649530926e-016
  11  -1.18836814275720e-014  -2.70265112570286e-018   1.17762584080951e-014
  12  -5.63052317854762e-020   2.45338716835314e-016
  13   4.71805224683296e-018 

 **************************************************************************
           Chebyshev    coefficients   for    1  component 

            approximate                 exact               absolute error 
   0   7.87693961448163e+000   7.87693961448163e+000   0.00000000000000e+000
   1  -9.17506088465858e-018   0.00000000000000e+000   9.17506088465858e-018
   2  -6.12080469173653e-002  -6.12080469173653e-002   6.93889390390723e-018
   3  -4.52527352070610e-019   0.00000000000000e+000   4.52527352070610e-019
   4   3.21472952728575e-004   3.21472952728575e-004   4.33680868994202e-019
   5   1.05591127204721e-018   0.00000000000000e+000  -1.05591127204721e-018
   6  -6.72136925723422e-007  -6.72136925723770e-007  -3.47389387492670e-019
   7  -2.80148584510188e-019   0.00000000000000e+000   2.80148584510188e-019
   8   7.51644631669679e-010   7.51644630959522e-010  -7.10157303343619e-019
   9  -1.03580012172317e-018   0.00000000000000e+000   1.03580012172317e-018
  10  -5.22635939953139e-013  -5.22635472164561e-013   4.67788578705549e-019
  11   4.10504892512396e-019   0.00000000000000e+000  -4.10504892512396e-019
  12   2.47598777203786e-016   2.47676511895987e-016   7.77346922007797e-020
  13   5.64498592059573e-020   0.00000000000000e+000  -5.64498592059573e-020

 -------------------------------------------------------------
           Chebyshev    coefficients   for    2  component 

            approximate                 exact               absolute error 
   0   4.00000000000000e+000   4.00000000000000e+000   0.00000000000000e+000
   1   4.84536915349748e-001   4.84536915349748e-001   0.00000000000000e+000
   2   4.52527352070610e-019   0.00000000000000e+000  -4.52527352070610e-019
   3  -5.12745998917449e-003  -5.12745998917449e-003   0.00000000000000e+000
   4   6.33559468722535e-019   0.00000000000000e+000  -6.33559468722535e-019
   5   1.61072544827145e-005   1.61072544827149e-005   4.23516473627150e-019
   6  -4.00824057621185e-019   0.00000000000000e+000   4.00824057621185e-019
   7  -2.40317346552771e-008  -2.40317346555260e-008  -2.48987981823362e-019
   8  -2.71762706239231e-019   0.00000000000000e+000   2.71762706239231e-019
   9   2.08935352220250e-011   2.08935351786580e-011  -4.33670309516695e-020
  10   2.72047090217653e-019   0.00000000000000e+000  -2.72047090217653e-019
  11  -1.18836814275720e-014  -1.18837079244649e-014  -2.64968928898489e-020
  12  -5.63052317854762e-020   0.00000000000000e+000   5.63052317854762e-020
  13   4.71805224683296e-018   4.76464654243101e-018   4.65942955980465e-020
 ------------------------------------------------------------

3) Начальное условие задачи Коши для системы (9)

xn = - 4/9,    y1(xn) = 3 + cos q (2 * xn - 1),    y2(xn) = 2 + sin q (2 * xn - 1), 

                   y'1(xn) = - 2q sin q (2 * xn - 1),      y'2(xn) = 2q cos q (2*xn-1).

Рассматриваются три элементарных сегмента: [- 4/9, - 1/3], [- 1/3, 0] и [0, 1]. Выполняются три обращения к функции de78d_c. Первое обращение осуществляется из начальной точки x = - 4/9 с шагом h = 1/9. Второе обращение осуществляется из точки x = - 1/3 с шагом h = 1/3. Третье обращение осуществляется из точки x = 0 с шагом h = 1. Приводятся фрагмент вызывающей функции и результаты счета, аналогичные результатам из примера 2. А именно, после каждого обращения к функции de78d_c приводятся: значения параметров xp, x (которые последовательно представляют границы каждого из трех элементарных сегментов); значения параметров h, k, iniapr, imax, при которых было вычислено приближенное значение решения yx в конце каждого элементарного сегмента; приближенные значения решения yx и производной dyx в конце каждого элементарного сегмента; точные значения решения yt и производной dyt в конце каждого элементарного сегмента; абсолютные погрешности dely и deldy приближенных значений yx и dyx. После этого для каждого элементарного сегмента приводятся приближенные значения коэффициентов Чебышёва au для компонент решения y1(x) и y2(x), приближенные значения коэффициентов Чебышёва adu их первых производных y'1(x), y'2(x) и приближенные значения коэффициентов Чебышёва ajk их вторых производных y''1(x), y''2(x). Даются также, как и в примерах 1 и 2, точные значения коэффициентов Чебышёва для компонент решения y1(x), y2(x) в (10), (11) и абсолютные погрешности приближенных значений au, вычисленных за три обращения к функции de78d_c.

/* Common Block Declarations */
struct {
    double wc1, wc2, wc3, hd2, hd4, hold;
    int lasn;
} com70d_;
struct {
    int jlas, jlas1, jlas2, jlas3;
    double h2d4, h2d8, h2d16, h2d32, h2d96;
} com80d_;
/*..........................................................................*/
 int main(void)
{
    /* System generated locals */
    int i__1, i__2, i__3;
    /* Builtin functions */
    double cos(double), sin(double);
    /* Local variables */
    extern int f_c();
    static double h__;
    static int j, l, m;
    static double p[65], q, s[77], u[22] /* was [2][11] */, x;
    static int k0;
    static double x1, au[28] /* was [2][14] */, du[22] /* was [2][11] */, xk;
    static int nx;
    static double xp, yp[2], yt[2], yx[2];
    static int kp1, kp2, kp3;
    static double xj0[11], ajk[24] /* was [2][12] */, adu[26] /* was [2][13] */,
            alp[11], zfi[28] /* was [2][14] */, dyp[2],
	    dyt[2], dyx[2], hd4i[12], rabc[11];
    extern int de78d_c( );
    static double cmar[121] /* was [11][11] */, ajkp[24] /* was [2][12] */,
                  alpn[11], dely[2];
    static int imax;
    static double rabc1[11], rabc2[13], rabc3[13], deldy[2];
    static int iniapr, jstarn, jstart;

/*       take into consideraitoin: all arraies depending on K */
/*              must be described for maximal value K */

    m = 2;
    q = .5;
    x = -.44444444444444442;
    x1 = x * 2. - 1.;
    yx[0] = cos(q * x1) + 3.;
    yx[1] = sin(q * x1) + 2.;
    dyx[0] = q * -2. * sin(q * x1);
    dyx[1] = q * 2. * cos(q * x1);
    h__ = .1111111111111111;
    xk = x + h__;
    x1 = xk * 2. - 1.;
    yt[0] = cos(q * x1) + 3.;
    yt[1] = sin(q * x1) + 2.;
    dyt[0] = q * -2. * sin(q * x1);
    dyt[1] = q * 2. * cos(q * x1);
    imax = 12;
    jstarn = 0;
    jstart = jstarn;
    iniapr = 1;
    k0 = 6;

    de78d_c((U_fp)f_c, &m, &k0, &iniapr, &imax, &jstart, yx, dyx, &x, xj0, au,
	    adu, ajk, &h__, hd4i, rabc, rabc1, rabc2, rabc3, alp, alpn, cmar,
	    p, s, &xp, yp, dyp, ajkp, u, du, zfi);

    dely[0] = yt[0] - yx[0];
    dely[1] = yt[1] - yx[1];
    deldy[0] = dyt[0] - dyx[0];
    deldy[1] = dyt[1] - dyx[1];

    /*  Операторы вывода на печать:  h, k0, iniapr, imax, jstart, yx, yt, dely, dyx, dyt, deldy, au, adu, ajk  */
    . . . . . . . . . . . . . . . . . . . . . . . . . . 
L1:
    h__ = .33333333333333331;
    xk = x + h__;
    x1 = xk * 2. - 1.;
    yt[0] = cos(q * x1) + 3.;
    yt[1] = sin(q * x1) + 2.;
    dyt[0] = q * -2. * sin(q * x1);
    dyt[1] = q * 2. * cos(q * x1);
    imax = 14;

/*  take into consideraitoin: parameter JSTART must be equal to zero */

    jstart = jstarn;
    k0 = 8;

    de78d_c((U_fp)f_c, &m, &k0, &iniapr, &imax, &jstart, yx, dyx, &x, xj0, au,
	    adu, ajk, &h__, hd4i, rabc, rabc1, rabc2, rabc3, alp, alpn, cmar,
	    p, s, &xp, yp, dyp, ajkp, u, du, zfi);

    dely[0] = yt[0] - yx[0];
    dely[1] = yt[1] - yx[1];
    deldy[0] = dyt[0] - dyx[0];
    deldy[1] = dyt[1] - dyx[1];

    /*  Операторы вывода на печать:  h, k0, iniapr, imax, jstart, yx, yt, dely, dyx, dyt, deldy, au, adu, ajk  */
    . . . . . . . . . . . . . . . . . . . . . . . . . . 
L2:
    h__ = 1.;
    xk = x + h__;
    x1 = xk * 2. - 1.;
    yt[0] = cos(q * x1) + 3.;
    yt[1] = sin(q * x1) + 2.;
    dyt[0] = q * -2. * sin(q * x1);
    dyt[1] = q * 2. * cos(q * x1);
    imax = 17;

/*  take into consideraitoin: parameter JSTART must be equal to zero */

    jstart = jstarn;
    k0 = 11;

    de78d_c((U_fp)f_c, &m, &k0, &iniapr, &imax, &jstart, yx, dyx, &x, xj0, au,
	    adu, ajk, &h__, hd4i, rabc, rabc1, rabc2, rabc3, alp, alpn, cmar,
	    p, s, &xp, yp, dyp, ajkp, u, du, zfi);

    dely[0] = yt[0] - yx[0];
    dely[1] = yt[1] - yx[1];
    deldy[0] = dyt[0] - dyx[0];
    deldy[1] = dyt[1] - dyx[1];

    /*  Операторы вывода на печать:  h, k0, iniapr, imax, jstart, yx, yt, dely, dyx, dyt, deldy, au, adu, ajk  */
    . . . . . . . . . . . . . . . . . . . . . . . . . . 
    return 0;
} /* main */
/*..........................................................................*/
 int f_c(double *x, double *y, double *dy,
	double *z__, int *m)
{
    /* System generated locals */
    double d__1;
    /* Builtin functions */
    double exp(double);
    /* Local variables */
    static real q;
    /* Parameter adjustments */
    --z__;
    --dy;
    --y;

    /* Function Body */
    q = .5f;
/* Computing 2nd power */
    d__1 = (1. - exp(3. - y[1] + dy[2] * .5 / q)) / (*x + 1.);
    z__[1] = q * -2. * dy[2] - d__1 * d__1;
/* Computing 2nd power */
    d__1 = dy[2] - q * 2. * (y[1] - 3.);
    z__[2] = q * 2. * dy[1] - d__1 * d__1;
    return 0;
} /* f_c */
 
Результаты:
----------------------------------------------------------------
 RESULTS AFTER THE FIRST CALL FUNCTION DE78D_C

 xp= -4.44444444444444e-001  x= -3.33333333333333e-001 
 h=1.11111111111111e-001 k0 =  6 iniapr= 1 imax= 12 jstart(initial)=  0

 yx= 3.67241224408306e+000 1.25982314680396e+000
 yt= 3.67241224408306e+000 1.25982314680396e+000
 dely=   0.00000000000000e+000 4.44089209850063e-016

 dyx= 7.40176853196037e-001 6.72412244083057e-001
 dyt= 7.40176853196037e-001 6.72412244083057e-001
 deldy=  -1.11022302462516e-016 0.00000000000000e+000 

 ------------------------------------------------------------
  COEFFICIENTS AU, ADU AND AJK ON   1 SEGMENT 

                  Chebyshev coefficients  for 1 component 

              for Y                   for  Y'                for Y'' 
   0   7.25957764277103e+000   1.55154596888058e+000  -1.25957764277103e+000
   1   4.31151352007386e-002  -3.50017733622692e-002  -4.31151352007385e-002
   2  -4.86198270661226e-004  -5.98898346005549e-004   4.86198270661245e-004
   3  -5.54571166090718e-006   4.50212533913976e-006   5.54571166101887e-006
   4   3.12659655791373e-008   3.85133724271155e-008  -3.12659658487305e-008
   5   2.13968683455316e-010  -1.73704256012294e-010  -2.13968485758434e-010
   6  -8.04201071953641e-013  -9.90594841474232e-013   8.00233482451708e-013
   7  -3.93093191061203e-015   3.17552969226868e-015
   8   1.10261447648218e-017 

                  Chebyshev coefficients  for 2 component 

              for Y                   for  Y'                for Y'' 
   0   2.44845403111942e+000   1.25957764277103e+000   1.55154596888058e+000
   1   3.50017733622692e-002   4.31151352007386e-002  -3.50017733622692e-002
   2   5.98898346005549e-004  -4.86198270661224e-004  -5.98898346005466e-004
   3  -4.50212533913708e-006  -5.54571166091028e-006   4.50212533898354e-006
   4  -3.85133724278731e-008   3.12659655800956e-008   3.85133728450707e-008
   5   1.73704276487934e-010   2.13968703439012e-010  -1.73704550228285e-010
   6   9.90614106418698e-013  -8.04187732538357e-013  -9.93773951476172e-013
   7  -3.19122116086649e-015  -3.94354742649275e-015
   8  -1.36928730086554e-017 

 ------------------------------------------------------------
 RESULTS AFTER THE SECOND CALL FUNCTION DE78D_C

 xp= -3.33333333333333e-001  x= 0.00000000000000e+000 
 h=3.33333333333333e-001 k0 =  8 iniapr= 1 imax= 14 jstart(initial)=  0

 yx= 3.87758256189037e+000 1.52057446139580e+000
 yt= 3.87758256189037e+000 1.52057446139580e+000
 dely=   0.00000000000000e+000 2.22044604925031e-016

 dyx= 4.79425538604203e-001 8.77582561890373e-001
 dyt= 4.79425538604203e-001 8.77582561890373e-001
 deldy=   0.00000000000000e+000 0.00000000000000e+000 

 ------------------------------------------------------------
  COEFFICIENTS AU, ADU AND AJK ON   2 SEGMENT 

                  Chebyshev coefficients  for 1 component 

              for Y                   for  Y'                for Y'' 
   0   7.56087835592219e+000   1.22816603568965e+000  -1.56087835592219e+000
   1   1.02704194891648e-001  -1.30526940340579e-001  -1.02704194891648e-001
   2  -5.44492816475258e-003  -4.28430301012670e-003   5.44492816475268e-003
   3  -1.19077351393356e-004   1.51335613483295e-004   1.19077351393024e-004
   4   3.15392064583599e-006   2.48164003412726e-006  -3.15392064594377e-006
   5   4.13702452650215e-008  -5.25775168324540e-008  -4.13702450847690e-008
   6  -7.30364065483074e-010  -5.74681774027316e-010   7.30364003474824e-010
   7  -6.84229823311261e-012   8.69588232732198e-012   6.84264519774502e-012
   8   9.05907989362457e-014   7.12775541431772e-014  -9.01120202219906e-014
   9   6.59977353177567e-016  -8.34370557611024e-016
  10  -6.95308798009187e-018 

                  Chebyshev coefficients  for 2 component 

              for Y                   for  Y'                for Y'' 
   0   2.77183396431035e+000   1.56087835592219e+000   1.22816603568965e+000
   1   1.30526940340579e-001   1.02704194891648e-001  -1.30526940340578e-001
   2   4.28430301012671e-003  -5.44492816475259e-003  -4.28430301012678e-003
   3  -1.51335613483290e-004  -1.19077351393362e-004   1.51335613483623e-004
   4  -2.48164003413808e-006   3.15392064584662e-006   2.48164003424765e-006
   5   5.25775168318569e-008   4.13702452661160e-008  -5.25775170145628e-008
   6   5.74681771885390e-010  -7.30364064797186e-010  -5.74681719305882e-010
   7  -8.69588862492119e-012  -6.84230963209043e-012   8.69565083458888e-012
   8  -7.12810310624499e-014   9.05796961936342e-014   7.22897897141930e-014
   9   8.38700890681798e-016   6.69349904761047e-016
  10   5.57791587300872e-018 

 ------------------------------------------------------------
 RESULTS AFTER THE THIRD CALL FUNCTION DE78D_C

 xp= 0.00000000000000e+000  x= 1.00000000000000e+000 
 h=1.00000000000000e+000 k0 = 11 iniapr= 1 imax= 17 jstart(initial)=  0

 yx= 3.87758256189037e+000 2.47942553860420e+000
 yt= 3.87758256189037e+000 2.47942553860420e+000
 dely=   0.00000000000000e+000 0.00000000000000e+000

 dyx= -4.79425538604203e-001 8.77582561890373e-001
 dyt= -4.79425538604203e-001 8.77582561890373e-001
 deldy=   1.11022302462516e-016 0.00000000000000e+000 

 ------------------------------------------------------------
  COEFFICIENTS AU, ADU AND AJK ON   3 SEGMENT 

                  Chebyshev coefficients  for 1 component 

              for Y                   for  Y'                for Y'' 
   0   7.87693961448163e+000  -4.65122731996281e-017  -1.87693961448163e+000
   1  -9.82558218814988e-018  -4.84536915349748e-001  -2.16305109674436e-016
   2  -6.12080469173653e-002  -7.24043763312976e-018   6.12080469173654e-002
   3  -3.51977482623631e-019   5.12745998917449e-003  -1.58381608609398e-016
   4   3.21472952728575e-004  -3.01670784164619e-018  -3.21472952728519e-004
   5   1.03582953258939e-018  -1.61072544827020e-005  -1.10114283143059e-016
   6  -6.72136925723428e-007  -2.37332984934340e-017   6.72136925521814e-007
   7  -2.87335759904896e-019   2.40317346602281e-008   4.59484880699357e-016
   8   7.51644631658897e-010  -1.56878972160969e-017  -7.51644964571883e-010
   9  -1.01904651588727e-018  -2.08935528566316e-011   9.61497591614457e-016
  10  -5.22635945431047e-013   2.09977773558447e-017   5.22938266853731e-013
  11   4.19652848342742e-019   1.18849606103121e-014   1.21586497380671e-016
  12   2.47603346048168e-016   2.53305202876399e-018
  13   4.87125390146920e-020 

                  Chebyshev coefficients  for 2 component 

              for Y                   for  Y'                for Y'' 
   0   4.00000000000000e+000   1.87693961448163e+000  -1.14088561731576e-016
   1   4.84536915349748e-001  -1.20676783995320e-017  -4.84536915349748e-001
   2  -1.40790993049452e-018  -6.12080469173653e-002  -6.58178481334482e-017
   3  -5.12745998917449e-003  -8.04398955575834e-019   5.12745998917448e-003
   4   6.58702230040201e-019   3.21472952728573e-004  -5.61650606665381e-017
   5   1.61072544827145e-005  -1.13436346362190e-017  -1.61072544826898e-005
   6  -4.51102521648622e-019  -6.72136925717616e-007   1.70707632057843e-016
   7  -2.40317346553525e-008  -5.17174116652126e-019   2.40317345329652e-008
   8  -2.80141862026480e-019   7.51644632253964e-010   1.85188507324102e-016
   9   2.08935352471683e-011   8.44736546819522e-018  -2.08936991616209e-011
  10   2.78750414847452e-019  -5.22636644096334e-013  -1.18916649530926e-016
  11  -1.18836768555192e-014  -2.70265112570286e-018   1.17666022324964e-014
  12  -5.63052317854762e-020   2.45137546510341e-016
  13   4.71418358673733e-018 

 **************************************************************************
           Chebyshev    coefficients   for    1  component 

            approximate                 exact               absolute error 
   0   7.87693961448163e+000   7.87693961448163e+000   0.00000000000000e+000
   1  -5.14182880301250e-017   0.00000000000000e+000   5.14182880301250e-017
   2  -6.12080469173653e-002  -6.12080469173653e-002   0.00000000000000e+000
   3   1.25710277284379e-018   0.00000000000000e+000  -1.25710277284379e-018
   4   3.21472952728575e-004   3.21472952728575e-004   5.96311194867027e-019
   5   9.15127337605652e-019   0.00000000000000e+000  -9.15127337605652e-019
   6  -6.72136925723411e-007  -6.72136925723770e-007  -3.58824332280603e-019
   7  -2.01127469139255e-019   0.00000000000000e+000   2.01127469139255e-019
   8   7.51644631724742e-010   7.51644630959522e-010  -7.65220338577013e-019
   9  -1.05256431547091e-018   0.00000000000000e+000   1.05256431547091e-018
  10  -5.22635977735515e-013  -5.22635472164561e-013   5.05570953873736e-019
  11   4.10504892512396e-019   0.00000000000000e+000  -4.10504892512396e-019
  12   2.47607918101008e-016   2.47676511895987e-016   6.85937949783599e-020
  13   5.64498592059573e-020   0.00000000000000e+000  -5.64498592059573e-020

 ------------------------------------------------------------
           Chebyshev    coefficients   for    2  component 

            approximate                 exact               absolute error 
   0   4.00000000000000e+000   4.00000000000000e+000   4.44089209850063e-016
   1   4.84536915349748e-001   4.84536915349748e-001   5.55111512312578e-017
   2  -5.38010622968303e-018   0.00000000000000e+000   5.38010622968303e-018
   3  -5.12745998917449e-003  -5.12745998917449e-003  -8.67361737988404e-019
   4   4.27395178664786e-019   0.00000000000000e+000  -4.27395178664786e-019
   5   1.61072544827145e-005   1.61072544827149e-005   4.26904605416167e-019
   6  -3.87891679587213e-019   0.00000000000000e+000   3.87891679587213e-019
   7  -2.40317346553220e-008  -2.40317346555260e-008  -2.04091926896434e-019
   8  -2.66371274305121e-019   0.00000000000000e+000   2.66371274305121e-019
   9   2.08935352312477e-011   2.08935351786580e-011  -5.25896908848523e-020
  10   2.77186612353165e-019   0.00000000000000e+000  -2.77186612353165e-019
  11  -1.18836784007126e-014  -1.18837079244649e-014  -2.95237522936224e-020
  12  -6.05901276228176e-020   0.00000000000000e+000   6.05901276228176e-020
  13   4.71781062525506e-018   4.76464654243101e-018   4.68359171759487e-020
 ------------------------------------------------------------

Заметим, что во втором и третьем примерах последовательные обращения к функции de78d_c выполнялись при одном и том же значении k (k = 11). При этом изменялись только параметры h, imax, iniapr.

4) Четвертый пример служит иллюстрацией того, как последовательные обращения к функции de78d_c можно выполнять не только с разными значениями параметров h, iniapr, imax, но и с разными значениями параметра k.

Рассматриваются начальные условия и элементарные сегменты из примера 3. Выполняются также, как и в примере 3, три обращения к функции. Первое обращение осуществляется из начальной точки x = - 4/9 с шагом h = 1/9 со значением параметра k = 6. Второе обращение осуществляется из точки x = - 1/3 с шагом h = 1/3 и со значением параметра k = 8. Третье обращение осуществляется из точки x = 0 с шагом h = 1 и со значением k = 11. Дополнительным и существенным отличием является только то, что каждое последующее обращение к функции de78d_c с новым значением параметра k (отличным от значения k при предыдущем обращении) выполняется с нулевым значением параметра jstart. Приводится фрагмент вызывающей функции и результаты счета, полностью аналогичные результатам из примера 3. Обратим внимание на то, что в примере 4 дается значение параметра jstart на входе в функцию de78d_c после каждого обращения к ней, в то время как в первых трех примерах приведено значение параметра jstart на выходе из функции de78d_c, которое всегда полагается равным 1.

 int main(void)
{
    /* System generated locals */
    int i__1, i__2, i__3;
    /* Builtin functions */
    double cos(double), sin(double);
    /* Local variables */
    extern int f_c();
    static double h__;
    static int j, l, m;
    static double p[65], q, s[77], u[22] /* was [2][11] */, x;
    static int k0;
    static double x1, au[28] /* was [2][14] */, du[22] /* was [2][11] */, xk;
    static int nx;
    static double xp, yp[2], yt[2], yx[2];
    static int kp1, kp2, kp3;
    static double xj0[11], ajk[24] /* was [2][12] */, adu[26] /* was [2][13] */,
            alp[11], zfi[28] /* was [2][14] */, dyp[2],
	    dyt[2], dyx[2], hd4i[12], rabc[11];
    extern int de78d_c( );
    static double cmar[121] /* was [11][11] */, ajkp[24] /* was [2][12] */,
                  alpn[11], dely[2];
    static int imax;
    static double rabc1[11], rabc2[13], rabc3[13], deldy[2];
    static int iniapr, jstarn, jstart;

/*       take into consideraitoin: all arraies depending on K */
/*              must be described for maximal value K */

    m = 2;
    q = .5;
    x = -.44444444444444442;
    x1 = x * 2. - 1.;
    yx[0] = cos(q * x1) + 3.;
    yx[1] = sin(q * x1) + 2.;
    dyx[0] = q * -2. * sin(q * x1);
    dyx[1] = q * 2. * cos(q * x1);
    h__ = .1111111111111111;
    xk = x + h__;
    x1 = xk * 2. - 1.;
    yt[0] = cos(q * x1) + 3.;
    yt[1] = sin(q * x1) + 2.;
    dyt[0] = q * -2. * sin(q * x1);
    dyt[1] = q * 2. * cos(q * x1);
    imax = 12;
    jstarn = 0;
    jstart = jstarn;
    iniapr = 1;
    k0 = 6;

    de78d_c((U_fp)f_c, &m, &k0, &iniapr, &imax, &jstart, yx, dyx, &x, xj0, au,
	    adu, ajk, &h__, hd4i, rabc, rabc1, rabc2, rabc3, alp, alpn, cmar,
	    p, s, &xp, yp, dyp, ajkp, u, du, zfi);

    dely[0] = yt[0] - yx[0];
    dely[1] = yt[1] - yx[1];
    deldy[0] = dyt[0] - dyx[0];
    deldy[1] = dyt[1] - dyx[1];

/*  Операторы вывода на печать:  h0, k, iniapr, imax, jstart, y, yt, dely, dy, dyt, deldy, au, adu, ajk  */
    . . . . . . . . . . . . . . . . . . . . . . . . . . 
L1:
    h__ = .33333333333333331;
    xk = x + h__;
    x1 = xk * 2. - 1.;
    yt[0] = cos(q * x1) + 3.;
    yt[1] = sin(q * x1) + 2.;
    dyt[0] = q * -2. * sin(q * x1);
    dyt[1] = q * 2. * cos(q * x1);
    imax = 14;

/*  take into consideraitoin: parameter JSTART must be equal to zero */

    jstart = jstarn;
    k0 = 8;

    de78d_c((U_fp)f_c, &m, &k0, &iniapr, &imax, &jstart, yx, dyx, &x, xj0, au,
	    adu, ajk, &h__, hd4i, rabc, rabc1, rabc2, rabc3, alp, alpn, cmar,
	    p, s, &xp, yp, dyp, ajkp, u, du, zfi);

    dely[0] = yt[0] - yx[0];
    dely[1] = yt[1] - yx[1];
    deldy[0] = dyt[0] - dyx[0];
    deldy[1] = dyt[1] - dyx[1];

/*  Операторы вывода на печать:  h, k0, iniapr, imax, jstart, y, yt, dely, dy, dyt, deldy, au, adu, ajk  */
    . . . . . . . . . . . . . . . . . . . . . . . . . . 
L2:
    h__ = 1.;
    xk = x + h__;
    x1 = xk * 2. - 1.;
    yt[0] = cos(q * x1) + 3.;
    yt[1] = sin(q * x1) + 2.;
    dyt[0] = q * -2. * sin(q * x1);
    dyt[1] = q * 2. * cos(q * x1);
    imax = 17;

/*  take into consideraitoin: parameter JSTART must be equal to zero */

    jstart = jstarn;
    k0 = 11;

    de78d_c((U_fp)f_c, &m, &k0, &iniapr, &imax, &jstart, yx, dyx, &x, xj0, au,
	    adu, ajk, &h__, hd4i, rabc, rabc1, rabc2, rabc3, alp, alpn, cmar,
	    p, s, &xp, yp, dyp, ajkp, u, du, zfi);

    dely[0] = yt[0] - yx[0];
    dely[1] = yt[1] - yx[1];
    deldy[0] = dyt[0] - dyx[0];
    deldy[1] = dyt[1] - dyx[1];

/*  Операторы вывода на печать:  h, k0, iniapr, imax, jstart, y, yt, dely, dy, dyt, deldy, au, adu, ajk  */
    . . . . . . . . . . . . . . . . . . . . . . . . . . 
    return 0;
} /* main */
/*..........................................................................*/
 int f_c(double *x, double *y, double *dy,
	double *z__, int *m)
{
    /* System generated locals */
    double d__1;
    /* Builtin functions */
    double exp(double);
    /* Local variables */
    static real q;
    /* Parameter adjustments */
    --z__;
    --dy;
    --y;

    /* Function Body */
    q = .5f;
/* Computing 2nd power */
    d__1 = (1. - exp(3. - y[1] + dy[2] * .5 / q)) / (*x + 1.);
    z__[1] = q * -2. * dy[2] - d__1 * d__1;
/* Computing 2nd power */
    d__1 = dy[2] - q * 2. * (y[1] - 3.);
    z__[2] = q * 2. * dy[1] - d__1 * d__1;
    return 0;
} /* f_c */
 
Результаты:
----------------------------------------------------------------
 RESULTS AFTER THE FIRST CALL FUNCTION DE78D_C

 xp= -4.44444444444444e-001  x= -3.33333333333333e-001 
 h=1.11111111111111e-001 k0 =  6 iniapr= 1 imax= 12 jstart(initial)=  0

 yx= 3.67241224408306e+000 1.25982314680396e+000
 yt= 3.67241224408306e+000 1.25982314680396e+000
 dely=   0.00000000000000e+000 4.44089209850063e-016

 dyx= 7.40176853196037e-001 6.72412244083057e-001
 dyt= 7.40176853196037e-001 6.72412244083057e-001
 deldy=  -1.11022302462516e-016 0.00000000000000e+000 
 ------------------------------------------------------------
  COEFFICIENTS AU, ADU AND AJK ON   1 SEGMENT 

                  Chebyshev coefficients  for 1 component 

              for Y                   for  Y'                for Y'' 
   0   7.25957764277103e+000   1.55154596888058e+000  -1.25957764277103e+000
   1   4.31151352007386e-002  -3.50017733622692e-002  -4.31151352007385e-002
   2  -4.86198270661226e-004  -5.98898346005549e-004   4.86198270661245e-004
   3  -5.54571166090718e-006   4.50212533913976e-006   5.54571166101887e-006
   4   3.12659655791373e-008   3.85133724271155e-008  -3.12659658487305e-008
   5   2.13968683455316e-010  -1.73704256012294e-010  -2.13968485758434e-010
   6  -8.04201071953641e-013  -9.90594841474232e-013   8.00233482451708e-013
   7  -3.93093191061203e-015   3.17552969226868e-015
   8   1.10261447648218e-017 

                  Chebyshev coefficients  for 2 component 

              for Y                   for  Y'                for Y'' 
   0   2.44845403111942e+000   1.25957764277103e+000   1.55154596888058e+000
   1   3.50017733622692e-002   4.31151352007386e-002  -3.50017733622692e-002
   2   5.98898346005549e-004  -4.86198270661224e-004  -5.98898346005466e-004
   3  -4.50212533913708e-006  -5.54571166091028e-006   4.50212533898354e-006
   4  -3.85133724278731e-008   3.12659655800956e-008   3.85133728450707e-008
   5   1.73704276487934e-010   2.13968703439012e-010  -1.73704550228285e-010
   6   9.90614106418698e-013  -8.04187732538357e-013  -9.93773951476172e-013
   7  -3.19122116086649e-015  -3.94354742649275e-015
   8  -1.36928730086554e-017 

 ------------------------------------------------------------
 RESULTS AFTER THE SECOND CALL FUNCTION DE78D_C

 xp= -3.33333333333333e-001  x= 0.00000000000000e+000 
 h=3.33333333333333e-001 k0 =  8 iniapr= 1 imax= 14 jstart(initial)=  0

 yx= 3.87758256189037e+000 1.52057446139580e+000
 yt= 3.87758256189037e+000 1.52057446139580e+000
 dely=   0.00000000000000e+000 2.22044604925031e-016

 dyx= 4.79425538604203e-001 8.77582561890373e-001
 dyt= 4.79425538604203e-001 8.77582561890373e-001
 deldy=   0.00000000000000e+000 0.00000000000000e+000 

 ------------------------------------------------------------
  COEFFICIENTS AU, ADU AND AJK ON   2 SEGMENT 

                  Chebyshev coefficients  for 1 component 

              for Y                   for  Y'                for Y'' 
   0   7.56087835592219e+000   1.22816603568965e+000  -1.56087835592219e+000
   1   1.02704194891648e-001  -1.30526940340579e-001  -1.02704194891648e-001
   2  -5.44492816475258e-003  -4.28430301012670e-003   5.44492816475268e-003
   3  -1.19077351393356e-004   1.51335613483295e-004   1.19077351393024e-004
   4   3.15392064583599e-006   2.48164003412726e-006  -3.15392064594377e-006
   5   4.13702452650215e-008  -5.25775168324540e-008  -4.13702450847690e-008
   6  -7.30364065483074e-010  -5.74681774027316e-010   7.30364003474824e-010
   7  -6.84229823311261e-012   8.69588232732198e-012   6.84264519774502e-012
   8   9.05907989362457e-014   7.12775541431772e-014  -9.01120202219906e-014
   9   6.59977353177567e-016  -8.34370557611024e-016
  10  -6.95308798009187e-018 

                  Chebyshev coefficients  for 2 component 

              for Y                   for  Y'                for Y'' 
   0   2.77183396431035e+000   1.56087835592219e+000   1.22816603568965e+000
   1   1.30526940340579e-001   1.02704194891648e-001  -1.30526940340578e-001
   2   4.28430301012671e-003  -5.44492816475259e-003  -4.28430301012678e-003
   3  -1.51335613483290e-004  -1.19077351393362e-004   1.51335613483623e-004
   4  -2.48164003413808e-006   3.15392064584662e-006   2.48164003424765e-006
   5   5.25775168318569e-008   4.13702452661160e-008  -5.25775170145628e-008
   6   5.74681771885390e-010  -7.30364064797186e-010  -5.74681719305882e-010
   7  -8.69588862492119e-012  -6.84230963209043e-012   8.69565083458888e-012
   8  -7.12810310624499e-014   9.05796961936342e-014   7.22897897141930e-014
   9   8.38700890681798e-016   6.69349904761047e-016
  10   5.57791587300872e-018 

 ------------------------------------------------------------
 RESULTS AFTER THE THIRD CALL FUNCTION DE78D_C

 xp= 0.00000000000000e+000  x= 1.00000000000000e+000 
 h=1.00000000000000e+000 k0 = 11 iniapr= 1 imax= 17 jstart(initial)=  0

 yx= 3.87758256189037e+000 2.47942553860420e+000
 yt= 3.87758256189037e+000 2.47942553860420e+000
 dely=   0.00000000000000e+000 0.00000000000000e+000

 dyx= -4.79425538604203e-001 8.77582561890373e-001
 dyt= -4.79425538604203e-001 8.77582561890373e-001
 deldy=   1.11022302462516e-016 0.00000000000000e+000 

 ------------------------------------------------------------
  COEFFICIENTS AU, ADU AND AJK ON   3 SEGMENT 

                  Chebyshev coefficients  for 1 component 

              for Y                   for  Y'                for Y'' 
   0   7.87693961448163e+000  -4.65122731996281e-017  -1.87693961448163e+000
   1  -9.82558218814988e-018  -4.84536915349748e-001  -2.16305109674436e-016
   2  -6.12080469173653e-002  -7.24043763312976e-018   6.12080469173654e-002
   3  -3.51977482623631e-019   5.12745998917449e-003  -1.58381608609398e-016
   4   3.21472952728575e-004  -3.01670784164619e-018  -3.21472952728519e-004
   5   1.03582953258939e-018  -1.61072544827020e-005  -1.10114283143059e-016
   6  -6.72136925723428e-007  -2.37332984934340e-017   6.72136925521814e-007
   7  -2.87335759904896e-019   2.40317346602281e-008   4.59484880699357e-016
   8   7.51644631658897e-010  -1.56878972160969e-017  -7.51644964571883e-010
   9  -1.01904651588727e-018  -2.08935528566316e-011   9.61497591614457e-016
  10  -5.22635945431047e-013   2.09977773558447e-017   5.22938266853731e-013
  11   4.19652848342742e-019   1.18849606103121e-014   1.21586497380671e-016
  12   2.47603346048168e-016   2.53305202876399e-018
  13   4.87125390146920e-020 

                  Chebyshev coefficients  for 2 component 

              for Y                   for  Y'                for Y'' 
   0   4.00000000000000e+000   1.87693961448163e+000  -1.14088561731576e-016
   1   4.84536915349748e-001  -1.20676783995320e-017  -4.84536915349748e-001
   2  -1.40790993049452e-018  -6.12080469173653e-002  -6.58178481334482e-017
   3  -5.12745998917449e-003  -8.04398955575834e-019   5.12745998917448e-003
   4   6.58702230040201e-019   3.21472952728573e-004  -5.61650606665381e-017
   5   1.61072544827145e-005  -1.13436346362190e-017  -1.61072544826898e-005
   6  -4.51102521648622e-019  -6.72136925717616e-007   1.70707632057843e-016
   7  -2.40317346553525e-008  -5.17174116652126e-019   2.40317345329652e-008
   8  -2.80141862026480e-019   7.51644632253964e-010   1.85188507324102e-016
   9   2.08935352471683e-011   8.44736546819522e-018  -2.08936991616209e-011
  10   2.78750414847452e-019  -5.22636644096334e-013  -1.18916649530926e-016
  11  -1.18836768555192e-014  -2.70265112570286e-018   1.17666022324964e-014
  12  -5.63052317854762e-020   2.45137546510341e-016
  13   4.71418358673733e-018 

 **************************************************************************
           Chebyshev    coefficients   for    1  component 

            approximate                 exact               absolute error 
   0   7.87693961448163e+000   7.87693961448163e+000   0.00000000000000e+000
   1  -9.82558218814988e-018   0.00000000000000e+000   9.82558218814988e-018
   2  -6.12080469173653e-002  -6.12080469173653e-002   6.93889390390723e-018
   3  -3.51977482623631e-019   0.00000000000000e+000   3.51977482623631e-019
   4   3.21472952728575e-004   3.21472952728575e-004   4.33680868994202e-019
   5   1.03582953258939e-018   0.00000000000000e+000  -1.03582953258939e-018
   6  -6.72136925723428e-007  -6.72136925723770e-007  -3.41671915098703e-019
   7  -2.87335759904896e-019   0.00000000000000e+000   2.87335759904896e-019
   8   7.51644631658897e-010   7.51644630959522e-010  -6.99374487071108e-019
   9  -1.01904651588727e-018   0.00000000000000e+000   1.01904651588727e-018
  10  -5.22635945431047e-013  -5.22635472164561e-013   4.73266486184747e-019
  11   4.19652848342742e-019   0.00000000000000e+000  -4.19652848342742e-019
  12   2.47603346048168e-016   2.47676511895987e-016   7.31658478186350e-020
  13   4.87125390146920e-020   0.00000000000000e+000  -4.87125390146920e-020

 ------------------------------------------------------------
           Chebyshev    coefficients   for    2  component 

            approximate                 exact               absolute error 
   0   4.00000000000000e+000   4.00000000000000e+000   4.44089209850063e-016
   1   4.84536915349748e-001   4.84536915349748e-001   0.00000000000000e+000
   2  -1.40790993049452e-018   0.00000000000000e+000   1.40790993049452e-018
   3  -5.12745998917449e-003  -5.12745998917449e-003   0.00000000000000e+000
   4   6.58702230040201e-019   0.00000000000000e+000  -6.58702230040201e-019
   5   1.61072544827145e-005   1.61072544827149e-005   4.23516473627150e-019
   6  -4.51102521648622e-019   0.00000000000000e+000   4.51102521648622e-019
   7  -2.40317346553525e-008  -2.40317346555260e-008  -1.73559036125876e-019
   8  -2.80141862026480e-019   0.00000000000000e+000   2.80141862026480e-019
   9   2.08935352471683e-011   2.08935351786580e-011  -6.85103812046978e-020
  10   2.78750414847452e-019   0.00000000000000e+000  -2.78750414847452e-019
  11  -1.18836768555192e-014  -1.18837079244649e-014  -3.10689457307157e-020
  12  -5.63052317854762e-020   0.00000000000000e+000   5.63052317854762e-020
  13   4.71418358673733e-018   4.76464654243101e-018   5.04629556936791e-020
 ------------------------------------------------------------

В процессе интегрирования дифференциальных уравнений (1), (2) с помощью функции de78d_c нужно следить за тем, чтобы параметры h, k, iniapr, imax удовлетворяли приведенным выше условиям, сформулированным в "Математическом описании" и в "Замечаниях по использованию". Особенно это относится к длине элементарного сегмента h и параметру k, определяющему порядок используемых частичных сумм рядов Чебышёва для решения и его производных. Необходимо заботиться о том, чтобы указанные ряды Чебышёва на каждом сегменте [x, x + h] длиной h быстро сходились, при этом сумма ряда Чебышёва для решения должна хорошо приближаться частичной суммой порядка k + 2, сумма ряда Чебышёва для первой производной должна хорошо приближаться частичной суммой (k+1)-го порядка, а сумма ряда Чебышёва для второй производной должна хорошо приближаться частичной суммой k-го порядка.