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

Текст подпрограммы и версий (Паскаль)
de79e_p.zip

Тексты тестовых примеров ( Паскаль )
tde79e1_p.zip , tde79e2_p.zip , tde79e3_p.zip , tde79e4_p.zip

Подпрограмма: DE79E (модуль DE79E_p)

Назначение

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

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

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

          Y'' = F (X,Y,Y') ,
          Y = ( y₁, ... , y_M ) ,   Y' = ( y'₁, ... , y'_M ) ,
          F = ( f₁ (X, y₁, ... , y_M, y'₁, ... , y'_M), ... , f_M (X, y₁, ... , y_M, y'₁, ... , y'_M ) )

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

procedure DE79E  (F :Proc_F80E; 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) точностью);

M -	количество уравнений в системе (тип: целый);
XN, YN, DYN -	начальные значения аргумента, решения и его производной; в случае системы уравнений (т.е. когда M ≠ 1) YN и DYN представляют одномерные массивы длины M (тип: с расширенной (Extended) точностью);
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 * K² + 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, подбирая для каждой конкретной задачи наиболее подходящий набор их значений.

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

Так как при интегрировании системы уравнений с помощью подпрограммы DE79E используются глобальные записи (структуры данных) с именами _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/y₂ + x²/(y₁y₂²) ,     y₁(0)  =  1 ,       y'₁(0)  =  0 ,  
             y'₂  =  -1/y₁ + x²/(y₁²y₂) ,     y₂(0)  =  1/2 ,     y'₂(0)  =  0 ,    0 ≤ x ≤ 3√2 .
     
        Точное решение системы:
                 y₁(x)  =  e^x*x  ,     y₂(x)  =  (1/2)*e^-x*x .

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

Unit Filis79_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure Filis79(X :Extended; var Y :Array of Extended;
                var DY :Array of Extended; var Z :Array of Extended;
                M :Integer);
implementation
procedure Filis79(X :Extended; var Y :Array of Extended;
                var DY :Array of Extended; var Z :Array of Extended;
                M :Integer);
begin
Z[0] := 1.e0/Y[1]+X*X/(Y[0]*Y[1]*Y[1]);
Z[1] := -1.e0/Y[0]+X*X/(Y[0]*Y[0]*Y[1]);
end;
end.
. . . . . . . . . . . . . . . . . . . . . . . . . . 
function Tde79e1: String;
var
M,K0,IMAX,INIAPR,NCHEB,_i :Integer;
XN,H,XK,H0 :Extended;
YN :Array [0..1] of Extended;
Y :Array [0..1] of Extended;
RAB :Array [0..468] of Extended;
DYN :Array [0..1] of Extended;
DY :Array [0..1] of Extended;
YT :Array [0..1] of Extended;
DYT :Array [0..1] of Extended;
DELY :Array [0..1] of Extended;
RELERR :Array [0..1] of Extended;
DELDY :Array [0..1] of Extended;
RELDY :Array [0..1] of Extended;
const
K :Integer = 10;
label
_42;
begin
M := 2;
XN := 0.0;
YN[0] := 1.e0;
YN[1] := 0.5e0;
DYN[0] := 0.e0;
DYN[1] := 0.e0;
H0 := 0.1e0;
XK := 3.e0*Sqrt(2.e0);
YT[0] := Exp(XK*XK);
YT[1] := 0.5e0*Exp(-XK*XK);
DYT[0] := XK/YT[1];
DYT[1] := -XK/YT[0];
K0 := K;
IMAX := 15;
for INIAPR:=1 to 2 do
 begin
  H := H0;
  NCHEB := 0;
  _i := INIAPR;
  DE79E(Filis79,M,XN,YN,DYN,XK,K0,_i,IMAX,H,Y,DY,RAB);
  DELY[0] := YT[0]-Y[0];
  DELY[1] := YT[1]-Y[1];
  RELERR[0] := DELY[0]/YT[0];
  RELERR[1] := DELY[1]/YT[1];
  DELDY[0] := DYT[0]-DY[0];
  DELDY[1] := DYT[1]-DY[1];
  RELDY[0] := DELDY[0]/DYT[0];
  RELDY[1] := DELDY[1]/DYT[1];
// Операторы вывода на печать: H0, K, INIAPR, IMAX, Y, YT, DELY, RELERR 
. . . . . . . . . . . . . . . . . . . . . . . . . . 
end;

Результаты:
------------------------------------------------------------------------
 H0=  0,1000000000000000    K=10    INIAPR=1    IMAX= 15 
 Y=   6,56599691373309485E+007       7,61498987235623847E-009 
 YT=  6,56599691373305110E+007       7,61498987235631423E-009 
 DELTA Y = -4,37499693362042308E-007       7,57548807082255228E-023 
 RELERR  = -6,66311146822189034E-015       9,94812625860848471E-015 

 DY=  5,57143313154075224E+008      -6,46153317289199974E-008 
 DYT= 5,57143313154069952E+008      -6,46153317289203810E-008 
 DELTA DY = -5,27215888723731041E-006      -3,83592084374395443E-022 
 RELDY    = -9,46284153962262645E-015       5,93654902188961732E-015 
------------------------------------------------------------------------
 H0=  0,1000000000000000    K=10    INIAPR=2    IMAX= 15 
 Y=   6,56599691373309674E+007       7,61498987235623578E-009 
 YT=  6,56599691373305110E+007       7,61498987235631423E-009 
 DELTA Y  = -4,56406269222497940E-007       7,84480644604245406E-023 
 RELERR   = -6,95105823561850235E-015       1,03017949827095797E-014  

 DY=  5,57143313154075418E+008      -6,46153317289199792E-008 
 DYT= 5,57143313154069952E+008      -6,46153317289203810E-008 
 DELTA DY = -5,46628143638372421E-006      -4,01802982547633193E-022 
 RELDY    = -9,81126634265482584E-015       6,21838458143820288E-015 
-----------------------------------------------------------------------

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

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

Результаты:
------------------------------------------------------------------------
 H0=  0,5000000000000000    K=15    INIAPR=1    IMAX= 28 
 Y=   6,56599691373913436E+007       7,61498987234876199E-009 
 YT=  6,56599691373305110E+007       7,61498987235631423E-009 
 DELTA Y  = -6,08325681241694838E-005       7,55223492520443563E-021 
 RELERR   = -9,26478780349952334E-013       9,91759024213585687E-013 

 DY=  5,57143313154651375E+008      -6,46153317288558113E-008 
 DYT= 5,57143313154069952E+008      -6,46153317289203810E-008 
 DELTA DY = -5,81423868425190449E-004      -6,45696927664951980E-020 
 RELDY    = -1,04358044815016215E-012       9,99293682146264429E-013 
------------------------------------------------------------------------
 H0=  0,5000000000000000    K=15    INIAPR=2    IMAX= 28 
 Y=   6,56599691372989531E+007       7,61498987236201987E-009 
 YT=  6,56599691373305110E+007       7,61498987235631423E-009 
 DELTA Y  = 3,15579272864852101E-005      -5,70563660262362987E-021 
 RELERR   = 4,80626593358283711E-013      -7,49263846474181692E-013 

 DY=  5,57143313153673267E+008      -6,46153317289487217E-008 
 DYT= 5,57143313154069952E+008      -6,46153317289203810E-008 
 DELTA DY = 3,96684510633349419E-004       2,83406941262299094E-020 
 RELDY    = 7,11997256841619915E-013      -4,38606339515939480E-013 
----------------------------------------------------------------------------

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

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

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

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

Unit Filis79_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure Filis79(X :Extended; var Y :Array of Extended;
                var DY :Array of Extended; var Z :Array of Extended;
                M :Integer);
implementation
procedure Filis79(X :Extended; var Y :Array of Extended;
                var DY :Array of Extended; var Z :Array of Extended;
                M :Integer);
begin
Z[0] := 1.e0/Y[1]+X*X/(Y[0]*Y[1]*Y[1]);
Z[1] := -1.e0/Y[0]+X*X/(Y[0]*Y[0]*Y[1]);
end;
end.
. . . . . . . . . . . . . . . . . . . . . . . . . . 
function Tde79e3: String;
var
M,K0,IMAX,INIAPR,_i :Integer;
XN,H,XK,H0 :Extended;
YN :Array [0..1] of Extended;
Y :Array [0..1] of Extended;
RAB :Array [0..468] of Extended;
DYN :Array [0..1] of Extended;
DY :Array [0..1] of Extended;
YT :Array [0..1] of Extended;
DYT :Array [0..1] of Extended;
DELY :Array [0..1] of Extended;
RELERR :Array [0..1] of Extended;
DELDY :Array [0..1] of Extended;
const
K :Integer = 10;
label
_42;
begin
M := 2;
XN := -3.e0*Sqrt(2.e0);
YN[0] := Exp(XN*XN);
YN[1] := 0.5e0*Exp(-XN*XN);
DYN[0] := XN/YN[1];
DYN[1] := -XN/YN[0];
YT[0] := 1.e0;
YT[1] := 0.5e0;
DYT[0] := 0.e0;
DYT[1] := 0.e0;
H0 := 0.1e0;
XK := 0.e0;
K0 := K;
IMAX := 14;
for INIAPR:=1 to 2 do
 begin
  H := H0;
  _i := INIAPR;
  DE79E(Filis79,M,XN,YN,DYN,XK,K0,_i,IMAX,H,Y,DY,RAB);
  DELY[0] := YT[0]-Y[0];
  DELY[1] := YT[1]-Y[1];
  RELERR[0] := DELY[0]/YT[0];
  RELERR[1] := DELY[1]/YT[1];
  DELDY[0] := DYT[0]-DY[0];
  DELDY[1] := DYT[1]-DY[1];
// Операторы вывода на печать: H0, K, INIAPR, IMAX, Y, YT, DELY, RELERR 
. . . . . . . . . . . . . . . . . . . . . . . . . . 
end;

Результаты:
------------------------------------------------------------------------
 H0=  0,1000000000000000    K=10    INIAPR=1    IMAX= 14 
 Y=    9,99999999998739906E-001       5,00000000000459418E-001 
 YT=  1,00000000000000000E+000       5,00000000000000000E-001 
 DELTA Y  =  1,26009413407152104E-012      -4,59417497534336805E-013 
 RELERR   =  1,26009413407152104E-012      -9,18834995068673610E-013 

 DY=    -6,68592480370430037E-014      -3,13194616590036987E-014 
 DYT=   0,00000000000000000E+000       0,00000000000000000E+000 
 DELTA DY =  6,68592480370430037E-014       3,13194616590036987E-014 
------------------------------------------------------------------------
 H0=  0,1000000000000000    K=10    INIAPR=2    IMAX= 14 
 Y=     9,99999999998741649E-001       5,00000000000458774E-001 
 YT=   1,00000000000000000E+000       5,00000000000000000E-001 
 DELTA Y  =  1,25835079118827298E-012      -4,58774186175292531E-013 
 RELERR   =  1,25835079118827298E-012      -9,17548372350585062E-013 

 DY=  -6,67760050271186401E-014      -3,12795240555406584E-014 
 DYT= 0,00000000000000000E+000       0,00000000000000000E+000 
 DELTA DY =  6,67760050271186401E-014       3,12795240555406584E-014 
------------------------------------------------------------------------

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

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

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

Unit Filis79_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure Filis79(X :Extended; var Y :Array of Extended;
                var DY :Array of Extended; var Z :Array of Extended;
                M :Integer);
implementation
procedure Filis79(X :Extended; var Y :Array of Extended;
                var DY :Array of Extended; var Z :Array of Extended;
                M :Integer);
begin
Z[0] := 1.e0/Y[1]+X*X/(Y[0]*Y[1]*Y[1]);
Z[1] := -1.e0/Y[0]+X*X/(Y[0]*Y[0]*Y[1]);
end;
end.
. . . . . . . . . . . . . . . . . . . . . . . . . . 
function Tde79e4: String;
var
M,K0,IMAX,INIAPR,_i :Integer;
XN,H,XK,H0 :Extended;
YN :Array [0..1] of Extended;
Y :Array [0..1] of Extended;
RAB :Array [0..468] of Extended;
DYN :Array [0..1] of Extended;
DY :Array [0..1] of Extended;
YT :Array [0..1] of Extended;
DYT :Array [0..1] of Extended;
DELY :Array [0..1] of Extended;
RELERR :Array [0..1] of Extended;
DELDY :Array [0..1] of Extended;
const
K :Integer = 10;
label
_42;
begin
M := 2;
XN := 3.e0*Sqrt(2.e0);
YN[0] := Exp(XN*XN);
YN[1] := 0.5e0*Exp(-XN*XN);
DYN[0] := XN/YN[1];
DYN[1] := -XN/YN[0];
XK := 0.e0;
YT[0] := 1.e0;
YT[1] := 0.5e0;
DYT[0] := 0.e0;
DYT[1] := 0.e0;
H0 := -0.1e0;
K0 := K;
IMAX := 14;
for INIAPR:=1 to 2 do
 begin
  H := H0;
  _i := INIAPR;
  DE79E(Filis79,M,XN,YN,DYN,XK,K0,_i,IMAX,H,Y,DY,RAB);
  DELY[0] := YT[0]-Y[0];
  DELY[1] := YT[1]-Y[1];
  RELERR[0] := DELY[0]/YT[0];
  RELERR[1] := DELY[1]/YT[1];
  DELDY[0] := DYT[0]-DY[0];
  DELDY[1] := DYT[1]-DY[1];
// Операторы вывода на печать: H0, K, INIAPR, IMAX, Y, YT, DELY, RELERR 
. . . . . . . . . . . . . . . . . . . . . . . . . . 
end;

Результаты:
------------------------------------------------------------------------
 H0= -0,1000000000000000    K=10    INIAPR=1    IMAX= 14 
 Y=    9,99999999998739906E-001       5,00000000000459418E-001 
 YT=  1,00000000000000000E+000       5,00000000000000000E-001 
 DELTA Y  =  1,26009413407152104E-012      -4,59417497534336805E-013 
 RELERR   =  1,26009413407152104E-012      -9,18834995068673610E-013 

 DY=  6,68592480370430037E-014       3,13194616590036987E-014 
 DYT=   0,00000000000000000E+000       0,00000000000000000E+000 
 DELTA DY =  -6,68592480370430037E-014      -3,13194616590036987E-014 
------------------------------------------------------------------------
 H0= -0,1000000000000000    K=10    INIAPR=2    IMAX= 14 
 Y=    9,99999999998741649E-001       5,00000000000458774E-001 
 YT=  1,00000000000000000E+000       5,00000000000000000E-001 
 DELTA Y  =  1,25835079118827298E-012      -4,58774186175292531E-013 
 RELERR   =  1,25835079118827298E-012      -9,17548372350585062E-013 

 DY=   6,67760050271186401E-014       3,12795240555406584E-014 
 DYT=  0,00000000000000000E+000       0,00000000000000000E+000 
 DELTA DY =  -6,67760050271186401E-014      -3,12795240555406584E-014 
-----------------------------------------------------------------------

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