|
Текст подпрограммы и версий de65e_p.zip |
Тексты тестовых примеров tde65e_p.zip |
Вычисление коэффициентов разложения решения задачи Коши для канонической системы обыкновенных дифференциальных уравнений второго порядка и коэффициентов разложения его первой и второй производных в ряды по смещенным многочленам Чебышёва первого рода с встроенной пошаговой обработкой результатов (получение аналитического выражения для приближенного решения задачи Коши и его производных в виде частичных сумм смещенных рядов Чебышёва с контролем точности приближенного решения и его первой производной).
Решается задача Коши для канонической системы 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] в равномерно сходящиеся ряды по смещенным многочленам Чебышёва первого рода.
Приводимое здесь описание метода решения преследует следующие цели. Во-первых, дать пользователю общее представление об используемом подходе к интегрированию обыкновенных дифференциальных уравнений. Во-вторых, на основе данного знания помочь пользователю яснее усвоить смысл и назначение формальных параметров. С более подробным изложением различных аспектов математического метода решения можно ознакомиться по приводимому в конце данного раздела списку литературы.
Интервал интегрирования [XN, XK] автоматически разбивается на элементарные сегменты переменной длины H: [xs, xs + H], x0 = XN, s = 0, 1, ... . На каждом элементарном сегменте решение исходной задачи Коши приближенно представляется в виде (K + 2)-й частичной суммы смещенного ряда Чебышёва
K+2 1
(3) Y(xs + αH) ≈ UK+2(xs + αH) = ∑'ai*[UK+2]Ti*(α), 0 ≤ α ≤ 1, ai*[UK+2] = 2/π ∫ UK+2(xs+αH)Ti*(α) / √(α(1-α))dα,
i=0 0
а его производные - в виде частичных сумм (K + 1)-го и К-го порядков
K+1 1
(4) Y'(xs + αH) ≈ U'K+2(xs + αH) = ∑'ai*[U'K+2]Ti*(α), 0 ≤ α ≤ 1, ai*[U'K+2] = 2/π ∫ U'K+2(xs + αH)Ti*(α) / √(α(1-α))dα,
i=0 0
K 1
(5) Y''(xs + αH) ≈ U''K+2(xs + αH) = ∑'ai*[U''K+2]Ti*(α), 0 ≤ α ≤ 1, ai*[U''K+2] = 2/π ∫ U''K+2(xs + αH)Ti*(α) / √(α(1-α))dα.
i=0 0
Здесь: штрих у знака суммы означает, что слагаемое с индексом 0 берется с дополнительным множителем 1/2; Ti* (α) - смещенный многочлен Чебышёва первого рода на [0, 1]: Ti* (α) = Ti (2α-1); Ti (t) - многочлен Чебышёва первого рода i-го порядка на [-1, 1].
В этом случае аналитическое решение задачи (1), (2) состоит из совокупности частичных сумм рядов Чебышёва, построенных на этих элементарных сегментах. Порядок этих частичных сумм задается пользователем при обращении к подпрограмме.
Коэффициенты ai*[UK+2] ряда Чебышёва для решения UK+2(xs + αH) и коэффициенты ai*[U'K+2] ряда Чебышёва для его производной U'K+2(xs + αH) на сегменте [xs, xs + H] выражаются через коэффициенты ai*[Φ] = ai*[U''K+2(xs+αH)] ряда Чебышёва второй производной решения Φ(α) = F(xs + αH, UK+2(xs + αH), U'K+2(xs + αH)), 0 ≤ α ≤ 1, на [xs, xs + H], которые, в свою очередь, вычисляются приближенно итерационным способом, исходя из некоторого начального приближения, с помощью квадратурной формулы Маркова на [xs, xs + H] с K + 1 узлом. При этом один из узлов квадратурной формулы совпадает с xs, а остальные K узлов лежат внутри интервала (xs, xs + H). Количество итераций, которое предписывается выполнить в этом итерационном процессе, одинаково для всех сегментов и задается при обращении к подпрограмме. Если при выбранном H ряды Чебышёва (3), (4), (5) для Y(X) = Y(xs + αH), 0 ≤ α ≤ 1, и его первой и второй производных на элементарном сегменте [xs, xs + H] быстро сходятся, то для того, чтобы приближенное решение в конце одного такого сегмента имело максимальный порядок точности относительно H, необходимо выполнить не менее K итераций; при этом погрешность приближенного решения в конце элементарного сегмента является величиной порядка O(HK + 3) , а погрешность приближенного значения производной Y' - величиной порядка O(HK + 2) при H --> 0. Если H подобрано достаточно малым (или, вернее сказать, выбрано довольно удачным), то хорошая точность приближенного решения может быть получена и при меньшем числе итераций.
Начальное приближение коэффициентов ai*[Φ] ряда Чебышёва для второй производной на сегменте [xs, xs + H] может вычисляться двумя способами. В первом способе начальное приближение определяется только с использованием значения решения UK+2(X) и его первой производной U'K+2(X) в узле xs. При этом погрешность начального приближения для всех коэффициентов a0*[Φ], a1*[Φ], ... , aK*[Φ] является величиной O(H2) при H --> 0. Во втором способе начальное приближение определяется через коэффициенты ряда Чебышёва производной Φ(α) на предыдущем элементарном сегменте [xs - 1, xs]. В этом случае погрешности начального приближения для коэффициентов a0*[Φ], a1*[Φ], ... , aK*[Φ] имеют, соответственно, порядки O(H), O(H2), ... , O(HK + 1). Второй способ определения начального приближения в некоторых случаях может привести к более быстрой сходимости итерационного процесса и, тем самым, к меньшему числу выполняемых итераций. Второй способ может быть применен только начиная со второго элементарного сегмента [x1, x2]. На начальном элементарном сегменте [x0, x0 + H] всегда применяется исключительно первый способ. Способ выбора начального приближения задается пользователем при обращении к подпрограмме.
На каждом частичном сегменте [xs, xs + H] = [xs, xs+1] для оценки погрешности решения UK+2(xs + αH) и его производной U'K+2(xs + αH), полученных при некотором заданном значении K = K1, вычисляются на этом же сегменте специальным образом второе приближенное решение вместе с его производной при K = K2 > K1, которые также представляются указанными выше частичными суммами при K = K2. Значение K2 задается пользователем при обращении к подпрограмме. Второе решение UK2+2, как имеющее более высокий порядок точности O(HK2 + 3), используется для оценки погрешности первого решения UK1+2 на каждом элементарном сегменте. Производная второго решения U'K2+2, имеющая порядок точности O(HK2 + 2), используется для оценки погрешности производной первого решения U'K1+2. Заданные значения K1 и K2 используются для всех сегментов, на которые автоматически разбивается интервал интегрирования [XN, XK]. Второе приближенное решение на сегменте [xs, xs+1] вычисляется аналогично первому решению итерационным способом. Единственным отличием при этом является то, что начальным приближением в этом дополнителном итерационном процессе служит уже вычисленное на отрезке [xs, xs+1] первое решение UK1+2(xs + αH). Количество итераций, которое потребуется выполнить в этом дополнительном итерационном процессе, одинаково для всех сегментов и задается при обращении к подпрограмме.
Длина элементарного сегмента H выбирается в подпрограмме из условия, чтобы на каждый элементарный сегмент [xs, xs + H] приходилась приблизительно одинаковая погрешность. Поэтому длина элементерного сегмента H - величина переменная, определяемая автоматически подпрограммой. Она изменяется от одного элементарного сегмента к другому элементарному сегменту. При обращении к подпрограмме задается только начальное значение длины первого элементерного сегмента [XN, XN + H] = [x0, x0 + H], которое может быть уменьшено подпрограммой, если на нем не будет достигнута заданная точность.
Для системы уравнений (M > 1) с проверкой на точность могут вычисляться либо все компоненты решения, либо некоторые из них (в частности, одна компонента). В последнем случае номера проверяемых на точность компонент задаются при обращении к подпрограмме. Совершенно аналогично вычисляются с проверкой на точность компоненты производной решения.
Для оценки уклонения приближенного решения от точного (т.е. для оценки абсолютной погрешности приближенного решения) в подпрограмме используются два способа, или две формулы. Первый способ опирается на первую формулу для абсолютной погрешности приближенного решения. Она является асимптотической (т.е. справедливой при H --> 0) и представляет собою разность двух частичных сумм смещенного ряда Чебышёва на сегменте [xs, xs + H] с приближенными коэффициентами, а именно разность приближенных частичных сумм (K2+2)-го и (K1+2)-го порядков в точке xs + H. Второй способ опирается на несколько завышенную оценку абсолютной погрешности, которая представляет сумму абсолютных величин (модулей) разностей приближенных коэффициентов этих частичных сумм и поэтому несколько завышена по сравнению с первой. Назовем эту формулу второй. Более того, при использовании этой, второй, формулы в том случае, когда точность решения оценивается по относительной погрешности или по мере погрешности, то дополнительно используется также завышенная оценка сверху для относительной погрешности и завышенная оценка сверху для меры погрешности. Второй способ оценки погрешности, хотя и является более точным, накладывает более сильное ограничение на размер элементарного сегмента и приводит к меньшей длине сегментов. Используемый способ оценки абсолютной погрешности задается при обращении к подпрограмме. Совершенно аналогичный прием применяется также для оценки уклонения производной приближенного решения от производной точного решения.
Если для выбранной подпрограммой DE65D длины H сегмента [xs, xs + H] точность приближенного решения UK1+2 или его производной U'K1+2 не достигается, то данный сегмент исключается из рассмотрения и программа сокращает его длину H до тех пор, пока на новом, сокращенном, сегменте не будет достигнута заданная точность вновь вычисленного приближенного решения UK1+2 и его производной. При выполнении условия достижения на сегменте заданной точности приближенного решения UK1+2 и его производной этот сегмент принимается в качестве элементарного (частичного) сегмента, на который продлено решение дифференциального уравнения. В качестве такого решения на сегменте [xs, xs + H] принимается частичная сумма порядка K1 + 2 с коэффициентами, равными коэффициентам Чебышёва ai*[UK2+2], i = 0, 1, ... , K1 + 2, второго приближенного решения UK2+2, поскольку эти коэффициенты имеют более высокий порядок точности по сравнению с коэффициентами Чебышёва ai*[UK1+2] первого приближенного решения. В качестве значения решения в конце xs + H сегмента [xs, xs + H] принимается значение второго приближенного решения в точке xs + H, т.е. UK2+2(xs + H), как имеющее более высокий порядок точности O(HK2+3) по сравнению с первым приближенным решением UK1+2(xs + H). В качестве значения производной решения в конце xs + H сегмента [xs, xs + H] принимается значение производной второго приближенного решения в точке xs + H, т.е. U'K2+2(xs + H), поскольку оно имеет более высокий порядок точности O(HK2+2) по сравнению с производной первого приближенного решения U'K1+2(xs + H). В качестве коэффициентов Чебышёва производной решения на сегменте [xs, xs + H] принимаются по той же причине коэффициенты Чебышёва производной второго приближенного решения ai*[U'K2+2], i = 0, 1, ... , K1 + 1, а в качестве коэффициентов Чебышёва второй производной решения принимаются коэффициенты Чебышёва второй производной второго приближенного решения ai*[U''K2+2], i = 0, 1, ... , K1.
После вычисления коэффициентов Чебышёва для решения задачи Коши и его производных на каждом отдельном элементарном сегменте [xs, xs + H] подпрограмма DE65D передает управление специальной подпрограмме обработки результатов для дальнейшей работы с полученными на этом сегменте коэффициентами Чебышёва ai*[Y], ai*[Y'], ai*[Y'']. Характер и цель этой обработки определяет сам пользователь и выполняет ее с помощью составленной им же подпрограммы, имя которой задается при обращении к подпрограмме DE65D. Например, он может использовать их для каких-то необходимых ему на данном элементарном сегменте вычислений, или запомнить вычисленные коэффициенты в удобном для себя виде, или просто их проигнорировать, если эти коэффициенты на данном элементарном сегменте ему не нужны. Таким образом, в подпрограмме обработки результатов можно, в частности, непосредственно приступить к применению полученных на элементарном сегменте [xs, xs + H] коэффициентов Чебышёва решения и его производных (или, по-другому, немедленно начать использование построенных на сегменте [xs, xs + H] частичных сумм (3), (4), (5) рядов Чебышёва), не дожидаясь окончания процедуры интегрирования системы дифференциальных уравнений (1) на всем промежутке [XN, XK]. Иными словами, подпрограмма DE65D после вычисления коэффициентов Чебышёва на элементарном сегменте [xs, xs + H] = [xs, xs+1] предоставляет возможность приостановить процедуру интегрирования системы (1) на время обработки относящихся к этому сегменту результатов, а затем (после завершения такой обработки) продолжить дальнейшее интегрирование, т.е. перейти к вычислению коэффициентов на следующем элементарном сегменте [xs+1, xs+2] и т.д.
В дальнейшем при описании параметров подпрограммы коэффициенты ряда Чебышёва будем называть коэффициентами Чебышёва.
Залеткин С.Ф. Численное интегрирование обыкновенных дифференциальных уравнений с использованием ортогональных разложений // Математическое моделирование. 2010. 22. № 1. 69 - 85.
Арушанян О.Б., Волченскова Н.И., Залеткин С.Ф. О применении ортогональных разложений для приближенного интегрирования обыкновенных дифференциальных уравнений // Вестник Московского университета. Серия 1. Математика. Механика. 2010. № 4. 40 - 43.
Арушанян О.Б., Волченскова Н.И., Залеткин С.Ф. Вычисление коэффициентов разложения решения задачи Коши в ряд по многочленам Чебышёва. Вестник Московского университета. Серия 1. Математика. Механика. 2012. № 5, 24 - 30.
Арушанян О.Б., Залеткин С.Ф. Обоснование одного подхода к применению ортогональных разложений для приближенного интегрирования канонических систем обыкновенных дифференциальных уравнений второго порядка // Вестник Московского университета. Серия 1. Математика. Механика. 2018. № 3. 29 - 33.
Арушанян О.Б., Залеткин С.Ф. К теории вычисления ортогонального разложения решения задачи Коши для обыкновенных дифференциальных уравнений второго порядка // Вычислительные методы и программирование. 2018. 19. 178 - 184.
Арушанян О.Б., Залеткин С.Ф. Приближенное интегрирование канонических систем обыкновенных дифференциальных уравнений второго порядка методом рядов Чебышёва с оценкой погрешности решения и его производной // Вестник Московского университета. Серия 1. Математика. Механика. 2022. № 4. 27 - 34.
Залеткин С.Ф. Приближенное интегрирование обыкновенных дифференциальных уравнений методом рядов Чебышёва с контролем точности // Математическое моделирование. 2022. 34. № 6. 53 - 74.
procedure DE65E(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 IU :Array of Integer;
var EPS :Array of Extended; var THRESH :Array of Extended;
var MEXACT :Array of Integer;
var NUMBES :Array of Integer; var MESTER :Integer;
var H :Extended; var HMIN :Extended; var HMAX :Extended;
K2 :Integer; var IMAX2 :Integer; var NATTEM :Integer;
var Y :Array of Extended; var DY :Array of Extended;
var RAB :Array of Extended; var IERR :Integer);
Параметры
| 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. Для скалярного случая при M = 1 параметры Y и Z представляют собой массивы длины 1, в случае системы уравнений, т.е. когда M ≠ 1 , параметры Y и Z представляют массивы длины M (тип параметров X, Y и Z: Extended); |
| F2TREA - |
имя подпрограммы-функции обработки результатов. Первый оператор подпрограммы
должен соответствовать процедурному типу: |
Function (S :Integer; XI :Extended; XE :Extended; Y :Extended; DY :Extended;
var AY :Array of Extended; var ADY :Array of Extended;
var AD2Y :Array of Extended; M :Integer;
KP1 :Integer; KP2 :Integer; KP3 :Integer): String;
|
| Здесь: |
| S - |
номер элементарного сегмента: S = 1, 2, ... (тип: целый); |
| XI, XE - |
начало и конец элементарного сегмента с данным номером S: [xs-1, xs-1 + H],
x0 = XN, т.е. XI = xs-1, XE = xs-1 + H
(тип: с расширенной (Extended) точностью); |
| Y, DY - |
значения решения задачи Коши и его первой производной, вычисленные подпрограммой DE65E в конце XE
элементарного сегмента [xs-1, xs-1 + H], т.е. в точке
xs = xs-1 + H (тип: с расширенной (Extended) точностью); |
| AY - |
двумерный массив с измерениями M, KP3 (KP3 = K+3). Переменная с индексом AY (N, I+1) содержит
вычисленный подпрограммой DE65E 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) содержит
вычисленный подпрограммой DE65E 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) содержит
вычисленный подпрограммой DE65E 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. |
| Здесь K - порядок частичной суммы смещенного ряда Чебышёва, с помощью которой аппроксимируется вторая производная первого приближенного решения задачи Коши на каждом элементарном сегменте разбиения интервала интегрирования; при этом само первое приближенное решение задачи Коши представляется на каждом элементарном сегменте частичной суммой (K + 2)-го порядка, а его первая призводная - частичной суммой (K + 1)-го порядка (см. "Математическое описание" и "Замечания по использованию"). | |
| Таким образом, для системы уравнений (т.е. когда 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 должна заканчиваться оператором RETURN. При работе подпрограммы F2TREA значения всех ее параметров не должны изменяться. |
| M - | количество уравнений в системе (тип: целый); |
| XN, YN, DYN - | начальные значения аргумента, решения и его производной; в случае системы уравнений (т.е. когда M ≠ 1) YN и DYN представляют одномерные массивы длины M (тип: с расширенной (Extended) точностью); |
| XK - | конец интервала интегрирования. На отрезке [XN, XK] вычисляется приближенное аналитическое решение задачи Коши в виде одной частичной суммы ряда Чебышёва либо в виде совокупности частичных сумм. Это же относится и к производным приближенного аналитического решения. XK может быть больше, меньше или равно XN (тип: с расширенной (Extended) точностью); |
| K - | порядок частичной суммы смещенного ряда Чебышёва, с помощью которой аппроксимируется вторая производная первого приближенного решения задачи Коши на каждом элементарном сегменте разбиения интервала интегрирования; при этом само первое приближенное решение задачи Коши приближается на каждом элементарном сегменте частичной суммой (K + 2)-го порядка, а его производная - частичной суммой (K + 1)-го порядка; K≥2. Этот параметр относится только к первому приближенному решению, имеющему меньший порядок точности O(HK + 3), и его производной (см. "Математическое описание" и "Замечания по использованию"; тип: целый); |
| INIAPR - |
целый указатель способа выбора начального
приближения коэффициентов Чебышёва для второй производной первого приближенного
решения на каждом элементарном сегменте: |
| INIAPR=1 - |
для первого способа, когда начальное
приближение определяется только с
использованием значения решения и его первой производной в
начале каждого элементарного сегмента; |
| INIAPR=2 - | для второго способа, когда начальное приближение коэффициентов Чебышёва на текущем элементарном сегменте (начиная со второго) определяется через коэффициенты Чебышёва, вычисленные на предыдущем элементарном сегменте, т.е. путем экстраполяции коэффициентов с предыдущего сегмента на следующий (см. "Математическое описание"); |
| IMAX - | целая переменная, задающая количество итераций, которое предполагается выполнить в итерационном процессе вычисления коэффициентов Чебышёва для второй производной первого приближенного решения задачи Коши на каждом элементарном сегменте, исходя из некоторого начального приближения, способ определения которого задается параметром INIAPR; IMAX ≥ 1. Для получения максимального порядка точности первого приближенного решения необходимо выполнить не менее K итераций. Этот параметр относится только к первому приближенному решению, имеющему меньший порядок точности O(HK + 3), и его производной. Если правая часть дифференциального уравнения не зависит от переменных Y и Y', т.е. дифференциальное уравнение имеет вид Y'' = F(X), то число итераций при вычислении первого приближенного решения можно положить равным 1. В этом случае параметр IMAX = 1 (см. "Математическое описание" и "Замечания по использованию"); |
| IU - | целый указатель типа погрешности численного решения и его производной; задается одномерным массивом длины 2. Переменная с индексом IU(1) задает тип погрешности производной решения, а переменная с индексом IU(2) задает тип погрешности решения: |
| IU(1)=1 - | компоненты производной решения проверяются на точность по абсолютной погрешности; |
| IU(1)=2 - | компоненты производной решения проверяются на точность по относительной погрешности; |
| IU(1)=3 - | компоненты производной решения проверяются на точность по мере погрешности; |
| IU(2)=1 - | компоненты решения проверяются на точность по абсолютной погрешности; |
| IU(2)=2 - | компоненты решения проверяются на точность по относительной погрешности; |
| IU(2)=3 - | компоненты решения проверяются на точность по мере погрешности. |
| Контроль точности по мере погрешности заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше определенной наперед заданной положительной константы THRESH (называемой границей перехода), то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Таким образом, контроль точности по мере погрешности состоит в том, что на тех участках интервала интегрирования, где абсолютная величина компоненты решения меньше значения THRESH, контроль точности ведется для нее по абсолютной погрешности, а там, где абсолютная величина этой компоненты решения равна значению THRESH или превосходит значение THRESH, контроль точности для нее ведется по относительной погрешности; аналогично выполняется контроль точности по мере погрешности и для производной решения (см. "Замечания по использованию"); | |
| EPS - | допустимая погрешность, с которой требуется вычислить проверяемые на точность компоненты решения и его производной; тип погрешности специфицируется с помощью параметра IU; задается одномерным массивом длины 2. Переменная с индексом EPS(1) задает точность для производной решения, а переменная с индексом EPS(2) задает точность для решения (см. "Замечания по использованию"; тип: с расширенной (Extended) точностью); |
| THRESH - | граница перехода, используемая при оценке меры погрешности решения и его производной; задается одномерным массивом длины 2. Переменная с индексом THRESH(1) задает границу перехода для производной решения, а переменная THRESH(2) задает границу перехода для решения (см. "Замечания по использованию"; тип: с расширенной (Extended) точностью); |
| MEXACT - | количество компонент численного решения и его производной, проверяемых на точность; задается одномерным массивом длины 2. Переменная с индексом MEXACT(1) показывает, сколько компонент производной надо проверять на точность; переменная с индексом MEXACT(2) показывает, сколько компонент решения надо проверять на точность. Нулевое значение переменной с индексом MEXACT(1) или MEXACT(2) сообщает о том, что соответственно производная или решение на точность не проверяется. Например, элементы массива MEXACT/O, M/ сообщают о том, что производная на точность вообще не проверяется, а проверяются на точность все компоненты решения (здесь M - количество уравнений в системе). Элементы массива MEXACT/M, O/ указывают, что на точность надо проверять все компоненты производной, а решение на точность проверять не надо. Элементы массива MEXACT/M, M/ сообщают, что все компоненты производной и все компоненты решения надо проверять на точность. Элементы массива MEXACT/1, 1/ уведомляют, что на точность проверяется какая-то одна компонента производной и какая-то одна компонента решения (тип: целый); |
| NUMBES - | одномерный массив, содержащий номера проверяемых на точность компонент производной и номера проверяемых на точность компонент решения. Первые MEXACT(1) элементов массива NUMBES содержат номера проверяемых на точность компонент производной, а следующие MEXACT(2) элементов этого массива NUMBES содержат номера проверяемых на точность компонент решения. Следовательно, NUMBES представляет одномерный массив длины MEXACT(1) + MEXACT(2). Если MEXACT(1) = 0, то массив NUMBES содержит только номера проверяемых на точность компонент решения; если MEXACT(2) = 0, то массив NUMBES содержит только номера проверяемых на точность компонент производной. В том случае, если все компоненты производной и все компоненты решения проверяются на точность (этому случаю соответствуют значения MEXACT(1) = MEXACT(2) = M), этот массив в подпрограмме не используется. Таким образом, массив NUMBES используется в подпрограмме только в том случае, когда число проверяемых на точность компонент производной или число проверяемых на точность компонент решения меньше числа уравнений M (тип: целый); |
| MESTER - | целый указатель способа оценки абсолютной погрешности компоненты приближенного решения и его производной: |
| MESTER=1 - | для оценки абсолютной погрешности используется первая формула (асимптотическая); |
| MESTER=2 - | для оценки абсолютной погрешности используется вторая формула (завышенная оценка). |
| (См. "Математическое описание"); | |
| H - | переменная с расширенной (Extended) точностью, содержащая начальное значение длины первого элементарного сегмента (аналог шага интегрирования для разностных методов). Может задаваться с учетом направления интегрирования, т.е. положительным, если XK > XN, отрицательным, если XK < XN, или без такого учета в виде абсолютной величины; |
| HMIN, HMAX - | соответственно минимальное и максимальное значение длины элементарного сегмента, которое разрешается использовать при разбиении интервала интегрирования на элементарные (частичные) сегменты (тип: с расширенной (Extended) точностью); |
| K2 - | порядок частичной суммы смещенного ряда Чебышёва второй производной решения, с помощью которой вычисляется частичная сумма порядка (K2+2) для оценивающего решения более высокого порядка точности на каждом элементарном сегменте разбиения области интегрирования. Этот параметр относится только ко второму приближенному решению более высокого порядка точности O(HK2 + 3) , K2 > K (т.е. к оценивающему решению) и его второй производной (см. "Математическое описание" и "Замечания по пользованию"; тип: целый); |
| IMAX2 - | целая переменная, задающая количество итераций, которое предполагается выполнить в дополнительном итерационном процессе для вычисления второго, оценивающего, решения более высокого порядка точности. Если правая часть дифференциального уравнения не зависит от переменных Y и Y', т.е. дифференциальное уравнение имеет вид Y'' = F(X), то число итераций при вычислении второго приближенного решения можно положить равным 1. В этом случае параметр IMAX2 = 1 (см. "Математическое описание" и "Замечания по использованию"); |
| NATTEM - | целая переменная, значение которой ограничивает число последовательных сокращений длины элементарного сегмента, совершаемых в одной точке, если приближенное решение (или его производная) на этом элементарном сегменте не достигает требуемой точности; |
| Y, DY - | искомое решение задачи Коши и его первая производная, вычисленные подпрограммой при значении аргумента XK; в качестве такого решения принимается значение второго, оценивающего, приближенного решения, поскольку оно имеет более высокий порядок точности по сравнению с первым приближенным решением. По такой же причине в качестве первой производной решения принимается первая производная оценивающего решения. Для системы уравнений (когда M ≠ 1) задаются массивами длины M. В случае совпадения значений параметров XN и XK значения Y и DY полагаются равными начальным значениям YN и DYN (тип: с расширенной (Extended) точностью); |
| RAB - | одномерный рабочий массив длины 2 * K * K + 5 * K + 7 * M * K + K2 * K2 + 8 * K2 + 6 * K2 * M + 24 * M + 5 (тип: с расширенной (Extended) точностью); |
| IERR - | целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом: |
|
IERR=65 - IEER=66 | когда какая-нибудь компонента решения или производной не может быть вычислена с требуемой точностью EPS; при этом IERR=65 указывает, что требуемая точность не может быть достигнута при минимальной длине частичного сегмента, равной HMIN, а IERR=66 показывает, что требуемая точность не может быть достигнута, т.к. исчерпано заданное число сокращений NATTEM длины элементарного сегмента в одной и той же точке интервала интегрирования. |
Версии: нет
Вызываемые подпрограммы
| DE98E - | выполнение одного шага приближенного интегрирования системы обыкновенных дифференциальных уравнений второго порядка методом рядов Чебышёва с контролем точности |
|
Кроме того, используются рабочие подпрограммы DE70EH, DE80EH, DE80EI, DE70EF, DE70EQ, DE71EE, DE70EP, DE71ET, DE71EP, DE71EI, DE71EF, DE71ES, DE70EA, DE70EC, DE75EK, DE75EE, DE75EN, DE98EK, DE98E1, DE98E0. UTDE75 - подпрограмма выдачи диагностических сообщений. |
Замечания по использованию
|
В ообщем случае заданная точность не гарантируется. Разбиение промежутка интегрирования на элементарные сегменты выполняется для того, чтобы на каждом таком сегменте ряды Чебышёва для решения и его первой и второй производных были быстросходящимися рядами; другими словами, чтобы убывание коэффициентов этих рядов Чебышёва на элементарном сегменте происходило достаточно быстро, вследствие чего можно было бы считать частичные суммы этих рядов близкими к многочленам наилучшего равномерного приближения на элементарном сегменте для решения и его производных. Если при вычислении первого приближенного решения начальное приближение для коэффициентов Чебышёва правой части системы (т.е. функции Φ(α)) определяется первым способом (т.е. при INIAPR = 1), то для получения максимального порядка точности первого приближенного решения в конце элементарного сегмента необходимо выполнить в итерационном процессе не менее K итераций; тогда IMAX ≥ K. Если начальное приближение коэффициентов Чебышёва функции Φ(α) определяется вторым способом (т.е. при INIAPR = 2), то для получения максимального порядка точности первого приближенного решения необходимо выполнить в итерационном процессе не менее K + 1 итераций; в этом случае IMAX ≥ K + 1. Однако в некоторых случаях при втором способе определения начального приближения итерационный процесс может сойтись за значительно меньшее число итераций. Если длина сегмента оказалась достаточно малой (или, вернее сказать, сегмент выбран довольно удачно), то хорошая точность первого приближенного решения может быть получена и с существенно меньшим числом итераций при любом способе выбора начального приближения. Для получения максимального порядка точности второго приближенного (оценивающего) решения необходимо выполнить в дополнительном итерационном процессе не менее (K2 - K) итераций. Указанное в данном разделе число итераций IMAX и IMAX2 носит асимптотический характер (когда все наши оценки справедливы при H --> 0). На практике же количество итераций в каждом из двух итерационных процессов зависит от длины сегмента, от заданной точности, от используемых порядков частичных сумм и от конкретной системы дифференциальных уравнений. Поэтому числа итераций в обоих итерационных процессах, т.е. значения параметров IMAX и IMAX2, могут быть как меньше, так и больше рекомендуемых здесь значений. При работе подпрограммы значения параметров M, XN, YN, DYN, XK, K, INIAPR, IMAX, IU, EPS, THRESH, MEXACT, NUMBES, MESTER, HMIN, HMAX, K2, IMAX2, NATTEM сохраняются. Если после работы подпрограммы нет необходимости иметь начальные значения решения и производной, то параметры YN и Y, а в случае производной параметры DYN и DY при обращении к ней можно совместить. Так как при интегрировании системы уравнений с помощью подпрограммы DE65E используются глобальные записи (структуры данных) с именами _COM70D, _COM80D, _COM75D и _COM98D, пользователь не должен портить элементы этих структур. Структуры определены в модуле Lstruct и являются объектами типа: type type type type
COM70D = record COM75D = record COM80D = record COM98D = record
elm1: Extended; // WC1 elm1: Integer; // LASN0 elm1: Integer; //JLAS elm1: Integer; //JLASA2
elm2: Extended; // WC2 elm2: Integer; // LASN2 elm2: Integer; //JLAS1 elm2: Integer; //JLAS12
elm3: Extended; // WC3 elm3: Extended; // WC30 elm3: Integer; //JLAS2 elm3: Integer; //JLAS22
elm4: Integer; // LASN elm4: Extended; // WC12 elm4: Integer; //JLAS3 elm4: Integer; //JLAS32
elm5: Extended; // HD2 elm5: Extended; // WC22 elm5: Extended; //H2D4 elm5: Integer; //JLAS0
elm6: Extended; // HD4 elm6: Extended; // WC32 elm6: Extended; //H2D8 elm6: Integer; //JLAS10
elm7: Extended; // HOLD end; elm7: Extended; //H2D16 elm7: Integer; //JLAS20
end; var elm8: Extended; //H2D32 elm8: Integer; //JLAS30
var _COM75D : COM75D; elm9: Extended; //H2D96 end;
_COM70D : COM70D; end; var
var _COM98D :COM98D;
_COM80D :COM80D;
При работе подпрограммы DE98E используется глобальная запись (структура данных) _STAT76 с элементами elm1 (NACCEP) и elm2 (NREJEC) целого типа и пользователь не должен портить значения элементов этой записи. для сбора полезной статистической информации. Целая переменная NACCEZ после выхода из подпрограммы содержит число принятых, т.е. удовлетворяющих условию достижения заданной точности, элементарных сегментов, на которые автоматически разбивался весь промежуток интегрирования [XN, XK], а целая переменная NREJEC содержит число отклоненных элементарных сегментов, а именно тех сегментов, на которых вычисленная оценка погрешности превышает заданное значение EPS, т.е. число таких сегментов, на которых не достигается заданная точность и длины которых подвергались сокращению. Когда приближенное решение задачи Коши не может быть вычислено в конце интервала интегрирования XK и работа подпрограммы DE65E прерывается в некоторой точке интервала интегрирования при значении параметра IERR = 65 или IERR = 66, то значение переменной NACCEZ содержит число элементарных сегментов, считая от начала промежутка интегрирования XN до этой точки прерывания, на которых приближенное решение вычислено с требуемой точностью. Значение переменной NREJEC содержит в этом случае число отклоненных сегментов, считая от начала промежутка интегрирования и включая число отклоненных сегментов с началом в этой точке прерывания. |
Данный пример не только иллюстрирует правила использования подпрограммы DE65E, но и показывает способность реализованного в подпрограмме метода рядов Чебышёва вычислять приближенное решение задачи Коши со стабильно высокой точностью. Вычисления на Фортране проводились с 19-20 значащими цифрами.
y'' = qy', y(0) = eq, y'(0) = qeq, q = 4, 0 ≤ x ≤ 7.
Точное решение задачи Коши y(x) = eq(1+x). Решение вычисляется на отрезке [0, XK] при XK = 7. Программа обработки результатов F2TRЕА выполняет следующие действия. 1) Печатает номер s и границы XI = xs-1, XE = xs каждого элементарного сегмента [xs-1, xs], s = 1,2,3, ... . 2) Печатает приближенное значение Y решения задачи Коши в конце каждого элементарного сегмента, т.е. Y(XE); точное значение решения YT в ХЕ и относительную погрешность (YT − Y)/Y приближенного решения в точке XE. Аналогичная информация выдается и для первой производной решения в конце каждого элементарного сегмента, т.е. приближенное значение DY, точное значение DYT и относительная погрешность. 3) Печатает срединную точку MIDPNT элементарного сегмента [xs-1, xs] = [XI, XE] и значение частичной суммы SOL с коэффициентами Чебышёва AY для решения задачи Коши, вычисленное в точке MIDPNT; точное значение решения SOLT в точке MIDPNT и относительную погрешность (SOLT − SOL)/SOL. Печатается также значение частичной суммы DER с коэффициентами Чебышёва ADY для первой производной решения, вычисленное в точке MIDPNT; точное значение первой производной DERT в точке MIDPNT и относительная погрешность (DERT − DER)/DER. 4) Печатает вычисленные на каждом отдельном элементарном сегменте [XI, XE] = [xs-1, xs] коэффициенты Чебышёва AY для решения, коэффициенты Чебышёва ADY для его первой производной и коэффициенты Чебышёва AD2Y для второй производной решения (для правой части дифференциального уравнения).
Приводятся вызывающая программа, подпрограмма F вычисления значений правой части уравнения, подпрограмма F2TRЕА обработки результатов. Далее представлены результаты работы подпрограммы F2TRЕА для каждого элементарного сегмента [xs-1, xs]; s = 1,2,3, ... ; значение параметра IERR, число N вычислений правой части дифференциального уравнения; число элементарных сегментов NACCER, на которые автоматически разбивался интервал интегрирования [0, XK], и число отклоненных элементарных сегментов NREJEC; приближенное значение решения Y задачи Коши в конце XK интервала интегрирования, точное значение YT решения в ХК и относительная погрешность (YT − Y)/Y; приближенное значение первой производной DY в конце ХК интервала интегрирования, точное значение производной DYT в ХК и относительная погрешность (DYT − DY)/DY.
Unit ftde65e_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure ftde65e(X :Extended; var Y :Array of Extended; var DY :Array of Extended;
var Z :Array of Extended; M :Integer);
implementation
procedure ftde65e(X :Extended; var Y :Array of Extended; var DY :Array of Extended;
var Z :Array of Extended; M :Integer);
begin
Z[0] := _BLOCK.elm2*DY[0];
_BLOCK.elm1 := _BLOCK.elm1+1;
end;
function TDE65E: String;
var
M,K0,INIAPR,IMAX,MESTER,K20,IMAX2,NATTEM,IERR :Integer;
XN,YN,DYN,XK,H,HMIN,HMAX,Y,YT,DY,DYT :Extended;
RАВ :Array [0..1867] of Extended;
IU :Array [0..1] of Integer;
EPS :Array [0..1] of Extended;
THRESH :Array [0..1] of Extended;
МЕХАСТ :Array [0..1] of Integer;
NUMBES :Array [0..1] of Integer;
const
K :Integer = 18;
K2 :Integer = 25;
begin
Result := ''; { результат функции }
IU[0] := 2;
IU[1] := 2;
EPS[0] := 0.5e-12;
EPS[1] := 0.5e-12;
THRESH[0] := 1.e0;
THRESH[1] := 1.e0;
MEXACT[0] := 1;
MEXACT[1] := 1;
NUMBES[0] := 1;
NUMBES[1] := 1;
M := 1;
_BLOCK.elm2 := 4.0e0;
XN := 0.0;
YN := Exp(_BLOCK.elm2);
DYN := _BLOCK.elm2*YN;
ХК := 7.e0;
YT := Exp(_BLOCK.elm2*(XK+1.e0));
DYT := _BLOCK.elm2*YT;
K0 := K;
INIAPR := 1;
IМАХ := 28;
MESTER := 1;
H := 1.e0;
HMIN := 1.e-3;
НМАХ := XK;
K20 := K2;
IMAX2 := 3;
NАТТЕМ := 3;
_BLOCK.elm1 := 0;
DE65E(ftde65e,ftreat65e,M,XN,YN,DYN,XK,K0,INIAPR,IMAX,IU,EPS,THRESH,
MEXACT,M,MESTER,H,HMIN,HMAX,K20,IMAX2,NATTEM,Y,DY,RAB,
IERR);
// Операторы вывода на печать: IERR, Y,YT, DY,DYT, NACCEP,NREJEC
. . . . . . . . . . . . . . . . . . . . . . . . . .
UtRes('TDE65E',Result); { вывод результатов в файл TDE65E.txt }
exit;
end;
function ftreat65e(S :Integer; XI :Extended; XE :Extended; Y :Extended;
DY :Extended; var AY :Array of Extended;
var ADY :Array of Extended; var AD2Y :Array of Extended;
var M :Integer; var KP1 :Integer; var KP2 :Integer;
KP3 :Integer): String;
var
L,J :Integer;
YT,DYT,MIDPNT,ALPN,SOL,SOLT,DER,DERT,Q :Extended;
BK :Array [0..1] of Extended;
label
_10;
begin
Result := ''; { результат функции }
{ }
// Операторы вывода на печать: номер сегмента [XI,XE]
Q := 4.e0;
YT := Exp(Q*(XE+1.e0));
DYT := Q*YT;
. . . . . . . . . . . . . . . . . . . . . . . . . .
MIDPNT := 0.5e0*(XI+XE);
ALPN := 2.e0*(MIDPNT-XI)/(XE-XI)-1.e0;
DE70EC(M,KP2,ALPN,AY,SOL,BK);
SOLT := Exp(Q*(MIDPNT+1.e0));
. . . . . . . . . . . . . . . . . . . . . . . . . .
DE70EC(M,KP1,ALPN,ADY,DER,BK);
DERT := Q*SOLT;
. . . . . . . . . . . . . . . . . . . . . . . . . .
// Операторы вывода на печать: SOL,SOLT, DER,DERT, Chebyshev coefficients
. . . . . . . . . . . . . . . . . . . . . . . . . .
UtResAdd('ftreat65e',Result); { вывод результатов в файл ftreat65e.txt }
exit;
end;
Результаты:
----------------------------------------------------------
Coefficients AY, ADY and AD2Y on 1 segment:
XI= 0,00000000000000000E+000 XE= 1,00000000000000000E+000
Solution Y and YT at the end of the segment [XI,XE],i.e. Y(XE), YT(XE):
Y = 2,98095798704172828E+003
YT= 2,98095798704172828E+003
(YT-Y)/Y= -4,46926000078356152E-019
Derivative DY and DYT at the end of the segment [XI,XE], i.e. DY(XE), DYT(XE):
DY = 1,19238319481669131E+004
DYT= 1,19238319481669131E+004
(DYT-DY)/DY= -5,21413666758082177E-019
----------------------------------------------------------------
MIDPNT= 5,00000000000000000E-001
Solution SOL and SOLT at the midpoint 0.5*(XI+XE) of segment [XI, XE]:
SOL= 4,03428793492735123E+002
SOLT= 4,03428793492735123E+002
(SOLT-SOL)/SOL= 3,43995967359339353E-019
Derivative DER and DERT at the midpoint 0.5*(XI+XE) of segment [XI, XE]:
DER = 1,61371517397094049E+003
DERT = 1,61371517397094049E+003
(DERT-DER)/DER= 4,81594354303075094E-019
Number of Chebyshev coefficients for 1 component
------------------------------------------------------------------
coefficient for Y for Y' for Y''
-------------------------------------------------------------------------------
0 1,83930069637042288E+003 7,35720278548169154E+003 2,94288111419267661E+004
1 1,28341741430283353E+003 5,13366965721133410E+003 2,05346786288453364E+004
2 5,55883282067589359E+002 2,22353312827035743E+003 8,89413251308142974E+003
3 1,71650850167654808E+002 6,86603400670619233E+002 2,74641360268247693E+003
4 4,09307315646249336E+001 1,63722926258499734E+002 6,54891705033998937E+002
5 7,92792390915507391E+000 3,17116956366202957E+001 1,26846782546481182E+002
6 1,29111201884956405E+000 5,16444807539825621E+000 2,06577923015930236E+001
7 1,81251796057689625E-001 7,25007184230758464E-001 2,90002873692303327E+000
8 2,23494464457366862E-002 8,93977857829467032E-002 3,57591143131786599E-001
9 2,45622449179612735E-003 9,82489796718450664E-003 3,92995918687387702E-002
10 2,43426019571520271E-004 9,73704078286118536E-004 3,89481631314436000E-003
11 2,19642960809140037E-005 8,78571843236958029E-005 3,51428737294028765E-004
12 1,81876268148271868E-006 7,27505072590237459E-006 2,91002029017446706E-005
13 1,39143903152865316E-007 5,56575612525306509E-007 2,22630245071478439E-006
14 9,89194048787709822E-009 3,95677619533781452E-008 1,58271050428732130E-007
15 6,56736288586542544E-010 2,62694520418900919E-009 1,05077813256082564E-008
16 4,08961537934832197E-011 1,63584638185592518E-010 6,54338177391577958E-010
17 2,39784329047044730E-012 9,59136140608313530E-012 3,83644817303352426E-011
18 1,32825871686546612E-013 5,31294433602101811E-013 2,12560177792475713E-012
19 6,97276526742657699E-015 2,78986446517792112E-014
20 3,48271497793810597E-016
****************************************************************
Coefficients AY, ADY and AD2Y on 2 segment:
XI= 1,00000000000000000E+000 XE= 2,34442743069806723E+000
Solution Y and YT at the end of the segment [XI,XE],i.e. Y(XE), YT(XE):
Y = 6,45454321855636837E+005
YT= 6,45454321855636816E+005
(YT-Y)/Y= -3,24087660930651905E-017
Derivative DY and DYT at the end of the segment [XI,XE], i.e. DY(XE), DYT(XE):
DY = 2,58181728742254733E+006
DYT= 2,58181728742254726E+006
(DYT-DY)/DY= -2,48349783647945212E-017
----------------------------------------------------------------
MIDPNT= 1,67221371534903362E+000
Solution SOL and SOLT at the midpoint 0.5*(XI+XE) of segment [XI, XE]:
SOL= 4,38642475828112541E+004
SOLT= 4,38642475828112538E+004
(SOLT-SOL)/SOL= -6,96543069159061782E-018
Derivative DER and DERT at the midpoint 0.5*(XI+XE) of segment [XI, XE]:
DER = 1,75456990331244960E+005
DERT = 1,75456990331245015E+005
(DERT-DER)/DER= 3,17494050128316636E-016
Number of Chebyshev coefficients for 1 component
------------------------------------------------------------------
coefficient for Y for Y' for Y''
-------------------------------------------------------------------------------
0 3,34088067307084037E+005 1,33635226922833615E+006 5,34540907691334445E+006
1 2,61947790871362022E+005 1,04779116348544809E+006 4,19116465394179218E+006
2 1,39248401818147358E+005 5,56993607272589429E+005 2,22797442909035754E+006
3 5,47987865612589608E+004 2,19195146245035838E+005 8,76780584980143188E+005
4 1,69685703225809977E+004 6,78742812903239857E+004 2,71497125161295807E+005
5 4,31314506026576912E+003 1,72525802410630724E+004 6,90103209642521852E+004
6 9,27745203346951743E+002 3,71098081338780401E+003 1,48439232535511403E+004
7 1,72742169801512859E+002 6,90968679206049397E+002 2,76387471682414592E+003
8 2,83324414526634169E+001 1,13329765810652340E+002 4,53319063242576109E+002
9 4,15030203248103421E+000 1,66012081299233262E+001 6,64048325196732396E+001
10 5,49076252404938374E-001 2,19630500961929038E+000 8,78522003846561772E+000
11 6,62121078961101891E-002 2,64848431584191898E-001 1,05939372633048146E+000
12 7,33396849516432157E-003 2,93358739805164936E-002 1,17343495918666774E-001
13 7,51064828628876670E-004 3,00425931443792580E-003 1,20170372564922800E-002
14 7,15141973187764739E-005 2,86056789251018319E-004 1,14422715660555951E-003
15 6,36211011722251338E-006 2,54484404496791575E-005 1,01793761542623429E-004
16 5,31078718491696946E-007 2,12431486290884217E-006 8,49725991047822760E-006
17 4,17544656827060535E-008 1,67017884669105734E-007 6,68071808984649351E-007
18 3,10234024855234579E-009 1,24093769059571949E-008 4,96368568470018090E-008
19 2,18485360819894547E-010 8,73924202761155552E-010
20 1,46239204556943689E-011
****************************************************************
Coefficients AY, ADY and AD2Y on 3 segment:
XI= 2,34442743069806723E+000 XE= 3,73719668264729049E+000
Solution Y and YT at the end of the segment [XI,XE],i.e. Y(XE), YT(XE):
Y = 1,69571755347846709E+008
YT= 1,69571755347846691E+008
(YT-Y)/Y= -1,04780943367533229E-016
Derivative DY and DYT at the end of the segment [XI,XE], i.e. DY(XE), DYT(XE):
DY = 6,78287021391386816E+008
DYT= 6,78287021391386766E+008
(DYT-DY)/DY= -7,42305618451402507E-017
----------------------------------------------------------------
MIDPNT= 3,04081205667267886E+000
Solution SOL and SOLT at the midpoint 0.5*(XI+XE) of segment [XI, XE]:
SOL= 1,04618747055159456E+007
SOLT= 1,04618747055159450E+007
(SOLT-SOL)/SOL= -5,06826382516371075E-017
Derivative DER and DERT at the midpoint 0.5*(XI+XE) of segment [XI, XE]:
DER = 4,18474988220637548E+007
DERT = 4,18474988220637802E+007
(DERT-DER)/DER= 6,07583119623828434E-016
Number of Chebyshev coefficients for 1 component
------------------------------------------------------------------
coefficient for Y for Y' for Y''
-------------------------------------------------------------------------------
0 8,59938816313994066E+007 3,43975526525597627E+008 1,37590210610239037E+009
1 6,81750402542297978E+007 2,72700161016919190E+008 1,09080064406767662E+009
2 3,70446101287317982E+007 1,48178440514927190E+008 5,92713762059708624E+008
3 1,49794228511506729E+007 5,99176914046026875E+007 2,39670765618410630E+008
4 4,77920184910505155E+006 1,91168073964202023E+007 7,64672295856807109E+007
5 1,25366937858825908E+006 5,01467751435303329E+006 2,00587100574120578E+007
6 2,78566238311266149E+005 1,11426495324506238E+006 4,45705981298019503E+006
7 5,36160118669407520E+004 2,14464047467761487E+005 8,57856189871008773E+005
8 9,09440545501596569E+003 3,63776218200628728E+004 1,45510487280227573E+005
9 1,37818170354157434E+003 5,51272681416568972E+003 2,20509072566483305E+004
10 1,88669408276619315E+002 7,54677633106128285E+002 3,01871053241630130E+003
11 2,35466266070045432E+001 9,41865064278301825E+001 3,76746025706914878E+002
12 2,69969914164206526E+000 1,07987965664689199E+001 4,31951862636027784E+001
13 2,86212396672665026E-001 1,14484958664118392E+000 4,57939834558254688E+000
14 2,82148655744760381E-002 1,12859462279332026E-001 4,51437848691654153E-001
15 2,59892839593161895E-003 1,03957135682723382E-002 4,15828540376423916E-002
16 2,24639707703781221E-004 8,98558821950022629E-004 3,59423552788484812E-003
17 1,82889053182682798E-005 7,31556334808943159E-005 2,92622705842404685E-004
18 1,40718288925627613E-006 5,62874045825868883E-006 2,25146058880909550E-005
19 1,02630463032810196E-007 4,10512288237743903E-007
20 7,11412247152372234E-009
****************************************************************
Coefficients AY, ADY and AD2Y on 4 segment:
XI= 3,73719668264729049E+000 XE= 5,11593009771574037E+000
Solution Y and YT at the end of the segment [XI,XE],i.e. Y(XE), YT(XE):
Y = 4,21171303210922011E+010
YT= 4,21171303210921953E+010
(YT-Y)/Y= -1,36833256391983252E-016
Derivative DY and DYT at the end of the segment [XI,XE], i.e. DY(XE), DYT(XE):
DY = 1,68468521284368801E+011
DYT= 1,68468521284368781E+011
(DYT-DY)/DY= -1,17551000481542175E-016
----------------------------------------------------------------
MIDPNT= 4,42656339018151543E+000
Solution SOL and SOLT at the midpoint 0.5*(XI+XE) of segment [XI, XE]:
SOL= 2,67242880518108898E+009
SOLT= 2,67242880518108874E+009
(SOLT-SOL)/SOL= -9,29604920801041781E-017
Derivative DER and DERT at the midpoint 0.5*(XI+XE) of segment [XI, XE]:
DER = 1,06897152207243502E+010
DERT = 1,06897152207243549E+010
(DERT-DER)/DER= 4,45286855690171237E-016
Number of Chebyshev coefficients for 1 component
------------------------------------------------------------------
coefficient for Y for Y' for Y''
-------------------------------------------------------------------------------
0 2,14839416590055952E+010 8,59357666360223810E+010 3,43743066544089500E+011
1 1,69792776049850162E+010 6,79171104199400648E+010 2,71668441679760234E+011
2 9,16881429688095494E+009 3,66752571875238193E+010 1,46701028750095253E+011
3 3,67893368545253526E+009 1,47157347418101403E+010 5,88629389672405391E+010
4 1,16378523489244691E+009 4,65514093956978692E+009 1,86205637582791294E+010
5 3,02544248094426245E+008 1,21017699237770442E+009 4,84070796951080367E+009
6 6,66034128380148977E+007 2,66413651352059183E+008 1,06565460540822659E+009
7 1,26981671475525877E+007 5,07926685902100701E+007 2,03170674360833359E+008
8 2,13324836575162405E+006 8,53299346300631378E+006 3,41319738520208071E+007
9 3,20149222190434394E+005 1,28059688876162583E+006 5,12238755504381592E+006
10 4,34005618011406717E+004 1,73602247204498425E+005 6,94408988816463650E+005
11 5,36348241790096946E+003 2,14539296715694797E+004 8,58157186854548927E+004
12 6,08890871777299235E+002 2,43556348709079537E+003 9,74225394793089072E+003
13 6,39151618383266000E+001 2,55660647343832521E+002 1,02264258919294662E+003
14 6,23839754088950007E+000 2,49535901601568749E+001 9,98143605656514410E+001
15 5,68930710370103687E-001 2,27572283870776793E+000 9,10289131059107604E+000
16 4,86870697545423722E-002 1,94748277361511360E-001 7,78993155130592640E-001
17 3,92437242791693435E-003 1,56974920848217155E-002 6,27900010003941134E-002
18 2,98939154597905428E-004 1,19575866441219020E-003 4,78296329674776644E-003
19 2,15850782706080778E-005 8,63384739689538218E-005
20 1,48129263535869606E-006
****************************************************************
Coefficients AY, ADY and AD2Y on 5 segment:
XI= 5,11593009771574037E+000 XE= 6,49863199748470672E+000
Solution Y and YT at the end of the segment [XI,XE],i.e. Y(XE), YT(XE):
Y = 1,06281577856571496E+013
YT= 1,06281577856571477E+013
(YT-Y)/Y= -1,81346177985619704E-016
Derivative DY and DYT at the end of the segment [XI,XE], i.e. DY(XE), DYT(XE):
DY = 4,25126311426285976E+013
DYT= 4,25126311426285909E+013
(DYT-DY)/DY= -1,59272373045262237E-016
----------------------------------------------------------------
MIDPNT= 5,80728104760022354E+000
Solution SOL and SOLT at the midpoint 0.5*(XI+XE) of segment [XI, XE]:
SOL= 6,69049704081582312E+011
SOLT= 6,69049704081582223E+011
(SOLT-SOL)/SOL= -1,33454595867585380E-016
Derivative DER and DERT at the midpoint 0.5*(XI+XE) of segment [XI, XE]:
DER = 2,67619881632632772E+012
DERT = 2,67619881632632889E+012
(DERT-DER)/DER= 4,36533724800513174E-016
Number of Chebyshev coefficients for 1 component
------------------------------------------------------------------
coefficient for Y for Y' for Y''
-------------------------------------------------------------------------------
0 5,41242196770310655E+012 2,16496878708124263E+013 8,65987514832496984E+013
1 4,28136968067799080E+012 1,71254787227119632E+013 6,85019148908478460E+013
2 2,31604256633441053E+012 9,26417026533764201E+012 3,70566810613505612E+013
3 9,31345258607847281E+011 3,72538103443138892E+012 1,49015241377255495E+013
4 2,95335300154771086E+011 1,18134120061908415E+012 4,72536480247633158E+012
5 7,69736830596593793E+010 3,07894732238637363E+011 1,23157892895454562E+012
6 1,69901157281032472E+010 6,79604629124128772E+010 2,71841851649648734E+011
7 3,24796214638614749E+009 1,29918485855445132E+010 5,19673943421761569E+010
8 5,47139097708847188E+008 2,18855639083533856E+009 8,75422556334013461E+009
9 8,23392580382260345E+007 3,29357032152873390E+008 1,31742812861075943E+009
10 1,11932640711939065E+007 4,47730562847580317E+007 1,79092225138613053E+008
11 1,38714483823036988E+006 5,54857935291201581E+006 2,21943174114229251E+007
12 1,57918547268638376E+005 6,31674189069568679E+005 2,52669675616035238E+006
13 1,66234210379297255E+004 6,64936841491477290E+004 2,65974736544499174E+005
14 1,62710546866257221E+003 6,50842187363146563E+003 2,60336874732989818E+004
15 1,48809579245211419E+002 5,95238316153066352E+002 2,38095325378328562E+003
16 1,27707452924524734E+001 5,10829806968802796E+001 2,04331937480717897E+002
17 1,03229967760659111E+000 4,12919939151102293E+000 1,65168086271733046E+001
18 7,88593047537755331E-002 3,15437751450605697E-001 1,26173220388591290E+000
19 5,71030875356647659E-003 2,28407156311782255E-002
20 3,92991354428195271E-004
****************************************************************
Coefficients AY, ADY and AD2Y on 6 segment:
XI= 6,49863199748470672E+000 XE= 7,00000000000000000E+000
Solution Y and YT at the end of the segment [XI,XE],i.e. Y(XE), YT(XE):
Y = 7,89629601826807080E+013
YT= 7,89629601826806952E+013
(YT-Y)/Y= -1,62224837928231183E-016
Derivative DY and DYT at the end of the segment [XI,XE], i.e. DY(XE), DYT(XE):
DY = 3,15851840730722831E+014
DYT= 3,15851840730722781E+014
(DYT-DY)/DY= -1,59229620551354967E-016
----------------------------------------------------------------
MIDPNT= 6,74931599874235336E+000
Solution SOL and SOLT at the midpoint 0.5*(XI+XE) of segment [XI, XE]:
SOL= 2,89694804931689009E+013
SOLT= 2,89694804931688960E+013
(SOLT-SOL)/SOL= -1,67101748035130061E-016
Derivative DER and DERT at the midpoint 0.5*(XI+XE) of segment [XI, XE]:
DER = 1,15877921972675602E+014
DERT = 1,15877921972675584E+014
(DERT-DER)/DER= -1,58740076640149164E-016
Number of Chebyshev coefficients for 1 component
------------------------------------------------------------------
coefficient for Y for Y' for Y''
-------------------------------------------------------------------------------
0 7,34442833286130009E+013 2,93777133314452002E+014 1,17510853325780801E+015
1 3,28559331508849087E+013 1,31423732603539635E+014 5,25694930414158539E+014
2 7,91171446492150346E+012 3,16468578596860139E+013 1,26587431438744055E+014
3 1,29542500026685414E+012 5,18170000106741657E+012 2,07268000042696663E+013
4 1,60372174819895004E+011 6,41488699279580014E+011 2,56595479711832006E+012
5 1,59482565151173812E+010 6,37930260604695259E+010 2,55172104241878115E+011
6 1,32476418432188610E+009 5,29905673728754535E+009 2,11962269491501465E+010
7 9,44623638066185708E+007 3,77849455226473742E+008 1,51139782090585995E+009
8 5,89951140897937951E+006 2,35980456359158618E+007 9,43921825436587334E+007
9 3,27734083423462133E+005 1,31093633369371419E+006 5,24374533484554291E+006
10 1,63941457838934928E+004 6,55765831371135823E+004 2,62306332553386688E+005
11 7,45809308270786925E+002 2,98323723392003624E+003 1,19329488835334778E+004
12 3,11103498284058597E+001 1,24441398344270662E+002 4,97765524864196777E+002
13 1,19817279702010167E+000 4,79268964952094540E+000 1,91707992553710938E+001
14 4,28576974439898209E-002 1,71430710776158753E-001 6,85813903808593750E-001
15 1,43099136210932421E-003 5,72469598819265691E-003 2,29482650756835938E-002
16 4,48032295474605952E-005 1,80289337585852407E-004 7,24792480468750000E-004
17 1,33293167407931954E-006 5,53027228717633088E-006 -6,58035278320312500E-005
18 3,34497592584242188E-008 -4,94743932451001135E-007 -2,52723693847656250E-005
19 1,21132968716658737E-009 7,26649689555373620E-007
20 1,06719678768094280E-008
****************************************************************
IERR= 0
N= 4661
NACCEP= 6 NREJEC= 1
Y = 7,89629601826807080E+013
YT= 7,89629601826806952E+013
(YT-Y)/Y= -1,62224837928231183E-016
DY = 3,15851840730722831E+014
DYT= 3,15851840730722781E+014
(DYT-DY)/DY= -1,59229620551354967E-016
Из приведенных результатов наглядно видно быстрое убывание на каждом частичном сегменте коэффициентов Чебышёва решения и коэффициентов Чебышёва его первой и второй производных.