Текст подпрограммы и версий (Паскаль)
de89e_p.zip
Тексты тестовых примеров (Паскаль)
tde89e1_p.zip , tde89e2_p.zip , tde89e3_p.zip

Подпрограмма:  DE89E (модуль DE89E_p)

Назначение

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

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

Решается задача Коши для канонической системы 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)  состоит из совокупности частичных сумм рядов Чебышёва, построенных на этих элементарных сегментах. Порядок этих частичных сумм и длина элементарных сегментов задаются пользователем при обращении к подпрограмме. Число NX элементарных сегментов, на которые разбивается промежуток интегрирования, равно [(XK-XN) / H], если длина промежутка интегрирования является целым кратным H. Если длина интервала интегрирования не является целым кратным H, то число элементарных сегментов NX равно [(XK-XN) / H] + 1; в этом случае последний элементарный сегмент считается нестандартным (квадратные скобки означают целую часть числа).

После вычисления коэффициентов Чебышёва для решения задачи Коши и его производных на каждом отдельном элементарном сегменте [xs, xs + H] подпрограмма DE89D передает управление специальной подпрограмме обработки результатов для дальнейшей работы с полученными на этом сегменте коэффициентами Чебышёва ai*[Y], ai*[Y'], ai*[Y'']. Характер и цель этой обработки определяет сам пользователь и выполняет ее с помощью составленной им же подпрограммы, имя которой задается при обращении к подпрограмме DE89D. Например, он может использовать их для каких-то необходимых ему на данном элементарном сегменте вычислений, или запомнить вычисленные коэффициенты в удобном для себя виде, или просто их проигнорировать, если эти коэффициенты на данном элементарном сегменте ему не нужны. Таким образом, в подпрограмме обработки результатов можно, в частности, непосредственно приступить к применению полученных на элементарном сегменте [xs, xs + H] коэффициентов Чебышёва решения и его производных (или, по-другому, немедленно начать использование построенных на сегменте [xs, xs + H] частичных сумм (6), (7), (8) рядов Чебышёва), не дожидаясь окончания процедуры интегрирования системы дифференциальных уравнений (1) на всем промежутке [XN, XK]. Иными словами, подпрограмма DE89D после вычисления коэффициентов Чебышёва на элементарном сегменте [xs, xs + H] = [xs, xs+1] предоставляет возможность приостановить процедуру интегрирования системы (1) на время обработки относящихся к этому сегменту результатов, а затем (после завершения такой обработки) продолжить дальнейшее интегрирование, т.е. перейти к вычислению коэффициентов на следующем элементарном сегменте [xs+1, xs+2] и т.д.

При разбиении промежутка интегрирования на элементарные сегменты решение задачи сводится к определению нескольких наборов коэффициентов 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]. В этом случае погрешности начального приближения для коэффициентов a0*[Φ], a1*[Φ], ... , aK*[Φ] имеют, соответственно, порядки O(H), O(H2), ... , O(HK + 1). Второй способ определения начального приближения в некоторых случаях может привести к более быстрой сходимости итерационного процесса и, тем самым, к меньшему числу выполняемых итераций. Второй способ может быть применен только начиная со второго элементарного сегмента [x0 + H, x0 + 2H]. На начальном элементарном сегменте [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.

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

procedure DE89E(F :Proc_F80E; F2TREA :Proc_F2TREA; var M :Integer; XN :Extended;
                var YN :Array of Extended; var DYN :Array of Extended;
                XK :Extended; var K :Integer; var INIAPR :Integer;
                IMAX :Integer; var H :Extended; var Y :Array of Extended;
                var DY :Array of Extended; var RAB :Array of Extended);

Параметры

F - имя подпрограммы вычисления значений правой части дифференциального уравнения. Первый оператор подпрограммы должен соответствовать процедурному типу:
 
 
       Procedure (X :Extended; var Y :Array of Extended; var DY :Array of Extended;
                         var Z :Array of Extended; M :Integer);

Здесь: X, Y, DY - значения независимой, зависимой переменных и производной решения соответственно. Вычисленное значение правой части должно быть помещено в Z. B случае системы уравнений, т.е. когда M ≠ 1 , параметры Y, DY, Z представляют одномерные массивы длины M (тип параметров X, Y, DY и Z: с расширенной (Extended) точностью);

F2TREA - имя подпрограммы обработки результатов. Первый оператор подпрограммы должен соответствовать процедурному типу:
 
 
     Procedure (S :Integer; XI :Extended; XE :Extended;
                        var Y :Array of Extended; var DY :Array of Extended;
                        var AY :Array of Extended; var ADY :Array of Extended;
                        var AD2Y :Array of Extended; M :Integer;
                        KP1 :Integer; KP2 :Integer; KP3 :Integer);
   Здесь:
S - номер элементарного сегмента S = 1, 2, ... , NX (тип целый);
XI, XE - начало и конец элементарного сегмента с данным номером S: [xs-1, xs-1 + H], x0 = XN, т.е. XI = xs-1, XE = xs-1 + H (тип: с расширенной (Extended) точностью);
Y, DY - значения решения задачи Коши и его первой производной, вычисленные подпрограммой DE89E в конце XE элементарного сегмента [xs-1, xs-1 + H], т.е. в точке xs = xs-1 + H (тип: с расширенной (Extended) точностью);
AY - двумерный массив с измерениями M, KP3 (KP3 = K+3). Переменная с индексом AY (N, I+1) содержит вычисленный подпрограммой DE89E I-й коэффициент Чебышёва для N-й компоненты решения yN(X) на элементарном сегменте [XI, XE], I = 0, 1, ... , K + 2;  N = 1, ... , M (тип: с расширенной (Extended) точностью);
ADY - двумерный массив с измерениями M, KP2 (KP2 = K+2). Переменная с индексом ADY (N, I+1) содержит вычисленный подпрограммой DE89E I-й коэффициент Чебышёва для N-й компоненты первой производной решения y'N(X) на элементарном сегменте [XI, XE], I = 0, 1, ... , K + 1;  N = 1, ... , M (тип: с расширенной (Extended) точностью);
AD2Y - двумерный массив с измерениями M, KP1 (KP1 = K+1). Переменная с индексом AD2Y (N, I+1) содержит вычисленный подпрограммой DE89E I-й коэффициент Чебышёва для N-й компоненты второй производной решения y''N(X) на элементарном сегменте [XI, XE], I = 0, 1, ... , K;  N = 1, ... , M (тип: с расширенной (Extended) точностью);
M - количество уравнений в системе (1) (тип: целый);
KP1 - целый параметр, имеющий значение K + 1;
KP2 - целый параметр, имеющий значение K + 2;
KP3 - целый параметр, имеющий значение K + 3;
    Таким образом, для системы уравнений (т.е. когда M ≠ 1) параметры Y, DY, AY, ADY, AD2Y представляют массивы с регулируемыми измерениями, описатели которых имеют вид Y(M), DY(M), AY (M, KP3), ADY (M, KP2), AD2Y (M, KP1). В случае, когда интегрируется одно скалярное уравнение (т.е. при M = 1), парметры AY, ADY, AD2Y представляют массивы с регулируемыми измерениями, описатели которых могут иметь вид AY (KP3), ADY (KP2), AD2Y (KP1), а параметры Y и DY могут быть переменными (простыми переменными). Данная подпрограмма F2TREA составляется пользователем и выполняет нужную ему обработку содержащихся в AY, ADY, AD2Y коэффициентов (см. "Математическое описание"). Обработка результатов на каждом элементарном сегменте в подпрограмме F2TREA должна заканчиваться оператором exit . При работе подпрограммы F2TREA значения всех ее параметров не должны изменяться.
M - количество уравнений в системе (тип: целый);
XN, YN, DYN - начальные значения аргумента, решения и его производной; в случае системы уравнений (т.е. когда M ≠ 1) YN и DYN представляют одномерные массивы длины M (тип: с расширенной (Extended) точностью);
XK - конец интервала интегрирования. На отрезке [XN, XK] вычисляется приближенное аналитическое решение задачи Коши в виде одной частичной суммы ряда Чебышёва либо в виде совокупности частичных сумм. XK может быть больше, меньше или равно XN (тип: с расширенной (Extended) точностью);
K - порядок частичной суммы смещенного ряда Чебышёва, с помощью которой аппроксимируется вторая производная решения задачи Коши на каждом элементарном сегменте разбиения интервала интегрирования; при этом само решение задачи Коши приближается на каждом элементарном сегменте частичной суммой (K + 2) - го порядка, а его первая производная - частичной суммой (K + 1) - го порядка; K≥2 (см. "Математическое описание" и "Замечания по использованию"; тип: целый);
INIAPR - целый указатель способа выбора начального приближения коэффициентов Чебышёва для второй производной решения на каждом элементарном сегменте:
INIAPR=1 - для первого способа, когда начальное приближение определяется только с использованием значений решения и его первой производной в начале каждого элементарного сегмента;
INIAPR=2 - для второго способа, когда начальное приближение коэффициентов Чебышёва на текущем элементарном сегменте (начиная со второго) определяется через коэффициенты Чебышёва, вычисленные на предыдущем элементарном сегменте, т.е. путем экстраполяции коэффициентов с предыдущего сегмента на следующий (см. "Математическое описание");
IMAX - целая переменная, задающая количество итераций, которое предполагается выполнить в итерационном процессе вычисления коэффициентов Чебышёва для второй производной решения задачи Коши на каждом элементарном сегменте, исходя из некоторого начального приближения, способ определения которого задается параметром INIAPR; IMAX≥1. Для получения максимального порядка точности приближенного решения необходимо выполнить не менее K итераций (см. "Математическое описание" и "Замечания по использованию");
H - переменная с расширенной (Extended) точностью, содержащая значение длины элементарных сегментов, на которые разбивается интервал интегрирования (диаметр разбиения промежутка интегрирования или аналог шага интегрирования для разностных методов). Предполагается, что на каждом элементарном сегменте ряды Чебышёва для решения и его производных являются быстросходящимися рядами. Может задаваться с учетом направления интегрирования, т.е. положительным, если XK > XN, отрицательным, если XK < XN , или без такого учета в виде абсолютной величины;
Y, DY - на выходе из подпрограммы содержат значения решения задачи Коши и его первой производной соответственно, вычисленные подпрограммой при значении аргумента XK. Для системы уравнений (когда M ≠ 1) задаются одномерными массивами длины M. В случае совпадения значений параметров XN и XK значения Y и DY полагаются равными начальным значениям YN и DYN (тип: с расширенной (Extended) точностью);
RAB - одномерный рабочий массив длины 2 * K2 + 10 * K + 7 * M * K + 12 * М +5 (тип: с расширенной (Extended) точностью).

Версии: нет

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

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

Кроме того, используются рабочие подпрограммы DE70EK, DE70EH, DE80EK, DE80EH, DE80E0, DE80EI, DE70EF, DE70EQ, DE71EE, DE70EP, DE71ET, DE71EP, DE71EI, DE71EF, DE71ES, DE70EA, DE70EC.

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

 

Разбиение промежутка интегрирования на элементарные сегменты (шаги) выполняется для того, чтобы на каждом таком сегменте ряды Чебышёва для решения и его первой и второй производных были быстросходящимися рядами. Другими словами, длина элементарных сегментов, задаваемая параметром 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, подбирая для каждой конкретной задачи наиболее подходящий набор их значений.

Для вычисления решения при заданном значении аргумента X [XN, XK] необходимо найти такой элементарный сегмент [xS-1, xS] (это можно сделать в подпрограмме F2TREA), чтобы X  [xS-1, xS], и определить нормализованный аргумент α по формуле α = (X - xS-1) / H = (X - xS-1) / (xS - xS-1). Далее, зная α и S и пользуясь соответствующим набором коэффициентов из массива AY, найти значение решения Y(X) по формуле:

      Y(X) = (y1(X), ... , yM(X)),
                                              K+2
      yN(X) = yN (xS-1 + αH) =  ∑ ' AY (N, I + 1) * TI*(α),     N = 1, 2, ... , M.
                                               I=0

Аналогично можно вычислить первую производную решения Y'(X):

      Y'(X) = (y'1(X), ... , y'M(X)),
                                                K+1
      y'N(X) = y'N (xS-1 + αH) =  ∑ ' ADY (N, I + 1) * TI*(α),     N = 1, 2, ... , M; 
                                                 I=0

и вторую производную решения Y'' = F(X,Y(X),Y'(X)):

      Y''(X) = (y''1(X), ... , y''M(X)),
                                                   K
      y''N(X) = y''N (xS-1 + αH) =  ∑ ' AD2Y (N, I + 1) * TI*(α),     N = 1, 2, ... , M. 
                                                  I=0                                                                             

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

Так как при интегрировании системы уравнений с помощью подпрограммы DE89E используются глобальные записи (структуры данных) с именами _COM70D и _COM80D для хранения промежуточных значений, пользователь не должен портить элементы этих структур. Структуры определены в модуле Lstruct и являются объектами типа:

type                                                         type  
 COM70D = record                                 COM80D = record  
  elm1: Extended;   // WC1                          elm1: Integer;       //JLAS
  elm2: Extended;   // WC2                          elm2: Integer;       //JLAS1
  elm3: Extended;   // WC3                          elm3: Integer;       //JLAS2 
  elm4: Integer;       // LASN                        elm4: Integer;       //JLAS3
  elm5: Extended;   // HD2                           elm5: Extended;    //H2D4
  elm6: Extended;   // HD4                           elm6: Extended;    //H2D8
  elm7: Extended;   // HOLD                        elm7: Extended;    //H2D16
 end;                                                           elm8: Extended;    //H2D32
                                                                   elm9: Extended;    //H2D96
                                                                  end;
var                                                            var
_COM70D : COM70D;                          _COM80D : COM80D;

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

Вычисления в трех примерах на Паскале проводились с 19-20 значащими цифрами.
1) Решается задача Коши

       y''1 = -2qy'2 - ((1 -e3-y1+y'2/(2q) ) / (x+1))2 ,    y1(0) = cos q + 3 ,    y'1(0) = 2q sin q ,

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

Точное решение системы

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

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

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

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

Интервал интегрирования [0, 1] разбивается на два элементарных сегмента длиной H = 0,5.

Unit F88EX1_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure F88EX1(X :Extended; var Y :Array of Extended;
                var DY :Array of Extended; var Z :Array of Extended;
                M :Integer);
implementation
procedure F88EX1(X :Extended; var Y :Array of Extended;
                var DY :Array of Extended; var Z :Array of Extended;
                M :Integer);
var
Q :Real;
begin
Q := 0.5e0;
Z[0] := -2.e0*Q*DY[1]-IntPower((1.e0-Exp(3.e0-Y[0]+0.5e0*DY[1]/Q))/(X+1.e0),2);
Z[1] := 2.e0*Q*DY[0]-IntPower(DY[1]-2.e0*Q*(Y[0]-3.e0),2);
end;
end.
. . . . . . . . . . . . . . . . . . . . . . . . . . 
procedure F2TREA1(S :Integer; XI :Extended; XE :Extended;
                var Y :Array of Extended; var DY :Array of Extended;
                var AY :Array of Extended; var ADY :Array of Extended;
                var AD2Y :Array of Extended; M :Integer; KP1 :Integer;
                KP2 :Integer; KP3 :Integer);
var
_i,L,J :Integer;
label
_20;
begin
//  SOLUTION Y AT THE END OF THE SEGMENT [XI,XE],i.e. Y(XE)]
//  Операторы вывода на печать : Y, AY, ADY, AD2Y, M 
. . . . . . . . . . . . . . . . . . . . . . . . .
end;

Unit Tde89e1_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, F88EX1_p, F2TREAT_p, DE89E_p;
function Tde89e1: String;

implementation

function Tde89e1: String;
var
M,K0,IMAX,INIAP0,INIAPR,_i :Integer;
XN,H,XK,H0,X1,Q :Extended;
RАВ :Array [0..534] of Extended;
YN :Array [0..1] of Extended;
DYN :Array [0..1] of Extended;
Y :Array [0..1] of Extended;
DY :Array [0..1] of Extended;
YT :Array [0..1] of Extended;
DELY :Array [0..1] of Extended;
DELDY :Array [0..1] of Extended;
DYT :Array [0..1] of Extended;
const
K :Integer = 11;
label
_42;
begin
Result := '';  { результат функции }
M := 2;
Q := 0.5e0;
XN := 0.0;
YN[0] := Cos(Q)+3.e0;
YN[1] := -Sin(Q)+2.e0;
DYN[0] := 2.e0*Q*Sin(Q);
DYN[1] := 2.e0*Q*Cos(Q);
ХК := 1.e0;
H0 := 0.5e0;
X1 := (2)*XK-1.e0;
YT[0] := Cos(Q*X1)+3.e0;
YT[1] := Sin(Q*X1)+2.e0;
DYT[0] := -2.e0*Q*Sin(Q*X1);
DYT[1] := 2.e0*Q*Cos(Q*X1);
K0 := K;
IМАХ := 13;
for INIAP0:=1 to 2 do
 begin
  INIAPR := INIAP0;
  H := H0;
  DE89E(F88EX1,F2TREA1,M,XN,YN,DYN,XK,K0,INIAPR,IMAX,H,Y,DY,RAB);
  DELY[0] := YT[0]-Y[0];
  DELY[1] := YT[1]-Y[1];
  DELDY[0] := DYT[0]-DY[0];
  DELDY[1] := DYT[1]-DY[1];
// Операторы вывода на печать: H0, K, INIAPR, IMAX, IORDY,Y, YT, DELY 
. . . . . . . . . . . . . . . . . . . . . . . . . . 
end;

Результаты:
------------------------------------------------------------------------
 COEFFICIENTS AY,ADY AND AD2Y ON   1   SEGMENT:  
 XI=      0,00000000000000000E+000      XE=      5,00000000000000000E-001 
 SOLUTION Y AT THE END OF THE SEGMENT [XI,XE],i.e. Y[XE-1]:
  Y=      4,00000000000000000E+000       2,00000000000000000E+000 

 DERIVATIVE DY AT THE END OF THE SEGMENT [XI,XE],i.e. DY[XE-1]:
 DY=     -3,25260651745651330E-019       9,99999999999999999E-001 

  Number of       Chebyshev coefficients for 1  component
              ----------------------------------------------------------------------------------------------------------
coefficient             for Y                                               for Y'                                             for Y''
 -------------------------------------------------------------------------------------------------------------------
      0            7,90766440054602746E+000       4,87106693080399116E-001      -1,90766440054602746E+000 
      1            6,13690356801086329E-002      -2,40340620085585710E-001      -6,13690356801086359E-002 
      2           -1,50605601386582224E-002      -3,84559236046994694E-003       1,50605601386582203E-002 
      3           -1,60442087410516791E-004       6,28342132945848626E-004       1,60442087410515037E-004 
      4            1,96510520421457702E-005       5,01773738245604203E-006      -1,96510520421467409E-005 
      5            1,25508828077427817E-007      -4,91532402816021935E-007      -1,25508828078307508E-007 
      6           -1,02440704940060995E-008      -2,61574064107065074E-009       1,02440704941364805E-008 
      7           -4,67226936007049560E-011       1,82980896270838956E-010       4,67226930837274720E-011 
      8            2,85969720593749471E-012       7,30200568826793687E-013      -2,85969703050104745E-012 
      9            1,01434358026460271E-014      -3,97249091607051980E-014      -1,01433211873239171E-014 
     10           -4,96631880001673895E-016      -1,26808963720262402E-016       4,96429069726800343E-016 
     11           -1,44117618696049546E-018       5,64123942871364026E-018       1,39591029707508696E-018 
     12            5,87629107157670860E-020       1,45407322611988225E-020  
     13            1,39814733280757909E-022 

  Number of       Chebyshev coefficients for 2  component
              ----------------------------------------------------------------------------------------------------------
coefficient             for Y                                               for Y'                                             for Y''
 -------------------------------------------------------------------------------------------------------------------
      0            3,51289330691960088E+000       1,90766440054602746E+000       4,87106693080399110E-001 
      1            2,40340620085585710E-001       6,13690356801086327E-002      -2,40340620085585715E-001 
      2            3,84559236046994685E-003      -1,50605601386582226E-002      -3,84559236046995161E-003 
      3           -6,28342132945848677E-004      -1,60442087410516896E-004       6,28342132945845355E-004 
      4           -5,01773738245607149E-006       1,96510520421456999E-005       5,01773738245389077E-006 
      5            4,91532402815993164E-007       1,25508828077391696E-007      -4,91532402817042779E-007 
      6            2,61574064106254872E-009      -1,02440704940266486E-008      -2,61574064177707622E-009 
      7           -1,82980896272013395E-010      -4,67226936106424676E-011       1,82980896236352615E-010 
      8           -7,30200578800963773E-013       2,85969720610156533E-012       7,30200418901962023E-013 
      9            3,97249144212797852E-014       1,01434326192138670E-014      -3,97249541475661744E-014 
     10            1,26810909115377371E-016      -4,96632230579199287E-016      -1,26729681436399400E-016 
     11           -5,64421382982633760E-018      -1,44011001632272045E-018       5,62429876976855425E-018 
     12           -1,50011460033616713E-020       5,85864455184224401E-020  
     13            5,63331206907908078E-022 

 COEFFICIENTS AY,ADY AND AD2Y ON   2   SEGMENT:  
 XI=      5,00000000000000000E-001      XE=      1,00000000000000000E+000 
 SOLUTION Y AT THE END OF THE SEGMENT [XI,XE],i.e. Y[XE-1]:
  Y=      3,87758256189037272E+000       2,47942553860420300E+000 

 DERIVATIVE DY AT THE END OF THE SEGMENT [XI,XE],i.e. DY[XE-1]:
 DY=     -4,79425538604203000E-001       8,77582561890372715E-001 

  Number of       Chebyshev coefficients for 1  component
              ----------------------------------------------------------------------------------------------------------
coefficient             for Y                                               for Y'                                             for Y''
 -------------------------------------------------------------------------------------------------------------------
      0            7,90766440054602746E+000      -4,87106693080399117E-001      -1,90766440054602745E+000 
      1           -6,13690356801086330E-002      -2,40340620085585710E-001       6,13690356801086330E-002 
      2           -1,50605601386582224E-002       3,84559236046994686E-003       1,50605601386582228E-002 
      3            1,60442087410516788E-004       6,28342132945848681E-004      -1,60442087410516745E-004 
      4            1,96510520421457713E-005      -5,01773738245605939E-006      -1,96510520421455686E-005 
      5           -1,25508828077428145E-007      -4,91532402816000929E-007       1,25508828077155543E-007 
      6           -1,02440704940057677E-008       2,61574064106641558E-009       1,02440704944685174E-008 
      7            4,67226936009639817E-011       1,82980896275921154E-010      -4,67226940324043730E-011 
      8            2,85969720601102188E-012      -7,30200587567397645E-013      -2,85969698306720241E-012 
      9           -1,01434361205774910E-014      -3,97249087842461103E-014       1,01435719090763044E-014 
     10           -4,96631878183547620E-016       1,26813114181703948E-016       4,96449398517534446E-016 
     11            1,44123297666945910E-018       5,64147043769925507E-018      -1,47722546001149979E-018 
     12            5,87653170593672403E-020      -1,53877652084531228E-020  
     13           -1,47959280850510796E-022 

  Number of       Chebyshev coefficients for 2  component
              ----------------------------------------------------------------------------------------------------------
coefficient             for Y                                               for Y'                                             for Y''
 -------------------------------------------------------------------------------------------------------------------
      0            4,48710669308039912E+000       1,90766440054602745E+000      -4,87106693080399124E-001 
      1            2,40340620085585710E-001      -6,13690356801086332E-002      -2,40340620085585716E-001 
      2           -3,84559236046994686E-003      -1,50605601386582226E-002       3,84559236046994164E-003 
      3           -6,28342132945848677E-004       1,60442087410516671E-004       6,28342132945844927E-004 
      4            5,01773738245606700E-006       1,96510520421456950E-005      -5,01773738245844877E-006 
      5            4,91532402815993064E-007      -1,25508828077473432E-007      -4,91532402817313480E-007 
      6           -2,61574064106385366E-009      -1,02440704940275167E-008       2,61574064048849973E-009 
      7           -1,82980896271949478E-010       4,67226935915437818E-011       1,82980896007322272E-010 
      8            7,30200578562391526E-013       2,85969720165406693E-012      -7,30200637952047687E-013 
      9            3,97249143469390255E-014      -1,01434364492758893E-014      -3,97248985380118112E-014 
     10           -1,26810965047172530E-016      -4,96631325542908905E-016       1,26786395816345992E-016 
     11           -5,64420155747397682E-018       1,44075449791302263E-018       5,60750542090125160E-018 
     12            1,50078593532606524E-020       5,84115148010547042E-020  
     13            5,61649180779372156E-022 
 -------------------------------------------------------------------------------------
 H0=  0,5000000000000000    K=11    INIAPR=1    IMAX=13 
 Y  =
      3,87758256189037272E+000       2,47942553860420300E+000 
 YT =
      3,87758256189037272E+000       2,47942553860420300E+000 
 DELY =
      2,16840434497100887E-019       4,33680868994201774E-019 
 DY =
     -4,79425538604203000E-001       8,77582561890372715E-001 
 DYT=
     -4,79425538604203000E-001       8,77582561890372716E-001 
 DELDY=
      0,00000000000000000E+000       1,40946282423115576E-018 

 *************************************************************************************

 COEFFICIENTS AY,ADY AND AD2Y ON   1   SEGMENT:  
 XI=      0,00000000000000000E+000      XE=      5,00000000000000000E-001 
 SOLUTION Y AT THE END OF THE SEGMENT [XI,XE],i.e. Y[XE-1]:
  Y=      4,00000000000000000E+000       2,00000000000000000E+000 

 DERIVATIVE DY AT THE END OF THE SEGMENT [XI,XE],i.e. DY[XE-1]:
 DY=     -3,25260651745651330E-019       9,99999999999999999E-001 

  Number of       Chebyshev coefficients for 1  component
              ----------------------------------------------------------------------------------------------------------
coefficient             for Y                                               for Y'                                             for Y''
 -------------------------------------------------------------------------------------------------------------------
      0            7,90766440054602746E+000       4,87106693080399116E-001      -1,90766440054602746E+000 
      1            6,13690356801086329E-002      -2,40340620085585710E-001      -6,13690356801086359E-002 
      2           -1,50605601386582224E-002      -3,84559236046994694E-003       1,50605601386582203E-002 
      3           -1,60442087410516791E-004       6,28342132945848626E-004       1,60442087410515037E-004 
      4            1,96510520421457702E-005       5,01773738245604203E-006      -1,96510520421467409E-005 
      5            1,25508828077427817E-007      -4,91532402816021935E-007      -1,25508828078307508E-007 
      6           -1,02440704940060995E-008      -2,61574064107065074E-009       1,02440704941364805E-008 
      7           -4,67226936007049560E-011       1,82980896270838956E-010       4,67226930837274720E-011 
      8            2,85969720593749471E-012       7,30200568826793687E-013      -2,85969703050104745E-012 
      9            1,01434358026460271E-014      -3,97249091607051980E-014      -1,01433211873239171E-014 
     10           -4,96631880001673895E-016      -1,26808963720262402E-016       4,96429069726800343E-016 
     11           -1,44117618696049546E-018       5,64123942871364026E-018       1,39591029707508696E-018 
     12            5,87629107157670860E-020       1,45407322611988225E-020  
     13            1,39814733280757909E-022 

  Number of       Chebyshev coefficients for 2  component
              ----------------------------------------------------------------------------------------------------------
coefficient             for Y                                               for Y'                                             for Y''
 -------------------------------------------------------------------------------------------------------------------
      0            3,51289330691960088E+000       1,90766440054602746E+000       4,87106693080399110E-001 
      1            2,40340620085585710E-001       6,13690356801086327E-002      -2,40340620085585715E-001 
      2            3,84559236046994685E-003      -1,50605601386582226E-002      -3,84559236046995161E-003 
      3           -6,28342132945848677E-004      -1,60442087410516896E-004       6,28342132945845355E-004 
      4           -5,01773738245607149E-006       1,96510520421456999E-005       5,01773738245389077E-006 
      5            4,91532402815993164E-007       1,25508828077391696E-007      -4,91532402817042779E-007 
      6            2,61574064106254872E-009      -1,02440704940266486E-008      -2,61574064177707622E-009 
      7           -1,82980896272013395E-010      -4,67226936106424676E-011       1,82980896236352615E-010 
      8           -7,30200578800963773E-013       2,85969720610156533E-012       7,30200418901962023E-013 
      9            3,97249144212797852E-014       1,01434326192138670E-014      -3,97249541475661744E-014 
     10            1,26810909115377371E-016      -4,96632230579199287E-016      -1,26729681436399400E-016 
     11           -5,64421382982633760E-018      -1,44011001632272045E-018       5,62429876976855425E-018 
     12           -1,50011460033616713E-020       5,85864455184224401E-020  
     13            5,63331206907908078E-022 

 COEFFICIENTS AY,ADY AND AD2Y ON   2   SEGMENT:  
 XI=      5,00000000000000000E-001      XE=      1,00000000000000000E+000 
 SOLUTION Y AT THE END OF THE SEGMENT [XI,XE],i.e. Y[XE-1]:
  Y=      3,87758256189037272E+000       2,47942553860420300E+000 

 DERIVATIVE DY AT THE END OF THE SEGMENT [XI,XE],i.e. DY[XE-1]:
 DY=     -4,79425538604203000E-001       8,77582561890372715E-001 

  Number of       Chebyshev coefficients for 1  component
              ----------------------------------------------------------------------------------------------------------
coefficient             for Y                                               for Y'                                             for Y''
 -------------------------------------------------------------------------------------------------------------------
      0            7,90766440054602746E+000      -4,87106693080399117E-001      -1,90766440054602746E+000 
      1           -6,13690356801086330E-002      -2,40340620085585710E-001       6,13690356801086328E-002 
      2           -1,50605601386582224E-002       3,84559236046994685E-003       1,50605601386582226E-002 
      3            1,60442087410516788E-004       6,28342132945848676E-004      -1,60442087410516860E-004 
      4            1,96510520421457712E-005      -5,01773738245606214E-006      -1,96510520421456703E-005 
      5           -1,25508828077428214E-007      -4,91532402816002453E-007       1,25508828077128438E-007 
      6           -1,02440704940057818E-008       2,61574064106641558E-009       1,02440704944278598E-008 
      7            4,67226936009639817E-011       1,82980896275074121E-010      -4,67226940595094273E-011 
      8            2,85969720599631645E-012      -7,30200587567397645E-013      -2,85969697629093883E-012 
      9           -1,01434361123424484E-014      -3,97249086901313384E-014       1,01435448040219922E-014 
     10           -4,96631877007112971E-016       1,26812521258640870E-016       4,96449398517534446E-016 
     11            1,44122383256377852E-018       5,64147043769925507E-018      -1,45689666927739658E-018 
     12            5,87653170593672403E-020      -1,51760069716395477E-020  
     13           -1,45923143958072574E-022 

  Number of       Chebyshev coefficients for 2  component
              ----------------------------------------------------------------------------------------------------------
coefficient             for Y                                               for Y'                                             for Y''
 -------------------------------------------------------------------------------------------------------------------
      0            4,48710669308039912E+000       1,90766440054602745E+000      -4,87106693080399117E-001 
      1            2,40340620085585710E-001      -6,13690356801086330E-002      -2,40340620085585710E-001 
      2           -3,84559236046994686E-003      -1,50605601386582224E-002       3,84559236046994687E-003 
      3           -6,28342132945848674E-004       1,60442087410516789E-004       6,28342132945848675E-004 
      4            5,01773738245606932E-006       1,96510520421457718E-005      -5,01773738245607030E-006 
      5            4,91532402815994441E-007      -1,25508828077429393E-007      -4,91532402816023044E-007 
      6           -2,61574064106313343E-009      -1,02440704940058314E-008       2,61574064110543434E-009 
      7           -1,82980896271622189E-010       4,67226936010110848E-011       1,82980896256865543E-010 
      8            7,30200578691296978E-013       2,85969720501112550E-012      -7,30200551186411872E-013 
      9            3,97249143867236158E-014      -1,01434352319217743E-014      -3,97248638464884932E-014 
     10           -1,26810949704697730E-016      -4,96630832974836317E-016       1,26785511955879292E-016 
     11           -5,64419540211679671E-018       1,44074445404408286E-018       5,60279149841218419E-018 
     12            1,50077547296258631E-020       5,83624114417935853E-020  
     13            5,61177033094169090E-022 
 -------------------------------------------------------------------------------------
 H0=  0,5000000000000000    K=11    INIAPR=2    IMAX=13 
 Y  =
      3,87758256189037272E+000       2,47942553860420300E+000 
 YT =
      3,87758256189037272E+000       2,47942553860420300E+000 
 DELY =
      2,16840434497100887E-019       2,16840434497100887E-019 
 DY =
     -4,79425538604203000E-001       8,77582561890372715E-001 
 DYT=
     -4,79425538604203000E-001       8,77582561890372716E-001 
 DELDY=
      5,42101086242752217E-020       7,04731412115577882E-019 

 *************************************************************************************

2) Решается задача Коши из примера 1 также с разбиением интервала интегрирования [0, 1] на два элементарных сегмента длиной H = 0,5. В этом примере компоненты решения y1(x) и y2(x) рассматриваются как координаты. Компоненты производной решения y'1(x) и y'2(x) понимаются как компоненты скорости, а компоненты второй производной решения y''1(x) и y''2(x) отождествляются с компонентами ускорения. Для данного примера подпрограмма обработки результатов F2TREA2 устроена следующим образом. Коэффициенты Чебышёва для каждой координаты запоминаются (сохраняются) в отдельном массиве. Всего таких массивов четыре (по два массива на каждом сегменте). Коэффициенты Чебышёва для каждой компоненты скорости также запоминаются в отдельных массивах и таких массивов тоже четыре (по два массива на каждом сегменте). Подобным образом запоминаются в своих массивах компоненты ускорения. Все эти массивы помещаются в глобальную структуру _FIRST_SEG для первого сегмента и в глобальную структуру _SECOND_SEG для второго сегмента. В этих же структурах размещаются значения координат и компоненты скорости в конце каждого сегмента (интервала) и границы каждого сегмента. Поскольку глобальные структуры FIRST_SEG и SECOND_SEG описаны в подпрограмме обработки результатов и в главной (вызывающей) программе, то на выходе из подпрограммы DE89E вычисленные в ней коэффициенты Чебышёва будут доступны также в главной программе. Приводятся фрагмент вызывающей программы, подпрограмма F88EX1 вычисления значений правой части системы, подпрограмма F2TREA2 обработки результатов. Далее представлены значения параметров H, K, INIAPR, IMAX, приближенное решение Y и производня DY, вычисленные в конце интервала XK = xf, точные решение YT и производная DYT (в точке xf) и абсолютные погрешности приближенных значений Y и DY. Затем даются те же результаты, что и в примере 1, а именно: границы элементарных сегментов, значения координат и компоненты скорости в конце каждого сегмента, коэффициенты Чебышёва для координат, компонент скорости и ускорения на каждом сегменте при значении параметра INIAPR = 1.

  
type
 FIRST_SEG = record
  elm1: Array [0..12] of Extended;    // AY1
  elm2: Array [0..12] of Extended;    // AY2
  elm3: Array [0..11] of Extended;    // AV1
  elm4: Array [0..11] of Extended;    // AV2
  elm5: Array [0..1] of Extended;     // AW1
  elm6: Array [0..1] of Extended;     // AW2
  elm7: Array [0..1] of Extended;     // YA
  elm8: Array [0..1] of Extended;     // VA
  elm9: Array [0..1] of Extended;     // INA
 end;
type
 SECOND_SEG = record
  elm1: Array [0..12] of Extended;   // BY1
  elm2: Array [0..12] of Extended;   // BY2
  elm3: Array [0..11] of Extended;   // BV1
  elm4: Array [0..11] of Extended;   // BV2
  elm5: Array [0..1] of Extended;    // BW1
  elm6: Array [0..1] of Extended;    // BW2
  elm7: Array [0..1] of Extended;     // YB
  elm8: Array [0..1] of Extended;     // VB
  elm9: Array [0..1] of Extended;     // INB
 end;
var
_FIRST_SEG: FIRST_SEG;
_SECOND_SEG: SECOND_SEG;

Unit F88EX1_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure F88EX1(X :Extended; var Y :Array of Extended;
                var DY :Array of Extended; var Z :Array of Extended;
                M :Integer);
implementation
procedure F88EX1(X :Extended; var Y :Array of Extended;
                var DY :Array of Extended; var Z :Array of Extended;
                M :Integer);
var
Q :Real;
begin
Q := 0.5e0;
Z[0] := -2.e0*Q*DY[1]-IntPower((1.e0-Exp(3.e0-Y[0]+0.5e0*DY[1]/Q))/(X+1.e0),2);
Z[1] := 2.e0*Q*DY[0]-IntPower(DY[1]-2.e0*Q*(Y[0]-3.e0),2);
end;
end.
. . . . . . . . . . . . . . . . . . . . . . . . . . 
procedure F2TREA2(S :Integer; XI :Extended; XE :Extended;
                var Y :Array of Extended; var DY :Array of Extended;
                var AY :Array of Extended; var ADY :Array of Extended;
                var AD2Y :Array of Extended; M :Integer;
                KP1 :Integer; KP2 :Integer; KP3 :Integer);
var
I :Integer;
label
_10,_15,_20,_25,_30,_35,_40,_45;
begin
{  }
{     AY1 - CHEBYSHEV COEFFICIENTS FOR 1 COORDINATE AT 1 INTERVAL }
{     AY2 - CHEBYSHEV COEFFICIENTS FOR 2 COORDINATE AT 1 INTERVAL }
{  }
{     AV1 - CHEBYSHEV COEFFICIENTS FOR 1 VELOCITY COMPONENT AT 1 INTERVAL }
{     AV2 - CHEBYSHEV COEFFICIENTS FOR 2 VELOCITY COMPONENT AT 1 INTERVAL }
{  }
{     AW1 - CHEBYSHEV COEFFICIENTS FOR 1 ACCELERATION COMPONENT AT 1 INTE }
{     AW2 - CHEBYSHEV COEFFICIENTS FOR 2 ACCELERATION COMPONENT AT 1 INTE }
{  }
{     BY1 - CHEBYSHEV COEFFICIENTS FOR 1 COORDINATE AT 2 INTERVAL }
{     BY2 - CHEBYSHEV COEFFICIENTS FOR 2 COORDINATE AT 2 INTERVAL }
{  }
{     BV1 - CHEBYSHEV COEFFICIENTS FOR 1 VELOCITY COMPONENT AT 2 INTERVAL }
{     BV2 - CHEBYSHEV COEFFICIENTS FOR 2 VELOCITY COMPONENT AT 2 INTERVAL }
{  }
{     BW1 - CHEBYSHEV COEFFICIENTS FOR 1 ACCELERATION COMPONENT AT 2 INTE }
{     BW2 - CHEBYSHEV COEFFICIENTS FOR 2 ACCELERATION COMPONENT AT 2 INTE }
{  }
{      YA - COORDINATES AT THE END OF 1 INTERVAL }
{      VA - VELOCITY    AT THE END OF 1 INTERVAL }
{  }
{      YB - COORDINATES AT THE END OF 2 INTERVAL }
{      VB - VELOCITY    AT THE END OF 2 INTERVAL }
{  }
{     INA - BOUNDARIES OF 1 INTERVAL }
{     INB - BOUNDARIES OF 2 INTERVAL }
{  }
case S of
 1: goto _10;
 2: goto _30;
end;
{  }
{ --------SAVING CHEBYSHEV COEFFICIENTS FOR  COORDINATES, }
{ --------     VELOCITY COMPONENTS AND ACCELERATION COMPONENTS }
{ ---------             AT THE FIRST INTERVAL }
{  };
_10:
for I:=1 to KP3 do
 begin
  _FIRST_SEG.elm1[I-1] := AY[(I-1)*M];
_15:
  _FIRST_SEG.elm2[I-1] := AY[1+(I-1)*M];
 end;
{  }
for I:=1 to KP2 do
 begin
  _FIRST_SEG.elm3[I-1] := ADY[(I-1)*M];
_20:
  _FIRST_SEG.elm4[I-1] := ADY[1+(I-1)*M];
 end;
{  }
for I:=1 to KP1 do
 begin
  _FIRST_SEG.elm5[I-1] := AD2Y[(I-1)*M];
_25:
  _FIRST_SEG.elm6[I-1] := AD2Y[1+(I-1)*M];
 end;
{  }
{ --------SAVING COORDINATES AND VELOCITY AT THE END OF 1 INTERVAL }
_FIRST_SEG.elm7[0] := Y[0];
_FIRST_SEG.elm7[1] := Y[1];
_FIRST_SEG.elm8[0] := DY[0];
_FIRST_SEG.elm8[1] := DY[1];
{  }
{ --------SAVING BOUNDARIES OF 1 INTERVAL }
_FIRST_SEG.elm9[0] := XI;
_FIRST_SEG.elm9[1] := XE;
exit;
{  }
{ --------SAVING CHEBYSHEV COEFFICIENTS FOR  COORDINATES, }
{ --------     VELOCITY COMPONENTS AND ACCELERATION COMPONENTS }
{ ------------           AT THE SECOND INTERVAL }
{  };
_30:
for I:=1 to KP3 do
 begin
  _SECOND_SEG.elm1[I-1] := AY[(I-1)*M];
_35:
  _SECOND_SEG.elm2[I-1] := AY[1+(I-1)*M];
 end;
{  }
for I:=1 to KP2 do
 begin
  _SECOND_SEG.elm3[I-1] := ADY[(I-1)*M];
_40:
  _SECOND_SEG.elm4[I-1] := ADY[1+(I-1)*M];
 end;
{  }
for I:=1 to KP1 do
 begin
  _SECOND_SEG.elm5[I-1] := AD2Y[(I-1)*M];
_45:
  _SECOND_SEG.elm6[I-1] := AD2Y[1+(I-1)*M];
 end;
{  }
{ --------SAVING COORDINATES AND VELOCITY AT THE END OF 2 INTERVAL }
_SECOND_SEG.elm7[0] := Y[0];
_SECOND_SEG.elm7[1] := Y[1];
_SECOND_SEG.elm8[0] := DY[0];
_SECOND_SEG.elm8[1] := DY[1];
{  }
{ --------SAVING BOUNDARIES OF 2 INTERVAL }
_SECOND_SEG.elm9[0] := XI;
_SECOND_SEG.elm9[1] := XE;
end;

Unit Tde89e2_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, F88EX1_p, F2TREAT_p, DE89E_p;
function Tde89e2: String;

implementation

function Tde89e2: String;
var
M,K0,IMAX,INIAPR,_i,KP1,KP2,KP3,J :Integer;
XN,H,XK,H0,Q,X1 :Extended;
YT :Array [0..1] of Extended;
DYT :Array [0..1] of Extended;
DELY :Array [0..1] of Extended;
DELDY :Array [0..1] of Extended;
YN :Array [0..1] of Extended;
DYN :Array [0..1] of Extended;
Y :Array [0..1] of Extended;
DY :Array [0..1] of Extended;
RAB :Array [0..534] of Extended;
const
K :Integer = 11;
begin
Result := '';  { результат функции }
{  }
{     AY1 - CHEBYSHEV COEFFICIENTS FOR 1 COORDINATE AT 1 INTERVAL }
{     AY2 - CHEBYSHEV COEFFICIENTS FOR 2 COORDINATE AT 1 INTERVAL }
{  }
{     AV1 - CHEBYSHEV COEFFICIENTS FOR 1 VELOCITY COMPONENT AT 1 INTERVAL }
{     AV2 - CHEBYSHEV COEFFICIENTS FOR 2 VELOCITY COMPONENT AT 1 INTERVAL }
{  }
{     AW1 - CHEBYSHEV COEFFICIENTS FOR 1 ACCELERATION COMPONENT AT 1 INTE }
{     AW2 - CHEBYSHEV COEFFICIENTS FOR 2 ACCELERATION COMPONENT AT 1 INTE }
{  }
{     BY1 - CHEBYSHEV COEFFICIENTS FOR 1 COORDINATE AT 2 INTERVAL }
{     BY2 - CHEBYSHEV COEFFICIENTS FOR 2 COORDINATE AT 2 INTERVAL }
{  }
{     BV1 - CHEBYSHEV COEFFICIENTS FOR 1 VELOCITY COMPONENT AT 2 INTERVAL }
{     BV2 - CHEBYSHEV COEFFICIENTS FOR 2 VELOCITY COMPONENT AT 2 INTERVAL }
{  }
{     BW1 - CHEBYSHEV COEFFICIENTS FOR 1 ACCELERATION COMPONENT AT 2 INTE }
{     BW2 - CHEBYSHEV COEFFICIENTS FOR 2 ACCELERATION COMPONENT AT 2 INTE }
{  }
{      YA - COORDINATES AT THE END OF 1 INTERVAL }
{      VA - VELOCITY    AT THE END OF 1 INTERVAL }
{      YB - COORDINATES AT THE END OF 2 INTERVAL }
{      VB - VELOCITY    AT THE END OF 2 INTERVAL }
{  }
{     INA - BOUNDARIES OF 1 INTERVAL }
{     INB - BOUNDARIES OF 2 INTERVAL }
{  }
M := 2;
Q := 0.5e0;
XN := 0.0;
YN[0] := Cos(Q)+3.e0;
YN[1] := -Sin(Q)+2.e0;
DYN[0] := 2.e0*Q*Sin(Q);
DYN[1] := 2.e0*Q*Cos(Q);
XK := 1.e0;
H0 := 0.5e0;
X1 := 2.e0*XK-1.e0;
YT[0] := Cos(Q*X1)+3.e0;
YT[1] := Sin(Q*X1)+2.e0;
DYT[0] := -2.e0*Q*Sin(Q*X1);
DYT[1] := 2.e0*Q*Cos(Q*X1);
K0 := K;
IMAX := 13;
INIAPR := 1;
H := H0;
DE89E(F88EX1,F2TREA2,M,XN,YN,DYN,XK,K0,INIAPR,IMAX,H,Y,DY,RAB);
DELY[0] := YT[0]-Y[0];
DELY[1] := YT[1]-Y[1];
DELDY[0] := DYT[0]-DY[0];
DELDY[1] := DYT[1]-DY[1];
//  Операторы вывода на печать: H0, K, INIAPR, Y, YT, DELY; 
//  коэффициенты: AY1, AY2, AV1, AV2, AW1, AW2, BY1, BY2, BV1, BV2, BW1, BW2   
. . . . . . . . . . . . . . . . . . . . . . . . . . 
end;

Результаты:
------------------------------------------------------------------------
 H0=  0,5000000000000000    K=11    INIAPR=1    IMAX=13 
  Y =
      3,87758256189037272E+000       2,47942553860420300E+000 
  YT=
      3,87758256189037272E+000       2,47942553860420300E+000 
 DELY =
      2,16840434497100887E-019       4,33680868994201774E-019 
 DY =
     -4,79425538604203000E-001       8,77582561890372715E-001 
 DYT=
     -4,79425538604203000E-001       8,77582561890372716E-001 
 DELDY =
      0,00000000000000000E+000       1,40946282423115576E-018 

 ****************************************************************

 CHEBYSHEV COEFFICIENTS  AT 1 SEGMENT
 INA=
      0,00000000000000000E+000       5,00000000000000000E-001 

 COORDINATES AND VELOCITY AT THE END OF 1 SEGMENT:
 YA=
      4,00000000000000000E+000       2,00000000000000000E+000 
 VA=
     -3,25260651745651330E-019       9,99999999999999999E-001 

   Number of                      Coefficients for 1 component
              ----------------------------------------------------------------------------------------------------------
coefficient             coordinate                                      velocity                                         acceleration
 -------------------------------------------------------------------------------------------------------------------
      0            7,90766440054602746E+000       4,87106693080399116E-001      -1,90766440054602746E+000 
      1            6,13690356801086329E-002      -2,40340620085585710E-001      -6,13690356801086359E-002 
      2           -1,50605601386582224E-002      -3,84559236046994694E-003       1,50605601386582203E-002 
      3           -1,60442087410516791E-004       6,28342132945848626E-004       1,60442087410515037E-004 
      4            1,96510520421457702E-005       5,01773738245604203E-006      -1,96510520421467409E-005 
      5            1,25508828077427817E-007      -4,91532402816021935E-007      -1,25508828078307508E-007 
      6           -1,02440704940060995E-008      -2,61574064107065074E-009       1,02440704941364805E-008 
      7           -4,67226936007049560E-011       1,82980896270838956E-010       4,67226930837274720E-011 
      8            2,85969720593749471E-012       7,30200568826793687E-013      -2,85969703050104745E-012 
      9            1,01434358026460271E-014      -3,97249091607051980E-014      -1,01433211873239171E-014 
     10           -4,96631880001673895E-016      -1,26808963720262402E-016       4,96429069726800343E-016 
     11           -1,44117618696049546E-018       5,64123942871364026E-018       1,39591029707508696E-018 
     12            5,87629107157670860E-020       1,45407322611988225E-020 
     13            1,39814733280757909E-022 

   Number of                      Coefficients for 2 component
              ----------------------------------------------------------------------------------------------------------
coefficient             coordinate                                      velocity                                         acceleration
 -------------------------------------------------------------------------------------------------------------------
      0            3,51289330691960088E+000       1,90766440054602746E+000       4,87106693080399110E-001 
      1            2,40340620085585710E-001       6,13690356801086327E-002      -2,40340620085585715E-001 
      2            3,84559236046994685E-003      -1,50605601386582226E-002      -3,84559236046995161E-003 
      3           -6,28342132945848677E-004      -1,60442087410516896E-004       6,28342132945845355E-004 
      4           -5,01773738245607149E-006       1,96510520421456999E-005       5,01773738245389077E-006 
      5            4,91532402815993164E-007       1,25508828077391696E-007      -4,91532402817042779E-007 
      6            2,61574064106254872E-009      -1,02440704940266486E-008      -2,61574064177707622E-009 
      7           -1,82980896272013395E-010      -4,67226936106424676E-011       1,82980896236352615E-010 
      8           -7,30200578800963773E-013       2,85969720610156533E-012       7,30200418901962023E-013 
      9            3,97249144212797852E-014       1,01434326192138670E-014      -3,97249541475661744E-014 
     10            1,26810909115377371E-016      -4,96632230579199287E-016      -1,26729681436399400E-016 
     11           -5,64421382982633760E-018      -1,44011001632272045E-018       5,62429876976855425E-018 
     12           -1,50011460033616713E-020       5,85864455184224401E-020 
     13            5,63331206907908078E-022 

 CHEBYSHEV COEFFICIENTS  AT 2 SEGMENT
 INB=
      5,00000000000000000E-001       1,00000000000000000E+000 

 COORDINATES AND VELOCITY AT THE END OF 2 SEGMENT:
 YB=
      3,87758256189037272E+000       2,47942553860420300E+000 
 VB=
     -4,79425538604203000E-001       8,77582561890372715E-001 

   Number of                      Coefficients for 1 component
              ----------------------------------------------------------------------------------------------------------
coefficient             coordinate                                      velocity                                         acceleration
 -------------------------------------------------------------------------------------------------------------------
      0            7,90766440054602746E+000      -4,87106693080399117E-001      -1,90766440054602745E+000 
      1           -6,13690356801086330E-002      -2,40340620085585710E-001       6,13690356801086330E-002 
      2           -1,50605601386582224E-002       3,84559236046994686E-003       1,50605601386582228E-002 
      3            1,60442087410516788E-004       6,28342132945848681E-004      -1,60442087410516745E-004 
      4            1,96510520421457713E-005      -5,01773738245605939E-006      -1,96510520421455686E-005 
      5           -1,25508828077428145E-007      -4,91532402816000929E-007       1,25508828077155543E-007 
      6           -1,02440704940057677E-008       2,61574064106641558E-009       1,02440704944685174E-008 
      7            4,67226936009639817E-011       1,82980896275921154E-010      -4,67226940324043730E-011 
      8            2,85969720601102188E-012      -7,30200587567397645E-013      -2,85969698306720241E-012 
      9           -1,01434361205774910E-014      -3,97249087842461103E-014       1,01435719090763044E-014 
     10           -4,96631878183547620E-016       1,26813114181703948E-016       4,96449398517534446E-016 
     11            1,44123297666945910E-018       5,64147043769925507E-018      -1,47722546001149979E-018 
     12            5,87653170593672403E-020      -1,53877652084531228E-020 
     13           -1,47959280850510796E-022 

   Number of                      Coefficients for 2 component
              ----------------------------------------------------------------------------------------------------------
coefficient             coordinate                                      velocity                                         acceleration
 -------------------------------------------------------------------------------------------------------------------
      0            4,48710669308039912E+000       1,90766440054602745E+000      -4,87106693080399124E-001 
      1            2,40340620085585710E-001      -6,13690356801086332E-002      -2,40340620085585716E-001 
      2           -3,84559236046994686E-003      -1,50605601386582226E-002       3,84559236046994164E-003 
      3           -6,28342132945848677E-004       1,60442087410516671E-004       6,28342132945844927E-004 
      4            5,01773738245606700E-006       1,96510520421456950E-005      -5,01773738245844877E-006 
      5            4,91532402815993064E-007      -1,25508828077473432E-007      -4,91532402817313480E-007 
      6           -2,61574064106385366E-009      -1,02440704940275167E-008       2,61574064048849973E-009 
      7           -1,82980896271949478E-010       4,67226935915437818E-011       1,82980896007322272E-010 
      8            7,30200578562391526E-013       2,85969720165406693E-012      -7,30200637952047687E-013 
      9            3,97249143469390255E-014      -1,01434364492758893E-014      -3,97248985380118112E-014 
     10           -1,26810965047172530E-016      -4,96631325542908905E-016       1,26786395816345992E-016 
     11           -5,64420155747397682E-018       1,44075449791302263E-018       5,60750542090125160E-018 
     12            1,50078593532606524E-020       5,84115148010547042E-020 
     13            5,61649180779372156E-022 

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

3) Решается задача Коши из примера 1 также с разбиением интервала интегрирования [0, 1] на два элементарных сегмента длиной H = 0,5. Система уравнений в примере 1 описывает в трехмерном пространстве переменных x, y1, y2 движение точки по цилиндрической поверхности. Ось кругового цилиндра параллельна оси x и проходит через точку (0, 3, 2). Фазовая траектория этой системы на плоскости y1, y2 -- окружность с центром в точке (3, 2) радиуса 1. При изменении x точка (y1(x), y2(x)) совершает движение по дуге окружности. В примере 3 подпрограмма обработки результатов выполняет следующие действия.
а) По вычисленным подпрограммой DE89E коэффициентам Чебышёва AY1, AY2 находится решение системы y1(x), y2(x) для нескольких равноотстоящих значений аргумента x, принадлежащих элементарному сегменту [XI, XE] (число этих значений задается переменной NI).
б) Определяется расстояние от точки (y1(x), y2(x)) до точки (3, 2); все расстояния запоминаются в массиве R.
в) Вычисляется возможное отклонение этого расстояния от постоянного значения (равного 1). Это отклонение может быть следствием приближенного интегрирования системы уравнений из-за замены точного решения задачи частичной суммой ряда с приближенными коэффициентами. Все отклонения запоминаются в массиве DELTR.
Массивы R и DELTR помещаются в глобальную структуру _DELTR (описана в модуле Lstruct) вместе с переменной NI, представляющей число точек x, взятых на элементарном сегменте. Поскольку глобальная структура _DELTR используется в подпрограмме обработки результатов FTREAT3 и в головной (вызывающей) программе, то на выходе из подпрограммы DE89E вычисленные в подпрограмме FTREAT3 массивы R и DELTR будут доступны также и в главной программе. Приводятся фрагмент вызывающей программы, подпрограмма F88EX1 вычисления значений правой части системы, подпрограмма FTREAT3 обработки результатов. Далее представлены значения параметров H, K, INIAPR, IMAX, приближенное решение Y и производная DY, вычисленные в конце интервала XK = xf, точные решение YT и производная DYT (в точке xf) и абсолютные погрешности приближенных значений Y и DY. Затем даются расстояния от точек (y1(x), y2(x)) до точки (3, 2) на фазовой плоскости и отклонения этих расстояний от постоянного значения, равного 1, для каждого элементарного сегмента. Эти данные приводятся для двух значений параметра INIAPR (1 и 2).

 
type
 DELTR = record
  elm1: Array [0..31] of Extended;   // R
  elm2: Array [0..31] of Extended;   // DELTR
  elm3: Integer;                                // NI 
 end;
var
_DELTR  :DELTR;

Unit F88EX1_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure F88EX1(X :Extended; var Y :Array of Extended;
                var DY :Array of Extended; var Z :Array of Extended;
                M :Integer);
implementation
procedure F88EX1(X :Extended; var Y :Array of Extended;
                var DY :Array of Extended; var Z :Array of Extended;
                M :Integer);
var
Q :Real;
begin
Q := 0.5e0;
Z[0] := -2.e0*Q*DY[1]-IntPower((1.e0-Exp(3.e0-Y[0]+0.5e0*DY[1]/Q))/(X+1.e0),2);
Z[1] := 2.e0*Q*DY[0]-IntPower(DY[1]-2.e0*Q*(Y[0]-3.e0),2);
end;
end.
. . . . . . . . . . . . . . . . . . . . . . . . . . 
procedure FTREAT3(S :Integer; XI :Extended; XE :Extended; var Y :Array of Extended;
                var AY :Array of Extended; var ADY :Array of Extended;
                M :Integer; KP1 :Integer; KP2 :Integer);
var
J :Integer;
HI,X :Extended;
YI :Array [0..1] of Extended;
BK :Array [0..3] of Extended;
label
_10;
begin
HI := (XE-XI)/(_DELTR.elm3);
X := XI;
for J:=1 to _DELTR.elm3 do
 begin
  X := X+HI;
{  }
{ ------CALCULATION OF SOLUTION YI AT POINT X (BY ALGORITHM CLENSHAW) }
{  }
  DE70EC(M,KP1,X,AY,YI,BK);
{  }
  _DELTR.elm1[(J-1)+(S-1)*16] := Sqrt(IntPower(YI[0]-1.e0,2)+IntPower(YI[1]-1.e0,2));
  _DELTR.elm2[(J-1)+(S-1)*16] := 1.e0-_DELTR.elm1[(J-1)+(S-1)*16];
_10:
 end;
end;

function Tde89e3: String;
var
S,M,K0,IMAX,INIAP0,INIAPR,_i,J :Integer;
XN,H,XK,H0,Q,X1 :Extended;
DELY :Array [0..1] of Extended;
DELDY :Array [0..1] of Extended;
YN :Array [0..1] of Extended;
Y :Array [0..1] of Extended;
DY :Array [0..1] of Extended;
RAB :Array [0..534] of Extended;
YT :Array [0..1] of Extended;
DYN :Array [0..1] of Extended;
DYT :Array [0..1] of Extended;
const
K :Integer = 11;
label
_10,_42;
begin
Result := '';  { результат функции }
_DELTR.elm3 := 16;
M := 2;
Q := 0.5e0;
XN := 0.0;
YN[0] := Cos(Q)+3.e0;
YN[1] := -Sin(Q)+2.e0;
DYN[0] := 2.e0*Q*Sin(Q);
DYN[1] := 2.e0*Q*Cos(Q);
XK := 1.e0;
H0 := 0.5e0;
X1 := (2)*XK-1.e0;
YT[0] := Cos(Q*X1)+3.e0;
YT[1] := Sin(Q*X1)+2.e0;
DYT[0] := -2.e0*Q*Sin(Q*X1);
DYT[1] := 2.e0*Q*Cos(Q*X1);
K0 := K;
IMAX := 13;
for INIAP0:=1 to 2 do
 begin
  INIAPR := INIAP0;
  H := H0;
  DE89E(F88EX1,F2TREA3,M,XN,YN,DYN,XK,K0,INIAPR,IMAX,H,Y,DY,RAB);
  DELY[0] := YT[0]-Y[0];
  DELY[1] := YT[1]-Y[1];
  DELDY[0] := DYT[0]-DY[0];
  DELDY[1] := DYT[1]-DY[1];
//  Операторы вывода на печать: H0, K, INIAPR, IMAX,Y,YT, DELY 
//  коэффициенты фазового движения   
. . . . . . . . . . . . . . . . . . . . . . . . .
end;

Результаты:
------------------------------------------------------------------------
 H0=  0,5000000000000000    K=11    INIAPR=1    IMAX=13 
 Y  =
      3,87758256189037272E+000       2,47942553860420300E+000 
 YT =
      3,87758256189037272E+000       2,47942553860420300E+000 
 DELY =
      2,16840434497100887E-019       4,33680868994201774E-019 
 DY =
     -4,79425538604203000E-001       8,77582561890372715E-001 
 DYT=
     -4,79425538604203000E-001       8,77582561890372716E-001 
 DELDY=
      0,00000000000000000E+000       1,40946282423115576E-018 
 ----------------------------------------------------------------
   Number of        Phase motion  for 1  segment
              ---------------------------------------------------
   point                  distance                                 error
 ----------------------------------------------------------------
      1            1,00000000000000000E+000      -1,08420217248550443E-019  
      2            1,00000000000000000E+000       1,08420217248550443E-019  
      3            1,00000000000000000E+000       5,42101086242752217E-020  
      4            1,00000000000000000E+000       2,16840434497100887E-019  
      5            1,00000000000000000E+000       2,16840434497100887E-019  
      6            1,00000000000000000E+000      -1,08420217248550443E-019  
      7            1,00000000000000000E+000       2,16840434497100887E-019  
      8            1,00000000000000000E+000       5,42101086242752217E-020  
      9            1,00000000000000000E+000       5,42101086242752217E-020  
     10            1,00000000000000000E+000       0,00000000000000000E+000  
     11            1,00000000000000000E+000       1,08420217248550443E-019  
     12            1,00000000000000000E+000       0,00000000000000000E+000  
     13            1,00000000000000000E+000       5,42101086242752217E-020  
     14            1,00000000000000000E+000       2,16840434497100887E-019  
     15            1,00000000000000000E+000       0,00000000000000000E+000  
     16            1,00000000000000000E+000       5,42101086242752217E-020  
 
   Number of        Phase motion  for 2  segment
              ---------------------------------------------------
   point                  distance                                 error
 ----------------------------------------------------------------
      1            1,00000000000000000E+000       4,33680868994201774E-019  
      2            1,00000000000000000E+000       1,62630325872825665E-019  
      3            1,00000000000000000E+000       1,08420217248550443E-019  
      4            1,00000000000000000E+000       4,33680868994201774E-019  
      5            1,00000000000000000E+000       1,62630325872825665E-019  
      6            1,00000000000000000E+000       2,71050543121376108E-019  
      7            1,00000000000000000E+000       3,79470760369926552E-019  
      8            1,00000000000000000E+000       3,79470760369926552E-019  
      9            1,00000000000000000E+000       1,08420217248550443E-019  
     10            1,00000000000000000E+000       1,62630325872825665E-019  
     11            1,00000000000000000E+000       2,71050543121376108E-019  
     12            1,00000000000000000E+000       2,71050543121376108E-019  
     13            1,00000000000000000E+000       3,79470760369926552E-019  
     14            1,00000000000000000E+000       3,25260651745651330E-019  
     15            1,00000000000000000E+000       3,79470760369926552E-019  
     16            1,00000000000000000E+000       3,25260651745651330E-019  
 
 ****************************************************************

 H0=  0,5000000000000000    K=11    INIAPR=2    IMAX=13 
 Y  =
      3,87758256189037272E+000       2,47942553860420300E+000 
 YT =
      3,87758256189037272E+000       2,47942553860420300E+000 
 DELY =
      2,16840434497100887E-019       2,16840434497100887E-019 
 DY =
     -4,79425538604203000E-001       8,77582561890372715E-001 
 DYT=
     -4,79425538604203000E-001       8,77582561890372716E-001 
 DELDY=
      5,42101086242752217E-020       7,04731412115577882E-019 
 ----------------------------------------------------------------
   Number of        Phase motion  for 1  segment
              ---------------------------------------------------
   point                  distance                                 error
 ----------------------------------------------------------------
      1            1,00000000000000000E+000      -1,08420217248550443E-019  
      2            1,00000000000000000E+000       1,08420217248550443E-019  
      3            1,00000000000000000E+000       5,42101086242752217E-020  
      4            1,00000000000000000E+000       2,16840434497100887E-019  
      5            1,00000000000000000E+000       2,16840434497100887E-019  
      6            1,00000000000000000E+000      -1,08420217248550443E-019  
      7            1,00000000000000000E+000       2,16840434497100887E-019  
      8            1,00000000000000000E+000       5,42101086242752217E-020  
      9            1,00000000000000000E+000       5,42101086242752217E-020  
     10            1,00000000000000000E+000       0,00000000000000000E+000  
     11            1,00000000000000000E+000       1,08420217248550443E-019  
     12            1,00000000000000000E+000       0,00000000000000000E+000  
     13            1,00000000000000000E+000       5,42101086242752217E-020  
     14            1,00000000000000000E+000       2,16840434497100887E-019  
     15            1,00000000000000000E+000       0,00000000000000000E+000  
     16            1,00000000000000000E+000       5,42101086242752217E-020  
 
   Number of        Phase motion  for 2  segment
              ---------------------------------------------------
   point                  distance                                 error
 ----------------------------------------------------------------
      1            1,00000000000000000E+000       3,79470760369926552E-019  
      2            1,00000000000000000E+000       5,42101086242752217E-020  
      3            1,00000000000000000E+000       0,00000000000000000E+000  
      4            1,00000000000000000E+000       3,25260651745651330E-019  
      5            1,00000000000000000E+000       5,42101086242752217E-020  
      6            1,00000000000000000E+000       1,62630325872825665E-019  
      7            1,00000000000000000E+000       2,71050543121376108E-019  
      8            1,00000000000000000E+000       2,71050543121376108E-019  
      9            1,00000000000000000E+000       0,00000000000000000E+000  
     10            1,00000000000000000E+000       1,08420217248550443E-019  
     11            1,00000000000000000E+000       1,62630325872825665E-019  
     12            1,00000000000000000E+000       1,62630325872825665E-019  
     13            1,00000000000000000E+000       3,25260651745651330E-019  
     14            1,00000000000000000E+000       2,71050543121376108E-019  
     15            1,00000000000000000E+000       2,71050543121376108E-019  
     16            1,00000000000000000E+000       2,16840434497100887E-019  
 
 ****************************************************************