|
Текст подпрограммы и версий (Паскаль) de88e_p.zip |
Тексты тестовых примеров (Паскаль) tde88e1_p.zip , tde88e2_p.zip , tde88e3_p.zip , tde88e4_p.zip , tde88e5_p.zip tde88e6_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] в равномерно сходящиеся ряды по смещенным многочленам Чебышёва первого рода
∞ 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
(8) Y'' (xs + αH) ≈ ∑ ' ai*[Φ] Ti*(α) .
i=0
В этом случае аналитическое решение задачи (1),(2) состоит из совокупности частичных сумм рядов Чебышёва, построенных на этих элементарных сегментах. Порядок этих частичных сумм и длина элементарных сегментов задаются пользователем при обращении к подпрограмме. Число NX элементарных сегментов, на которые разбивается промежуток интегрирования, равно [(XK-XN) / H], если длина промежутка интегрирования является целым кратным H. Если длина интервала интегрирования не является целым кратным H, то число элементарных сегментов NX равно [(XK-XN) / H] + 1; в этом случае последний элементарный сегмент считается нестандартным (квадратные скобки означают целую часть числа).
При разбиении промежутка интегрирования на элементарные сегменты решение задачи сводится к определению нескольких наборов коэффициентов 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 DE88E(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; IORDY :Integer; var H :Extended;
var Y :Array of Extended; var DY :Array of Extended;
var NX :Integer; var XS :Array of Extended;
var AY :Array of Extended; var ADY :Array of Extended;
var AD2Y :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 - | конец интервала интегрирования. На отрезке [XN, XK] вычисляется приближенное аналитическое решение задачи Коши в виде одной частичной суммы ряда Чебышёва либо в виде совокупности частичных сумм. XK может быть больше, меньше или равно XN (тип: с расширенной (Extended) точностью); |
| K - | порядок частичной суммы смещенного ряда Чебышёва, с помощью которой аппроксимируется вторая производная решения задачи Коши на каждом элементарном сегменте разбиения интервала интегрирования; при этом само решение задачи Коши приближается на каждом элементарном сегменте частичной суммой (K + 2) - го порядка, а его первая производная - частичной суммой (K + 1) - го порядка; K≥2 (см. "Математическое описание" и "Замечания по использованию"; тип: целый); |
| INIAPR - |
целый указатель способа выбора начального
приближения коэффициентов Чебышёва для второй производной
решения на каждом элементарном сегменте: |
| INIAPR=1 - |
для первого способа, когда начальное
приближение определяется только с
использованием значения решения и его первой производной в
начале каждого элементарного сегмента; |
| INIAPR=2 - | для второго способа, когда начальное приближение коэффициентов Чебышёва на текущем элементарном сегменте (начиная со второго) определяется через коэффициенты Чебышёва, вычисленные на предыдущем элементарном сегменте, т.е. путем экстраполяции коэффициентов с предыдущего сегмента на следующий (см. "Математическое описание"); |
| IMAX - | целая переменная, задающая количество итераций, которое предполагается выполнить в итерационном процессе вычисления коэффициентов Чебышёва для второй производной решения задачи Коши на каждом элементарном сегменте, исходя из некоторого начального приближения, способ определения которого задается параметром INIAPR; IMAX≥1. Для получения максимального порядка точности приближенного решения необходимо выполнить не менее K итераций (см. "Математическое описание" и "Замечания по использованию"); |
| IORDY - | целая переменная, принимающая следующие значения: |
| IORDY=0 - |
если требуется получить приближенное аналитическое представление в виде
частичной суммы ряда Чебышёва только для решения задачи Коши; в этом
случае вычисляются коэффициенты разложения только для решения задачи
Коши; |
| IORDY=1 - | если требуется получить аналитическое представление в виде частичной суммы как для решения задачи Коши, так и для его первой производной; в этом случае вычисляются коэффициенты разложения для решения и его первой производной; |
| IORDY=2 - | если требуется получить аналитическое представление в виде частичной суммы как для решения задачи Коши, так и для его первой и второй производных; в этом случае вычисляются коэффициенты разложения для решения и его первой и второй производных; |
| для решения Y(X) на каждом элементарном сегменте строится частичная сумма (K+2)-го порядка, для его производной Y'(X) - частичная сумма (K+1)-го порядка, для второй производной Y''(X) - частичная сумма K-го порядка; |
| H - | переменная с расширенной (Extended) точностью, содержащая значение длины элементарных сегментов, на которые разбивается интервал интегрирования (диаметр разбиения промежутка интегрирования или аналог шага интегрирования для разностных методов). Предполагается, что на каждом элементарном сегменте ряды Чебышёва для решения и его производных являются быстросходящимися рядами. Может задаваться с учетом направления интегрирования, т.е. положительным, если XK > XN, отрицательным, если XK < XN , или без такого учета в виде абсолютной величины; |
| Y, DY - | на выходе из подпрограммы содержат значения решения задачи Коши и его первой производной соответственно, вычисленные подпрограммой при значении аргумента XK. Для системы уравнений (когда M ≠ 1) задаются одномерными массивами длины M. В случае совпадения значений параметров XN и XK значения Y и DY полагаются равными начальным значениям YN и DYN (тип: с расширенной (Extended) точностью); |
| NX - | целая переменная; на выходе из подпрограммы содержит число элементарных сегментов, на которые разбивается интервал интегрирования [XN, XK] (см. "Математическое описание"); |
| XS - | одномерный массив длиной NX + 1. На выходе из подпрограммы содержит точки x0 = XN, xs = x0 + s*H, s = 1, ... , NX - 1, xNX = XK, задающие разбиение интервала [XN, XK] интегрирования и являющиеся границами элементарных сегментов. При этом XS(S + 1) = x0 + S*H, S=0,1, ... , NX - 1, XS(NX + 1) = XK. Если XN < XK, то H > 0 и в массиве XS точки деления располагаются в порядке возрастания. Если XK < XN, то при работе подпрограммы значение H будет выбрано отрицательным; в этом случае нумерация точек разбиения интервала интегрирования ведется справа налево и в массиве XS точки деления будут расположены в порядке убывания: XS(1) = XN, XS(S+1) = x0 - S*|H|, S=1, ... , NX - 1, XS(NX + 1) = XK (тип: с расширенной (Extended) точностью); |
| AY - |
двумерный массив с измерениями M, (K + 3) * NX. На выходе из
подпрограммы содержит коэффициенты Чебышёва для решения Y(X). При этом переменная с
индексом AY (N, (S - 1) * (K + 3) + I + 1)
содержит I-й коэффициент Чебышёва для N-й компоненты решения yN(X)
на S-м элементарном сегменте, I = 0, 1, ... , K + 2;
S = 1, 2, ... , NX. Массив представляет последовательность групп
столбцов по K + 3 столбца в каждой группе; число групп равно NX, а длина всех
столбцов одинакова и равна M. Номер группы - это номер элементарного сегмента, на которые
разбивается интервал интегрирования, а номер столбца в каждой группе (уменьшенный на
единицу) - это номер коэффициента Чебышёва для решения Y(X) на элементарном сегменте,
соответствующем этой группе. Заметим, что удобно передавать подпрограмме DE88E в качестве
фактического параметра, соответствующего данному формальному параметру AY,
массив, описываемый в вызывающей программной единице как трехмерный массив BY с
измерениями M, K + 3, NX. При этом переменная с индексом BY (N, I + 1, S)
будет содержать I-й коэффициент Чебышёва N-й компоненты решения yN(X)
на S-м элементарном сегменте [xS-1, xS]. Другими словами, индексы этого
массива (индексные выражения N, I+ 1, S) имеют следующие значения: N - номер компоненты решения, N = 1, 2, ... , M; I - номер коэффициента Чебышёва, I = 0, 1, ... , K + 2; S - номер элементарного сегмента, S = 1, 2, ... , NX. Если интегрируется одно скалярное уравнение (т.е. M = 1), то первое измерение массива BY равно 1 и в описании массива его можно опустить. В этом случае переменная с индексом, служащая для обращения к элементам массива BY, имеет вид BY (I + 1, S) и представляет I-й коэффициент Чебышёва решения Y(X) на S-м элементарном сегменте [xS-1, xS]. Если интервал интегрирования [XN, XK] не разбивается на элементарные сегменты (т.е. H = XK - XN), то третье измерение массива BY равно 1 и в описании массива его можно опустить. В этом случае переменная с индексом, служащая для обращения к элементам массива BY, имеет вид BY (N, I + 1) и представляет I-й коэффициент Чебышёва для N-й компоненты решения yN(X) на сегменте [XN, XK]. Если интегрируется одно скалярное уравнение и интервал интегрирования [XN, XK] не разбивается на элементарные сегменты (т.е. M = 1 и H = XK - XN), то одновременно первое и третье измерения массива BY равны 1 и эти измерения в описании массива можно опустить. В этом случае переменная с индексом, служащая для обращения к элементам массива BY, имеет вид BY (I + 1) и представляет I-й коэффициент Чебышёва решения Y(X) на [XN, XK] (см. "Математическое описание", "Замечания по использованию" и "Примеры использования"; тип: с расширенной (Extended) точностью); |
| ADY - |
двумерный массив с измерениями M, (K + 2) * NX. На выходе из подпрограммы содержит коэффициенты
Чебышёва для первой производной решения Y'(X). При этом переменная с индексом
ADY (N, (S - 1) * (K + 2) + I + 1)
содержит I-й коэффициент Чебышёва для N-й компоненты первой производной решения y'N(X) на S-м
элементарном сегменте, I = 0, 1, ... , K + 1; S = 1, 2, ... , NX.
Массив представляет последовательность групп столбцов по K + 2 столбца в каждой группе; число
групп равно NX, а длина всех столбцов одинакова и равна M. Номер группы - это номер элементарного сегмента, на
которые разбивается интервал интегрирования, а номер столбца в каждой группе (уменьшенный на единицу) - это
номер коэффициента Чебышёва для первой производной решения Y'(X) на элементарном сегменте, соответствующем этой
группе. Как и для формального параметра AY, удобно передавать подпрограмме DE88E в качестве фактического
параметра, соответствующего данному формальному параметру ADY, массив, описываемый в вызывающей
программной единице, как трехмерный BDY с измерениями M, K + 2, NX. При этом
переменная с индексом BDY (N, I + 1, S) будет содержать I-й коэффициент Чебышёва
N-й компоненты первой производной решения y'N(X) на S-м элементарном сегменте [xS-1, xS].
Другими словами, индексы этого массива (индексные выражения N, I + 1, S) имеют следующее значение: N - номер компоненты производной решения, N = 1, 2, ... , M; I - номер коэффициента Чебышёва, I = 0, 1, ... , K + 1; S - номер элементарного сегмента, S = 1, 2, ... , NX. Если при обращении к подпрограмме DE88E параметр IORDY принимает нулевое значение (IORDY = 0), т.е. когда требуется получить аналитическое представление для решения Y(X) и не требуется аналитическое решение для производной Y'(X), то массив ADY, представляющий собой массив переменной длины, является рабочим массивом и может иметь размер M * (K + 2). Обращение с массивом BDY при интегрировании скалярных уравнений (M = 1) или при интегрировании уравнений без разбиения интервала интегрирования на элементарные сегменты (H = XK - XN) полностью аналогично обращению с массивом BY (соответствующим формальному параметру AY) и изложено в описании параметра AY (см. "Математическое описание", "Замечания по использованию" и "Примеры использования"; тип: с расширенной (Extended) точностью); |
| AD2Y - |
двумерный массив с измерениями M, (K + 1) * NX. На выходе из подпрограммы содержит коэффициенты
Чебышёва для второй производной решения Y''(X). При этом переменная с индексом
AD2Y (N, (S - 1) * (K + 1) + I + 1)
содержит I-й коэффициент Чебышёва для N-й компоненты второй производной решения y''N(X) на S-м
элементарном сегменте, I = 0, 1, ... , K; S = 1, 2, ... , NX.
Массив представляет последовательность групп столбцов по K + 1 столбцу в каждой группе; число
групп равно NX, а длина всех столбцов одинакова и равна M. Номер группы - это номер элементарного сегмента, на
которые разбивается интервал интегрирования, а номер столбца в каждой группе (уменьшенный на единицу) - это
номер коэффициента Чебышёва для второй производной решения Y''(X) на элементарном сегменте, соответствующем этой
группе. Как и для формальных параметров AY, ADY, удобно передавать подпрограмме DE88E в качестве фактического
параметра, соответствующего данному формальному параметру AD2Y, массив, описываемый в вызывающей
программной единице, как трехмерный массив BD2Y с измерениями M, K + 1, NX. При этом
переменная с индексом BD2Y (N, I + 1, S) будет содержать I-й коэффициент Чебышёва
N-й компоненты второй производной решения y''N(X) на S-м элементарном сегменте [xS-1, xS].
Другими словами, индексы этого массива (индексные выражения N, I + 1, S) имеют следующее значение: N - номер компоненты решения, N = 1, 2, ... , M; I - номер коэффициента Чебышёва, I = 0, 1, ... , K; S - номер элементарного сегмента, S = 1, 2, ... , NX. Если при обращении к подпрограмме DE88E параметр IORDY принимает нулевое значение или единичное значение (IORDY = 0 или IORDY = 1), т.е. когда не требуется получать аналитическое представление для второй производной решения Y''(X), то массив AD2Y, представляющий собой массив переменной длины, является рабочим массивом и может иметь размер M * (K + 1). Обращение с массивом BD2Y при интегрировании скалярных уравнений (M = 1) или при интегрировании уравнений без разбиения интервала интегрирования на элементарные сегменты (H = XK - XN) полностью аналогично обращению с массивом BY (соответствующим формальному параметру AY) и изложено в описании параметра AY (см. "Математическое описание", "Замечания по использованию" и "Примеры использования"; тип: с расширенной (Extended) точностью); |
| RAB - | одномерный рабочий массив длины 2 * K2 + 10 * K + 4 * M * K + 6 * М +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], чтобы X ∈ [xS-1, xS], и определить нормализованный аргумент α по формуле α = (X - xS-1) / H = (X - xS-1) / (xS - xS-1). Далее, зная α и S и пользуясь соответствующим набором коэффициентов из массива BY, найти значение решения Y(X) по формуле: Y(X) = (y1(X), ... , yM(X)),
K+2
yN(X) = yN (xS-1 + αH) = ∑ ' BY (N, I + 1, S) * 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) = ∑ ' BDY (N, I + 1, S) * 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) = ∑ ' BD2Y (N, I + 1, S) * TI*(α), N = 1, 2, ... , M.
I=0
При работе подпрограммы значения параметров M, XN, YN, DYN, XK, K, INIAPR, IMAX, IORDY сохраняются. Значение параметра H сохраняется, если он задан с учетом направления интегрирования, иначе его знак меняется на противоположный. Если после работы подпрограммы нет необходимости иметь начальные значения решения YN и/или производной DYN, то параметры YN и Y, а в случае производной параметры DYN и DY при обращении к ней можно совместить. Так как при интегрировании системы уравнений с помощью подпрограммы DE88E используются глобальные записи (структуры данных) с именами _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] в смещенные ряды Чебышёва, коэффициенты которых выражаются через цилиндрические функции:
∞
(9) y1(x) = 3 + J0 (q) + 2 ∑ (-1)i J2i (q) T2i* (x) ,
i=1
∞
(10) y2(x) = 2 + 2 ∑ (-1)i J2i+1 (q) T2i+1* (x) .
i=0
Интервал интегрирования не разбивается на элементарные сегменты, поэтому параметр H полагается равным 1. Вычисляются коэффициенты Чебышёва обеих компонент решения системы y1(x), y2(x), коэффициенты Чебышёва их первых y'1(x), y'2(x) и вторых производных y''1(x), y''2(x) на отрезке [0, 1]. Определяются также решение Y системы и его первая производная DY в конце интервала интегрирования xf. Приводятся подпрограмма вычисления значений правой части системы и фрагмент вызывающей программы, а также результаты счета, включая точное значение YT решения в точке xf, абсолютную погрешность DELY приближенного решения Y, вычисленного в точке xf, точное значение первой производной DYT решения в точке xf, абсолютную погрешность DELDY приближенного значения DYT, значения параметров H, K, INIAPR, IMAX, IORDY, NX и XS. Даются также точные значения коэффициентов Чебышёва для компонент решения y1(x), y2(x) в (9), (10) и абсолютные погрешности их приближенных значений, вычисленных программой DE88E.
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.
. . . . . . . . . . . . . . . . . . . . . . . . . .
function tde88e1: String;
var
M,IORDY,K0,IMAX,INIAPR,NXP1,NX,I,_i,KP3,KP2,KP1,L,J :Integer;
XN,H,XK,H0,X1,Q :Extended;
XS :Array [0..4] of Extended;
AY :Array [0..27] of Extended;
ADY :Array [0..25] of Extended;
AD2Y :Array [0..23] of Extended;
RАВ :Array [0..456] 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
_20;
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;
IORDY := 2;
H0 := 1.e0;
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МАХ := 16;
INIAPR := 1;
H := H0;
DE88E(F88EX1,M,XN,YN,DYN,XK,K0,INIAPR,IMAX,IORDY,H,Y,DY,NX,XS,AY,
ADY,AD2Y,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];
NXP1 := NX+1;
// Операторы вывода на печать: H0, K, INIAPR, IMAX, IORDY,Y, YT, DELY
. . . . . . . . . . . . . . . . . . . . . . . . . .
end;
Результаты:
------------------------------------------------------------------------
NX= 1
XS= 0,00000000000000000E+000
1,00000000000000000E+000
H0= 1,0000000000000000 K=11 INIAPR=1 IMAX=16 IORDY=2
Y = 3,87758256189037272E+000 2,47942553860420300E+000
YT = 3,87758256189037272E+000 2,47942553860420300E+000
DELY = 4,33680868994201774E-019 -4,33680868994201774E-019
DY = -4,79425538604203010E-001 8,77582561890372720E-001
DYT= -4,79425538604203000E-001 8,77582561890372716E-001
DELDY= 9,26992857475106291E-018 -3,68628738645071508E-018
COEFFICIENTS AY,ADY AND AD2Y ON 1 SEGMENT
-------------------------------------------------------------------------------------
Number of Chebyshev coefficients
coefficient --------------------------------------------------------------------
for Y for Y' for Y''
----------------------------------------------------------------------------------------------------------------
FOR 1 COMPONENT
0 7,87693961448162581E+000 -1,62969139051727385E-018 -1,87693961448162585E+000
1 -1,43995601033231058E-020 -4,84536915349747775E-001 -3,57041327926632679E-017
2 -6,12080469173652827E-002 -1,57971644662927013E-018 6,12080469173652526E-002
3 -6,29274893731007294E-020 5,12745998917448695E-003 -2,30664012196291068E-017
4 3,21472952728575153E-004 -8,24586574152061380E-019 -3,21472952728590798E-004
5 -2,66744792306166748E-020 -1,61072544827154893E-005 -9,87301603319612475E-018
6 -6,72136925723786440E-007 -2,91096989539727883E-019 6,72136925718989410E-007
7 -7,68127794727337242E-021 2,40317346553852059E-008 -2,88668828424265556E-018
8 7,51644630955338995E-010 -7,60212070160734554E-020 -7,51644631796354479E-010
9 1,69977160255607756E-019 -2,08935351856419877E-011 -4,54009659728304982E-019
10 -5,22635331281252579E-013 -6,19519897621795268E-018 5,22634886757078476E-013
11 -2,57918323980800914E-019 1,18780656081154199E-014 2,47353949388989802E-016
12 2,47459700169071248E-016 5,15320727893728755E-018
13 9,91001399795632220E-020
FOR 2 COMPONENT
0 4,00000000000000000E+000 1,87693961448162581E+000 1,76216734346784643E-017
1 4,84536915349747773E-001 9,29195143137967472E-019 -4,84536915349747756E-001
2 4,61280025858904393E-020 -6,12080469173652819E-002 1,39048928621265944E-017
3 -5,12745998917448812E-003 5,60171122450843958E-019 5,12745998917449884E-003
4 1,95735196894681237E-020 3,21472952728575581E-004 7,18283939271646688E-018
5 1,61072544827149609E-005 2,46994807419353979E-019 -1,61072544827104512E-005
6 7,57388627003220220E-021 -6,72136925723637221E-007 2,24294324432938730E-018
7 -2,40317346555225512E-008 6,52215369385811261E-020 2,40317346568421025E-008
8 -3,58518459271872262E-021 7,51644630994214278E-010 4,16740210049115767E-019
9 2,08935351719506413E-011 1,79947443905580250E-019 -2,08935349727543736E-011
10 7,94264505818028649E-021 -5,22635196008807519E-013 -6,06136777055177323E-018
11 -1,18836942594794924E-014 -1,37758358421631210E-019 1,18728675979271059E-014
12 -2,86996580045065020E-021 2,47351408290148039E-016
13 4,75675785173361614E-018
*************************************************************************************
Number of Chebyshev coefficients
coefficient --------------------------------------------------------------------
approximate exact absolute error
----------------------------------------------------------------------------------------------------------------
FOR 1 COMPONENT
0 7,87693961448162581E+000 7,87693961448162581E+000 4,33680868994201774E-019
1 -1,43995601033231058E-020 0,00000000000000000E+000 1,43995601033231058E-020
2 -6,12080469173652827E-002 -6,12080469173652826E-002 8,13151629364128326E-020
3 -6,29274893731007294E-020 0,00000000000000000E+000 6,29274893731007294E-020
4 3,21472952728575153E-004 3,21472952728575194E-004 4,11340375010369602E-020
5 -2,66744792306166748E-020 0,00000000000000000E+000 2,66744792306166748E-020
6 -6,72136925723786440E-007 -6,72136925723769759E-007 1,66806140400146487E-020
7 -7,68127794727337242E-021 0,00000000000000000E+000 7,68127794727337242E-021
8 7,51644630955338995E-010 7,51644630959521991E-010 4,18299571115677884E-021
9 1,69977160255607756E-019 0,00000000000000000E+000 -1,69977160255607756E-019
10 -5,22635331281252579E-013 -5,22635472164560617E-013 -1,40883308038781575E-019
11 -2,57918323980800914E-019 0,00000000000000000E+000 2,57918323980800914E-019
12 2,47459700169071248E-016 2,47676511895986538E-016 2,16811726915289745E-019
13 9,91001399795632220E-020 0,00000000000000000E+000 -9,91001399795632220E-020
FOR 2 COMPONENT
0 4,00000000000000000E+000 4,00000000000000000E+000 0,00000000000000000E+000
1 4,84536915349747773E-001 4,84536915349747773E-001 -8,13151629364128326E-020
2 4,61280025858904393E-020 0,00000000000000000E+000 -4,61280025858904393E-020
3 -5,12745998917448812E-003 -5,12745998917448815E-003 -2,83756037330190614E-020
4 1,95735196894681237E-020 0,00000000000000000E+000 -1,95735196894681237E-020
5 1,61072544827149609E-005 1,61072544827149482E-005 -1,27369270720915201E-020
6 7,57388627003220220E-021 0,00000000000000000E+000 -7,57388627003220220E-021
7 -2,40317346555225512E-008 -2,40317346555260458E-008 -3,49450527708696004E-021
8 -3,58518459271872262E-021 0,00000000000000000E+000 3,58518459271872262E-021
9 2,08935351719506413E-011 2,08935351786579598E-011 6,70731849333656507E-021
10 7,94264505818028649E-021 0,00000000000000000E+000 -7,94264505818028649E-021
11 -1,18836942594794924E-014 -1,18837079244649228E-014 -1,36649854303880382E-020
12 -2,86996580045065020E-021 0,00000000000000000E+000 2,86996580045065020E-021
13 4,75675785173361614E-018 4,76464654243100702E-018 7,88869069739088286E-021
----------------------------------------------------------------
2) Решается задача Коши из примера 1 с разбиением интервала интегрирования [0, 1] на два элементарных сегмента длиной H = 0,5. Так как значение H здесь уменьшено по сравнению с примером 1, то и число итераций IMAX также уменьшается. Приводятся фрагмент вызывающей программы и результаты счета, включающие те же данные, что и в примере 1, а именно: значения параметров NX, XS, H, K, INIAPR, IMAX, IORDY, точное значение YT решения в точке xf, абсолютную погрешность DELY приближенного решения Y, вычисленного подпрограммой DE88E в точке xf, точное значение первой производной DYT в точке xf, абсолютную погрешность DELDY приближенного значения DYT. Далее для каждого из двух элементарных сегментов приводятся коэффициенты Чебышёва для обеих компонент решения y1(x), y2(x), коэффициенты Чебышёва для их первых производных y'1(x), y'2(x) и вторых производных y''1(x), y''2(x). Все перечисленные данные приводятся при двух способах определения начального приближения коэффициентов Чебышёва второй производной решения, т.е. при INIAPR = 1 и при INIAPR = 2.
function tde88e2: String;
var
M,IORDY,K0,IMAX,INIAP0,INIAPR,NXP1,NX,I,_i,KP3,KP2,KP1,L,J :Integer;
XN,H,XK,H0,X1,Q :Extended;
XS :Array [0..2] of Extended;
AY :Array [0..55] of Extended;
ADY :Array [0..51] of Extended;
AD2Y :Array [0..47] of Extended;
RAB :Array [0..456] 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
_20,_30,_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);
XK := 1.e0;
IORDY := 2;
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;
DE88E(F88EX1,M,XN,YN,DYN,XK,K0,INIAPR,IMAX,IORDY,H,Y,DY,NX,XS,AY,
ADY,AD2Y,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];
NXP1 := NX+1;
// Операторы вывода на печать: H0, K, INIAPR, IMAX, IORDY,Y, YT, DELY
. . . . . . . . . . . . . . . . . . . . . . . . . .
end;
Результаты:
------------------------------------------------------------------------
NX= 2
XS= 0,00000000000000000E+000
5,00000000000000000E-001
1,00000000000000000E+000
H0= 0,5000000000000000 K=11 INIAPR=1 IMAX=13 IORDY=2
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
-------------------------------------------------------------------------------------
Number of Chebyshev coefficients
coefficient --------------------------------------------------------------------
for Y for Y' for Y''
----------------------------------------------------------------------------------------------------------------
FOR 1 COMPONENT
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
FOR 2 COMPONENT
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
-------------------------------------------------------------------------------------
Number of Chebyshev coefficients
coefficient --------------------------------------------------------------------
for Y for Y' for Y''
----------------------------------------------------------------------------------------------------------------
FOR 1 COMPONENT
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
FOR 2 COMPONENT
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
-------------------------------------------------------------------------------------
NX= 2
XS= 0,00000000000000000E+000
5,00000000000000000E-001
1,00000000000000000E+000
H0= 0,5000000000000000 K=11 INIAPR=2 IMAX=13 IORDY=2
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
COEFFICIENTS AY,ADY AND AD2Y ON 1 SEGMENT
-------------------------------------------------------------------------------------
Number of Chebyshev coefficients
coefficient --------------------------------------------------------------------
for Y for Y' for Y''
----------------------------------------------------------------------------------------------------------------
FOR 1 COMPONENT
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
FOR 2 COMPONENT
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
-------------------------------------------------------------------------------------
Number of Chebyshev coefficients
coefficient --------------------------------------------------------------------
for Y for Y' for Y''
----------------------------------------------------------------------------------------------------------------
FOR 1 COMPONENT
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
FOR 2 COMPONENT
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
----------------------------------------------------------------
3) Решается задача Коши
(11) y'' = - 4q(tg y)y' / (1+tg2 y), y(0) = - arctg q, y'(0) = 2q / (1 + q2), 0 ≤ x ≤ xf, xf = 1, q = 1 / 16.
Точное решение задачи y(x) = arctg (q (2x -1)) = arctg (x / 8 - 1 /16) раскладывается на отрезке [0, 1] в смещенный ряд Чебышёва
∞
(12) y(x) = arctg (q (2x - 1)) = 2 ∑ (-1)i (p2i+1) / (2i+1) T*2i+1 (x), p = (√(1+q2) - 1) / q, q ≠ 0.
i=0
Интервал интегрирования не разбивается на элементарные сегменты, поэтому параметр H полагается равным 1. Вычисляются коэффициенты Чебышёва решения y(x) и коэффициенты Чебышёва его производных y'(x), y''(x) на отрезке [0, 1]. Определяется также решение Y задачи и его первая производная DY в конце интервала интегрирования xf. Приводятся подпрограмма вычисления значений правой части и фрагмент вызывающей программы, а также результаты счета, включая точные значения YT, DYT решения и производной в точке xf, абсолютные погрешности DELY, DELDY приближенного решения Y и производной DY, вычисленных в точке xf, значения параметров NX, XS, H, K, INIAPR, IMAX, IORDY. Даются также точные значения коэффициентов Чебышёва решения y(x) в (12) и абсолютные погрешности их приближенных значений, вычисленных подпрограммой DE88E.
Unit FATAN2M_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure FATAN2M(X :Extended; var Y :Array of Extended; var DY :Array of Extended; var Z :Array of Extended;
M :Integer);
implementation
procedure FATAN2M(X :Extended; var Y :Array of Extended; var DY :Array of Extended; var Z :Array of Extended;
M :Integer);
var
R :Extended;
begin
{ ---------------------------------------- }
R := Tan(Y[0]);
Z[0] := -4.e0*_BLOCKQ.elm1*R*DY[0]/(1.e0+R*R);
end;
end.
. . . . . . . . . . . . . . . . . . . . . . . . . .
function tde88e3: String;
var
M,IORDY,K0,IMAX,INIAPR,NXP1,NX,I,KP1,KP2,KP3,J,KD2,INCH :Integer;
YN,DYN,Y,DY,YT,DYT,DELY,DELDY,H0,H,XN,XK,R,W,P :Extended;
AYT :Array [0..12] of Extended;
DELAY :Array [0..12] of Extended;
XS :Array [0..1] of Extended;
AY :Array [0..12] of Extended;
ADY :Array [0..11] of Extended;
AD2Y :Array [0..10] of Extended;
RAB :Array [0..350] of Extended;
const
K :Integer = 10;
label
_31,_32,_34;
begin
Result := ''; { результат функции }
{ PROGRAM TDEATA }
_BLOCKQ.elm1 := 0.0625e0;
M := 1;
XN := 0.e0;
YN := -ArcTan(_BLOCKQ.elm1);
DYN := 2.e0*_BLOCKQ.elm1/(1.e0+_BLOCKQ.elm1*_BLOCKQ.elm1);
IORDY := 2;
H0 := 1.e0;
XK := 1.e0;
YT := ArcTan(_BLOCKQ.elm1);
DYT := DYN;
K0 := K;
IMAX := 5;
INIAPR := 1;
H := H0;
DE88E(FATAN2M,M,XN,YN,DYN,XK,K0,INIAPR,IMAX,IORDY,H,Y,DY,NX,XS,AY,
ADY,AD2Y,RAB);
DELY := YT -Y;
DELDY := DYT-DY;
NXP1 := NX+1;
// Операторы вывода на печать: H0, K, INIAPR, IMAX, IORDY,Y, YT, DELY
. . . . . . . . . . . . . . . . . . . . . . . . . .
end;
Результаты:
------------------------------------------------------------------------
NX= 1
XS= 0,00000000000000000E+000
1,00000000000000000E+000
H0= 1,0000000000000000 K= 10 INIAPR=1 IMAX= 5 IORDY= 2
Y= 6,24188099959573485E-002
YT= 6,24188099959573485E-002
DELTA Y = -6,77626357803440271E-021
DY= 1,24513618677042802E-001
DYT= 1,24513618677042802E-001
DELTA DY = -1,21972744404619249E-019
COEFFICIENTS AY,ADY AND AD2Y ON 1 SEGMENT
Number of Chebyshev coefficients
coefficient --------------------------------------------------------------------
for Y for Y' for Y''
----------------------------------------------------------------------------------------------------------------
0 8,33798057453451896E-021 2,49513144620722138E-001 1,27504928341373898E-019
1 6,24390837627947298E-002 4,82411733240925740E-021 -1,94173653401340173E-003
2 1,54131320805714157E-022 -2,43190430456780992E-004 1,08208459011736868E-019
3 -2,02856215326620923E-005 3,59106676596354415E-021 3,78690964084620576E-006
4 1,06444358492032111E-022 2,37027935164115572E-007 6,51156578201743386E-020
5 1,18629478381439031E-008 1,88795703009103037E-021 -5,53732177964338148E-009
6 6,09986617428389837E-023 -2,31021598762490090E-010 2,73565172183537313E-020
7 -8,25881307954036154E-012 4,23989148262894757E-022 7,19659065638068354E-012
8 -8,90195362915827252E-021 2,25167464640033060E-013 1,54848210669926781E-020
9 6,26074083173728942E-015 2,85286505281327615E-019 -8,76821210037437461E-015
10 1,29587702280904640E-020 -2,19205302509359365E-016 -1,02548293690608015E-017
11 -4,98193869339453103E-018 -2,33064303842290943E-019
12 -4,85550633004772797E-021
Number of Chebyshev coefficients for Y
coefficient --------------------------------------------------------------------
approximate exact absolute error
----------------------------------------------------------------------------------------------------------------
0 8,33798057453451896E-021 0,00000000000000000E+000 -4,16899028726725948E-021
1 6,24390837627947298E-002 6,24390837627947296E-002 -1,99899775552014880E-019
2 1,54131320805714157E-022 0,00000000000000000E+000 -1,54131320805714157E-022
3 -2,02856215326620923E-005 -2,02856215326620922E-005 6,12113653289240479E-023
4 1,06444358492032111E-022 0,00000000000000000E+000 -1,06444358492032111E-022
5 1,18629478381439031E-008 1,18629478381438272E-008 -7,59358264672214976E-023
6 6,09986617428389837E-023 0,00000000000000000E+000 -6,09986617428389837E-023
7 -8,25881307954036154E-012 -8,25881307956705560E-012 -2,66940574902166506E-023
8 -8,90195362915827252E-021 0,00000000000000000E+000 8,90195362915827252E-021
9 6,26074083173728942E-015 6,26074793977291928E-015 7,10803562985504982E-021
10 1,29587702280904640E-020 0,00000000000000000E+000 -1,29587702280904640E-020
11 -4,98193869339453103E-018 -4,99262670434113154E-018 -1,06880109466005119E-020
12 -4,85550633004772797E-021 0,00000000000000000E+000 4,85550633004772797E-021
--------------------------------------------------------------------------
Заметим, что в четвертой колонке последней таблицы (absolute error) для нулевого коэффициента Чебышёва a0*[y] решения задачи Коши (11) разность между его точным значением, равным нулю, и приближенным значением, вычисленным подпрограммой DE88E, приводится с дополнительным множителем 1/2 (т.е. в таблице в четвертой колонке указана величина -1/2a0*[y].
4) Решается обыкновенное дифференциальное уравнение из примера 3, но начальное условие задается в точке XN = xf = 1:
(13) y'' = - 4q(tg y)y' / (1+tg2y), y(1) = arctg q, y'(1) = 2q / (1 + q2), q = 1 / 16,
а концом интервала интегрирования теперь является точка XK = 0. В силу указанных в (11) и (13) начальных условий решение задачи Коши (13) совпадает с решением задачи Коши (11); решением этих двух задач на [0, 1] является одна и та же функция
(14) y(x) = arctg (x / 8 - 1 / 16), 0 ≤ x ≤ 1.
Так как y(x) удовлетворяет соотношению y(1 - α) = - y(α), 0 ≤ α ≤ 1, то коэффициенты Чебышёва функции y(x) = y(1 - α), x = 1 - α, на [0, 1] отличаются знаком от коэффициентов Чебышёва функции y(x) = y(α), x = α, на том же отрезке [0, 1] и равны им по абсолютной величине. Можно показать, что аналогичным свойством обладают коэффициенты Чебышёва и для вторых производных этих функций. Кроме того, функция y(x) = y(α), x = α, 0 ≤ α ≤ 1, (определяемая по формуле (14)) и функция y(x) = y(1 - α), x = 1 - α, 0 ≤ α ≤ 1, имеют одинаковые коэффициенты Чебышёва для первых производных этих функций на [0, 1]:
a*i {[y (1 - α)]'x=1-α} = a*i {[y (α)]'x=α}, i = 0, 1, ... .
В примере 4 интегрирование выполняется справа налево без разбиения интервала интегрирования [XN, XK] на элементарные сегменты, поэтому параметр H = - 1. При этом вычисляются коэффициенты Чебышёва функции y(x) = y(1 - α), x = 1 - α, на [0, 1] и коэффициенты Чебышёва ее производных y'(x), y''(x). Определяются также решение Y задачи (13) и его первая производная DY в конце интервала интегрирования XK. Приводятся фрагмент вызывающей программы и результаты счета, включая точные значения YT, DYT решения и производной в точке XK = 0, абсолютные погрешности DELY, DELDY приближенного решения Y и производной DY, вычисленных в точке XK = 0, значения параметров NX, XS, H, K, INIAPR, IMAX, IORDY.
function tde88e4: String;
var
M,IORDY,K0,IMAX,INIAPR,NXP1,NX,I,KP1,KP2,KP3,J :Integer;
YN,DYN,Y,DY,YT,DYT,DELY,DELDY,H0,H,XN,XK :Extended;
XS :Array [0..1] of Extended;
AY :Array [0..12] of Extended;
ADY :Array [0..11] of Extended;
AD2Y :Array [0..10] of Extended;
RAB :Array [0..350] of Extended;
const
K :Integer = 10;
begin
Result := ''; { результат функции }
{ PROGRAM TDEATA }
_BLOCKQ.elm1 := 0.0625e0;
M := 1;
XN := 1.e0;
YN := ArcTan(_BLOCKQ.elm1);
DYN := 2.e0*_BLOCKQ.elm1/(1.e0+_BLOCKQ.elm1*_BLOCKQ.elm1);
IORDY := 2;
H0 := -1.e0;
XK := 0.e0;
YT := -ArcTan(_BLOCKQ.elm1);
DYT := DYN;
K0 := K;
IMAX := 5;
INIAPR := 1;
H := H0;
DE88E(FATAN2M,M,XN,YN,DYN,XK,K0,INIAPR,IMAX,IORDY,H,Y,DY,NX,XS,AY,
ADY,AD2Y,RAB);
DELY := YT -Y;
DELDY := DYT-DY;
NXP1 := NX+1;
// Операторы вывода на печать: H0, K, INIAPR, IMAX, IORDY,Y, YT, DELY
. . . . . . . . . . . . . . . . . . . . . . . . . .
end;
Результаты:
------------------------------------------------------------------------
NX= 1
XS= 1,00000000000000000E+000
0,00000000000000000E+000
H0= -1,0000000000000000 K= 10 INIAPR=1 IMAX= 5 IORDY= 2
Y= -6,24188099959573485E-002
YT= -6,24188099959573485E-002
DELTA Y = 6,77626357803440271E-021
DY= 1,24513618677042802E-001
DYT= 1,24513618677042802E-001
DELTA DY = -1,21972744404619249E-019
COEFFICIENTS AY,ADY AND AD2Y ON 1 SEGMENT
Number of Chebyshev coefficients
coefficient --------------------------------------------------------------------
for Y for Y' for Y''
----------------------------------------------------------------------------------------------------------------
0 -8,33798057453451896E-021 2,49513144620722138E-001 -1,27504928341373898E-019
1 -6,24390837627947298E-002 4,82411733240925740E-021 1,94173653401340173E-003
2 -1,54131320805714157E-022 -2,43190430456780992E-004 -1,08208459011736868E-019
3 2,02856215326620923E-005 3,59106676596354415E-021 -3,78690964084620576E-006
4 -1,06444358492032111E-022 2,37027935164115572E-007 -6,51156578201743386E-020
5 -1,18629478381439031E-008 1,88795703009103037E-021 5,53732177964338148E-009
6 -6,09986617428389837E-023 -2,31021598762490090E-010 -2,73565172183537313E-020
7 8,25881307954036154E-012 4,23989148262894757E-022 -7,19659065638068354E-012
8 8,90195362915827252E-021 2,25167464640033060E-013 -1,54848210669926781E-020
9 -6,26074083173728942E-015 2,85286505281327615E-019 8,76821210037437461E-015
10 -1,29587702280904640E-020 -2,19205302509359365E-016 1,02548293690608015E-017
11 4,98193869339453103E-018 -2,33064303842290943E-019
12 4,85550633004772797E-021
---------------------------------------------------------------------
Из приведенных результатов видно, что на выходе из подпрограммы DE88E в массиве XS точки XN, XK располагаются в порядке убывания. Вычисленные подпрограммой коэффициенты Чебышёва функции y(x) = y(1 - α), x = 1 - α, на [0, 1] отличаются знаком (как и должно быть) от коэффициентов Чебышёва функции y(x) из примера 3, а коэффициенты Чебышёва первой производной y'(x) совпадают с коэффициентами Чебышёва первой производной функции y(x) из примера 3 (как и должно быть). Аналогичным свойством обладают коэффициенты Чебышёва для вторых производных этих функций.
5) Решается задача Коши
(15) y''(x) = 1536x2 + 192x + 16, y(0) = 1, y'(0) = 2, 0 ≤ x ≤ 1.
Точное решение задачи y(x) = 128x4 + 32x3 + 8x2 +2x + 1 = 50 + 76T*1(x) + 35T*2(x) + 9T*3(x) + T*4(x). Интервал интегрирования не разбивается на элементарные сегменты, поэтому параметр H полагается равным 1. Вычисляются коэффициенты Чебышёва решения y(x), коэффициенты Чебышёва его первой производной y'(x) = 206 + 296T*1(x) + 108T*2(x) + 16T*3(x) на [0, 1] и коэффициенты Чебышёва второй производной y''(x) = 688 + 864T*1(x) + 192T*2(x) на [0, 1]. Приводятся подпрограмма вычисления значений правой части, фрагмент вызывающей программы и результаты счета, аналогичные результатам из примеров 1-4. Поскольку правая часть дифференциального уравнения не зависит от y и y', то параметр IMAX можно положить равным 1.
Unit FPOLIN4_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure FPOLIN4(X :Extended; var Y :Array of Extended; var DY :Array of Extended; var Z :Array of Extended;
M :Integer);
implementation
procedure FPOLIN4(X :Extended; var Y :Array of Extended; var DY :Array of Extended; var Z :Array of Extended;
M :Integer);
begin
Z[0] := 1536.e0*X*X+192.e0*X+16.e0;
end;
end.
. . . . . . . . . . . . . . . . . . . . . . . . . .
function tde88e5: String;
var
M,IORDY,K0,IMAX,INIAPR,NXP1,NX,I,KP1,KP2,KP3,J :Integer;
YN,DYN,Y,DY,YT,DYT,DELY,DELDY,H0,H,XN,XK :Extended;
XS :Array [0..1] of Extended;
AY :Array [0..4] of Extended;
ADY :Array [0..3] of Extended;
AD2Y :Array [0..2] of Extended;
RAB :Array [0..46] of Extended;
const
K :Integer = 2;
begin
Result := ''; { результат функции }
{ PROGRAM POL4 }
M := 1;
XN := 0.e0;
YN := 1.e0;
DYN := 2.e0;
IORDY := 2;
H0 := 1.e0;
XK := 1.e0;
YT := 171.e0;
DYT := 626.e0;
K0 := K;
IMAX := 1;
INIAPR := 1;
H := H0;
DE88E(FPOLIN4,M,XN,YN,DYN,XK,K0,INIAPR,IMAX,IORDY,H,Y,DY,NX,XS,AY,
ADY,AD2Y,RAB);
DELY := YT -Y;
DELDY := DYT-DY;
NXP1 := NX+1;
// Операторы вывода на печать: H0, K, INIAPR, IMAX, IORDY,Y, YT, DELY
. . . . . . . . . . . . . . . . . . . . . . . . . .
end;
Результаты:
------------------------------------------------------------------------
NX= 1
XS= 0,00000000000000000E+000
1,00000000000000000E+000
H0= 1,0000000000000000 K= 2 INIAPR=1 IMAX= 1 IORDY= 2
Y= 1,71000000000000000E+002
YT= 1,71000000000000000E+002
DELTA Y = -2,77555756156289135E-017
DY= 6,26000000000000000E+002
DYT= 6,26000000000000000E+002
DELTA DY = -5,55111512312578270E-017
COEFFICIENTS AY,ADY AND AD2Y ON 1 SEGMENT
Number of Chebyshev coefficients
coefficient --------------------------------------------------------------------
for Y for Y' for Y''
----------------------------------------------------------------------------------------------------------------
0 1,00000000000000000E+002 4,12000000000000000E+002 1,37600000000000000E+003
1 7,60000000000000000E+001 2,96000000000000000E+002 8,64000000000000000E+002
2 3,50000000000000000E+001 1,08000000000000000E+002 1,92000000000000000E+002
3 9,00000000000000000E+000 1,60000000000000000E+001
4 1,00000000000000000E+000
---------------------------------------------------------------------
6) Решается обыкновенное дифференциальное уравнение из примера 5, но начальное условие задается в точке XN = 1:
(16) y''(x) = 1536x2 + 192x + 16, y(1) = 171, y'(1) =626, 0 ≤ x ≤ 1,
а концом интервала интегрирования является точка XK = 0. В силу указанных в (15) и (16) начальных условий решение задачи Коши (16) совпадает с решением задачи Коши (15); решением этих двух задач является одна и та же функция
y(x) = 50 + 76T*1(x) + 35T*2(x) + 9T*3(x) + T*4(x), y'(x) = 206 + 296T*1(x) + 108T*2(x) + 16T*3(x).
Так как T*i (1-α) = - T*i (α), если i - нечетное, и T*i (1-α) = T*i (α), если i - четное (0 ≤ α ≤ 1), то коэффициенты Чебышёва функции y(x) = y(1 - α), x = 1 - α, 0 ≤ α ≤ 1, легко получаются из коэффициентов Чебышёва функции y(x) = y(α), x = α, на [0, 1], а коэффициенты Чебышёва ее первой и второй производных (по x) y'(x) = y'x=1-α(1 - α), y''(x) = y''x=1-α(1-α) легко получаются из коэффициентов Чебышёва первой и второй производных функции y(x) = y(α), x = α, 0 ≤ α ≤ 1:
(17) y(x) = y(1-α) = 50 - 76T*1(α) + 35T*2(α) - 9T*3(α) + T*4(α),
(18) y'(x) = y'(1 - α) = 206 - 296T*1(α) + 108T*2(α) - 16T*3(α).
(19) y''(x) = y''(1 - α) = 688 - 864T*1(α) + 192T*2(α).
В примере 6 интегрирование выполняется справа налево без разбиения интервала интегрирования [XN, XK] на элементарные сегменты, поэтому параметр H = - 1. При этом вычисляются коэффициенты Чебышёва функции y(x) = y(1 - α), x = 1 - α, на [0, 1] и коэффициенты Чебышёва ее производных y'(x), y''(x). Приводятся фрагмент вызывающей программы и результат счета, аналогичные результатам счета из примера 5. Поскольку правая часть дифференциального уравнения не зависит от y и y', то параметр IMAX (как и в примере 5) следует положить равным 1.
function tde88e6: String;
var
M,IORDY,K0,IMAX,INIAPR,NXP1,NX,I,KP1,KP2,KP3,J :Integer;
YN,DYN,Y,DY,YT,DYT,DELY,DELDY,H0,H,XN,XK :Extended;
XS :Array [0..1] of Extended;
AY :Array [0..4] of Extended;
ADY :Array [0..3] of Extended;
AD2Y :Array [0..2] of Extended;
RAB :Array [0..46] of Extended;
const
K :Integer = 2;
begin
Result := ''; { результат функции }
{ PROGRAM POL4 }
M := 1;
XN := 1.e0;
YN := 171.e0;
DYN := 626.e0;
IORDY := 2;
H0 := -1.e0;
XK := 0.e0;
YT := 1.e0;
DYT := 2.e0;
K0 := K;
IMAX := 1;
INIAPR := 1;
H := H0;
DE88E(FPOLIN4,M,XN,YN,DYN,XK,K0,INIAPR,IMAX,IORDY,H,Y,DY,NX,XS,AY,
ADY,AD2Y,RAB);
DELY := YT -Y;
DELDY := DYT-DY;
NXP1 := NX+1;
// Операторы вывода на печать: H0, K, INIAPR, IMAX, IORDY,Y, YT, DELY
. . . . . . . . . . . . . . . . . . . . . . . . . .
end;
Результаты:
------------------------------------------------------------------------
NX= 1
XS= 1,00000000000000000E+000
0,00000000000000000E+000
H0= -1,0000000000000000 K= 2 INIAPR=1 IMAX= 1 IORDY= 2
Y= 9,99999999999999941E-001
YT= 1,00000000000000000E+000
DELTA Y = 5,89805981832114412E-017
DY= 2,00000000000000007E+000
DYT= 2,00000000000000000E+000
DELTA DY = -6,93889390390722838E-017
COEFFICIENTS AY,ADY AND AD2Y ON 1 SEGMENT
Number of Chebyshev coefficients
coefficient --------------------------------------------------------------------
for Y for Y' for Y''
----------------------------------------------------------------------------------------------------------------
0 9,99999999999999999E+001 4,12000000000000000E+002 1,37600000000000000E+003
1 -7,60000000000000000E+001 -2,96000000000000000E+002 -8,64000000000000000E+002
2 3,50000000000000000E+001 1,08000000000000000E+002 1,92000000000000000E+002
3 -9,00000000000000000E+000 -1,60000000000000000E+001
4 1,00000000000000000E+000
---------------------------------------------------------------------
Из приведенных результатов видно, что на выходе из подпрограммы DE88E в массиве XS точки XN, XK располагаются в порядке убывания. Вычисленные подпрограммой коэффициенты Чебышёва функции y(x) = y(1 - α), x = 1 - α, 0 ≤ α ≤ 1, и коэффициенты Чебышёва ее производных y'(x), y''(x) находятся в полном согласии с представлениями (17), (18), (19) (как и должно быть).