|
Текст подпрограммы и версий (Си) de53d_c.zip |
Тексты тестовых примеров (Си) tde53d1_c.zip , tde53d2_c.zip , tde53d3_c.zip |
Вычисление коэффициентов разложения решения задачи Коши для обыкновенных дифференциальных уравнений первого порядка и его производной в ряд по смещенным многочленам Чебышёва первого рода с встроенной пошаговой обработкой результатов (получение аналитического выражения для приближенного решения задачи Коши и его производной в виде частичной суммы смещенного ряда Чебышёва).
Решается задача Коши для системы M обыкновенных дифференциальных уравнений первого порядка
(1) Y ' = F (X, Y) ,
Y = ( y1, ... , yM ) , F = ( f1 (X, y1, ... , yM), ... , fM (X, y1, ... , yM) ) , XN ≤ X ≤ XK
с начальными условиями, заданными в точке XN :
(2) Y(XN) = YN , YN = ( y10, ... , yM0 ) ,
при условии, что правая часть системы (1) имеет непрерывные ограниченные частные производные по переменным X, Y. Предполагается, что на отрезке [XN, XK] задача (1),(2) имеет единственное решение. Тогда решение задачи Коши и его производная Y'(XN + αΔ) = F(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) Φ(α) = ∑ ' 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)) является быстросходящимся, то его сумма на [XN, XK] (и сумма ряда (4)) хорошо приближается частичной суммой некоторого порядка. Эта частичная сумма принимается в качестве приближенного аналитического решения задачи (1), (2) на промежутке [XN, XK]. В противном случае, т.е. при медленной сходимости ряда (3) на интервале [XN, XK], получение аналитического решения в виде одной частичной суммы на всем отрезке интегрирования [XN, XK] может быть затруднено. Поэтому целесообразно использовать разбиение промежутка интегрирования [XN, XK] на такие элементарные сегменты некоторой длины H: [xs, xs+H], x0 = XN, s = 0, 1, ... , на каждом из которых ряды Чебышёва для решения Y(X) и его производной Y'(X) будут сходиться значительно быстрее. На каждом подобном сегменте решение исходной задачи Коши приближенно представляется в виде (K + 1) - й частичной суммы смещенного ряда Чебышёва
K+1 1
(5) Y (xs + αH) ≈ ∑ ' ai*[Y] Ti*(α) , 0 ≤ α ≤ 1 , ai*[Y] = 2/π ∫ Y(xs + αH) Ti* (α) / √ (α (1 - α) ) dα ,
i=0 0
а его производная - в виде частичной суммы K-го порядка
K 1
(6) Y' (xs + αH) ≈ ∑ ' ai*[Y ' ] Ti*(α) , 0 ≤ α ≤ 1 , ai*[Y ' ] = 2/π ∫ Y ' (xs + αH) Ti* (α) / √ (α (1 - α) ) dα .
i=0 0
В этом случае аналитическое решение задачи (1),(2) состоит из совокупности частичных сумм рядов Чебышёва, построенных на этих элементарных сегментах. Порядок этих частичных сумм и длина элементарных сегментов задаются пользователем при обращении к подпрограмме. Число NX элементарных сегментов, на которые разбивается промежуток интегрирования, равно [(XK-XN) / H], если длина промежутка интегрирования является целым кратным H. Если длина интервала интегрирования не является целым кратным H, то число элементарных сегментов NX равно [(XK-XN) / H] + 1; в этом случае последний элементарный сегмент считается нестандартным (квадратные скобки означают целую часть числа).
После вычисления коэффициентов Чебышёва для решения задачи Коши и его производной на каждом отдельном элементарном сегменте [xs, xs+H] подпрограмма de53d_c передает управление специальной подпрограмме обработки результатов для дальнейшей работы с полученными на этом сегменте коэффициентами Чебышёва ai*[Y], ai*[Y']. Характер и цель этой обработки определяет сам пользователь и выполняет ее с помощью составленной им же подпрограммы, имя которой задается при обращении к подпрограмме de53d_c. Например, он может использовать их для каких-то необходимых ему на данном элементарном сегменте вычислений, или запомнить вычисленные коэффициенты в удобном для себя виде, или просто их проигнорировать, если эти коэффициенты на данном элементарном сегменте ему не нужны. Таким образом, в подпрограмме обработки результатов можно, в частности, непосредственно приступить к применению полученных на элементарном сегменте [xs, xs+H] коэффициентов Чебышёва (или, по-другому, немедленно начать использование построенных на сегменте [xs, xs+H] частичных сумм (5), (6) рядов Чебышёва) решения и его производной, не дожидаясь окончания процедуры интегрирования системы дифференциальных уравнений (1) на всем промежутке [XN, XK]. Иными словами, подпрограмма de53d_c после вычисления коэффициентов Чебышёва на элементарном сегменте [xs, xs+H]=[xs, xs+1] предоставляет возможность приостановить процедуру интегрирования системы (1) на время обработки относящихся к этому сегменту результатов, а затем (после завершения такой обработки) продолжить дальнейшее интегрирование, т.е. перейти к вычислению коэффициентов на следующем элементарном сегменте [xs+1, xs+2] и т.д.
При разбиении промежутка интегрирования на элементарные сегменты решение задачи сводится к определению нескольких наборов коэффициентов ai*[Y], i = 0, 1, ... , K + 1. Коэффициенты ai*[Y] ряда Чебышёва для решения на сегменте [xs, xs + H] выражаются через коэффициенты ai*[Φ] ряда Чебышёва его производной Φ(α) = F(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 + 2) при H --> 0. Если H подобрано достаточно малым, то хорошая точность приближенного решения может быть получена и при меньшем числе итераций. Вообще, число итераций зависит от K и H. С увеличением H число итераций может также возрастать.
Начальное приближение коэффициентов ai*[Φ] ряда Чебышёва для производной на сегменте [xs, xs + H] может вычисляться двумя способами. В первом способе начальное приближение определяется только с использованием значения решения 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] всегда применяется исключительно первый способ. Способ выбора начального приближения задается пользователем при обращении к подпрограмме.
В дальнейшем при описании параметров подпрограммы коэффициенты ряда Чебышёва будем называть коэффициентами Чебышёва.
О.Б.Арушанян, С.Ф.Залеткин. Приближенное решение задачи Коши для обыкновенных дифференциальных уравнений методом рядов Чебышёва. Вычислительные методы и программирование: Новые вычислительные технологии (Электронный научный журнал) (17), 121 - 131, 2016.
О.Б.Арушанян, С.Ф.Залеткин. Использование рядов Чебышёва для приближенного аналитического решения обыкновенных дифференциальных уравнений. Вестник Московского университета. Серия 1: Математика. Механика. 5, 52 - 56, 2016.
О.Б.Арушанян, Н.И.Волченскова, С.Ф.Залеткин. Вычисление коэффициентов разложения решения задачи Коши в ряд по многочленам Чебышёва. Вестник Московского университета. Серия 1: Математика. Механика. 5, 24 - 30, 2012.
int de53d_c (U_fp f, S_fp ftreat, int *m, double *xn, double *yn, double *xk,
int *k, int *iniapr, int *imax, double *h, double *y, double *rab)
Параметры
| f - |
имя функции вычисления значений правой
части дифференциального уравнения. Первый
оператор функции должен иметь вид: int f (x, y , dy, m). Здесь: x, y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в dy. В случае системы уравнений, т.е. когда m ≠ 1 , параметры y и dy представляют массивы длины m (тип параметров x, y и dy: с двойной точностью); |
| ftreat - |
имя функции обработки результатов. Первый оператор функции
должен иметь вид: int ftreat (s, xi, xe, y, ay, ady, m, kp1, kp2). Здесь: |
| s - | номер элементарного сегмента, s=1, 2, ... , nx (тип: целый); |
| xi, xe - | начало и конец элементарного сегмента с данным номером s: [xs-1, xs-1+h], x0=xn, т.е. xi=xs-1, xe=xs-1+h (тип: с двойной точностью); |
| y - | значение решения задачи Коши, вычисленное функцией de53d_c в конце xe элементарного сегмента [xs-1, xs-1+h], т.е. в точке xs= xs-1+h (тип: с двойной точностью); |
| ay - | двумерный массив с измерениями m, kp2 (kp2 = k+2). Переменная с индексом ay (n, i+1) содержит вычисленный функцией de53d_c i-й коэффициент Чебышёва для N-й компоненты решения yN(X) на элементарном сегменте [xi, xe], i = 0, 1, ... , k+1; n = 1, ... , m (тип: с двойной точностью); |
| ady - | двумерный массив с измерениями m, kp1 (kp1 = k+1). Переменная с индексом ady (n, i+1) содержит вычисленный функцией de53d_c i-й коэффициент Чебышёва для N-й компоненты производной решения y'N(X) на элементарном сегменте [xi, xe], i = 0, 1, ... , k; n = 1, ... , m (тип: с двойной точностью); |
| m - | количество уравнений в системе (1) (тип: целый); |
| kp1 - | целый параметр, имеющий значение k+1; |
| kp2 - | целый параметр, имеющий значение k+2. |
| Таким образом, для системы уравнений (т.е. когда m ≠ 1) параметры y, ay, ady представляют массивы с регулируемыми измерениями, описатели которых имеют вид y(m), ay(m, kp2), ady(m, kp1). В случае, когда интегрируется одно скалярное уравнение (т.е. при m = 1), параметры ay, ady представляют массивы с регулируемыми измерениями, описатели которых могут иметь вид ay(kp2), ady(kp1), а параметр y может быть переменной (простой переменной). Данная функция ftreat составляется пользователем и выполняет нужную ему обработку содержащихся в ay, ady коэффициентов (см. "Математическое описание"). Обработка результатов на каждом элементарном сегменте в функции ftreat должна заканчиваться оператором return. При работе функции ftreat значения всех ее параметров не должны изменяться. | |
| m - | количество уравнений в системе (тип: целый); |
| xn, yn - | начальные значения аргумента и решения; в случае системы уравнений (т.е. когда m ≠ 1) yn представляет массив длины m (тип: с двойной точностью); |
| xk - | конец интервала интегрирования. На отрезке [xn, xk] вычисляется приближенное аналитическое решение задачи Коши в виде одной частичной суммы ряда Чебышёва либо в виде совокупности частичных сумм. xk может быть больше, меньше или равно xn (тип: с двойной точностью); |
| k - | порядок частичной суммы смещенного ряда Чебышёва, с помощью которой аппроксимируется производная решения задачи Коши на каждом элементарном сегменте разбиения интервала интегрирования; при этом само решение задачи Коши приближается на каждом элементарном сегменте частичной суммой (k + 1) - го порядка; k≥2 (см. "Математическое описание" и "Замечания по использованию"; тип: целый); |
| iniapr - |
целый указатель способа выбора начального
приближения коэффициентов Чебышёва для производной
решения на каждом элементарном сегменте: |
| iniapr=1 - |
для первого способа, когда начальное
приближение определяется только с
использованием значения решения в
начале каждого элементарного сегмента; |
| iniapr=2 - | для второго способа, когда начальное приближение коэффициентов Чебышёва на текущем элементарном сегменте (начиная со второго) определяется через коэффициенты Чебышёва, вычисленные на предыдущем элементарном сегменте, т.е. путем экстраполяции коэффициентов с предыдущего сегмента на следующий (см. "Математическое описание"); |
| imax - | целая переменная, задающая количество итераций, которое предполагается выполнить в итерационном процессе вычисления коэффициентов Чебышёва для производной решения задачи Коши на каждом элементарном сегменте, исходя из некоторого начального приближения, способ определения которого задается параметром iniapr; imax≥1. Для получения максимального порядка точности приближенного решения необходимо выполнить не менее k итераций (см. "Математическое описание" и "Замечания по использованию"); |
| h - | переменная с двойной точностью, содержащая значение длины элементарных сегментов, на которые разбивается интервал интегрирования (диаметр разбиения промежутка интегрирования или аналог шага интегрирования для одношаговых методов). Предполагается, что на каждом элементарном сегменте ряды Чебышёва для решения и его производной являются быстросходящимися рядами. Может задаваться с учетом направления интегрирования, т.е. положительным, если xk > xn, отрицательным, если xk < xn , или без такого учета в виде абсолютной величины; |
| y - | на выходе из функции содержит значение решения задачи Коши, вычисленное функцией при значении аргумента xk. Для системы уравнений (когда m ≠ 1) задается массивом длины m. В случае совпадения значений параметров xn и xk значение y полагается равным начальному значению yn (тип: с двойной точностью); |
| rab - | одномерный рабочий массив длины 2 * k2 + 7 * k + 5 * m * k + 8 * m + 1. (тип: с удвоенной точностью). |
Версии: нет
Вызываемые подпрограммы
| de44d_c - | выполнение одного шага приближенного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом рядов Чебышёва |
|
Кроме того, используются рабочие функции de70dk_c, de70dh_c, de70d0_c, de70di_c, de70df_c, de70dq_c, de71de_c, de70dp_c, de71dt_c, de71dp_c, de71di_c, de71df_c, de71ds_c, de70da_c, de70dc_c. |
Замечания по использованию
|
Разбиение промежутка интегрирования на элементарные сегменты (шаги) выполняется для того, чтобы на каждом таком сегменте ряды Чебышёва для решения и его производной были быстросходящимися рядами. Другими словами, длина элементарных сегментов, задаваемая параметром h, подбирается таким образом, чтобы убывание коэффициентов этих рядов Чебышёва на элементарном сегменте происходило достаточно быстро, вследствие чего можно было бы считать частичные суммы этих рядов близкими к многочленам наилучшего равномерного приближения на элементарном сегменте для решения и его производной. Порядок этих частичных сумм задается параметром k. Если начальное приближение для коэффициентов Чебышёва функции Φ(α) определяется первым способом (т.е. при iniapr = 1), то для получения максимального порядка точности приближенного решения в конце элементарного сегмента необходимо выполнить в итерационном процессе не менее k итераций; тогда imax≥k. Если начальное приближение коэффициентов Чебышёва функции Φ(α) определяется вторым способом (т.е. при iniapr = 2), то для получения максимального порядка точности приближенного решения необходимо выполнить в итерационном процессе не менее k + 1 итераций; в этом случае imax≥k + 1. Однако в некоторых случаях при втором способе определения начального приближения итерационный процесс может сойтись за значительно меньшее число итераций. Если диаметр разбиения h подобран достаточно малым, то хорошая точность приближенного решения может быть получена и с существенно меньшим числом итераций при любом способе выбора начального приближения. Вообще, число итераций зависит от k и h. С увеличением h число итераций может также возрастать. Если правая часть дифференциального уравнения не зависит от переменной Y, т.е. дифференциальное уравнение имеет вид Y'= F(X), то число итераций можно положить равным 1 при любых h и k, удовлетворяющих описанным выше условиям. В этом случае параметр imax = 1. Как следует из вышеописанного, управлять точностью приближенного решения задачи Коши можно с помощью четырех параметров h, k, imax, iniapr, подбирая для каждой конкретной задачи наиболее подходящий набор их значений. При работе функции значения параметров m, xn, yn, xk, k, iniapr, imax, сохраняются. Значение параметра h сохраняется, если он задан с учетом направления интегрирования, иначе его знак меняется на противоположный. Если после работы функции нет необходимости иметь начальное значение решения, то параметры yn и y при обращении к ней можно совместить. При интегрировании системы уравнений с помощью функции de53d_c используется внешняя структура, элементы которой нельзя портить: struct { double wc1, wc2, wc3; double hd2, hd4, hold; int lasn; }com70d_; |
Вычисления в трех примерах на Си проводились с 15-16 значащими цифрами.
1) Решается задача Коши для системы уравнений
y'1 = -2q (y2-1) + (1 -e1-y1+cos q(2x-1) ) / (x+1), y1(0) = 1+ cos q, 0 ≤ x ≤ xf, xf = 1,
y'2 = 2q (y1-1) + ( 1 -e1-y2+sin q(2x-1) ) / (x+1), y2(0) = 1- sin q, q = 0,5.
Точное решение системы
y1(x) = 1 + cos q(2x - 1), y2(x) = 1 + sin q(2x - 1).
Функции y1(x), y2(x) разлагаются на отрезке [0, 1] в смещенные ряды Чебышёва, коэффициенты которых выражаются через цилиндрические функции:
∞
(7) y1(x) = 1 + J0 (q) + 2 ∑ (-1)i J2i (q) T2i* (x),
i=1
∞
(8) y2(x) = 1 + 2 ∑ (-1)i J2i+1 (q) T2i+1* (x).
i=0
Интервал интегрирования [0, 1] разбивается на два элементарных сегмента длиной h = 0,5. В примере 1 функция обработки результатов ftreat_c выполняет следующие действия: 1) печатает границы каждого элементарного сегмента xi, xe; 2) печатает значение решения задачи Коши в конце каждого элементарного сегмента, т.е. Y(XE); 3) печатает вычисленные на каждом отдельном элементарном сегменте коэффициенты Чебышёва решения и его производной. Приводятся фрагмент вызывающей программы, функция f_c вычисления значений правой части системы, функция ftreat_c обработки результатов. Далее представлены результаты работы функции ftreat_c для каждого из двух элементарных сегментов, значения параметров h, k, iniapr, imax, а также вычисленное функцией de53d_c значение y решения задачи Коши в точке xk = xf, точное значение yt решения в точке xk = xf и абсолютная погрешность dely приближенного решения y. Все перечисленные данные приводятся при двух способах определения начального приближения коэффициентов Чебышёва правой части (производной решения), т.е. при iniapr = 1 и при iniapr = 2.
int main(void)
{
/* Builtin functions */
double cos(double), sin(double);
/* Local variables */
extern int f_c();
static double h__;
static int k, m;
static double q, y[2], h0, xk, xn, yn[2], yt[2], rab[446];
extern int de53d_c( );
static double dely[2];
static int imax, iniap0, iniapr;
extern int ftreat_c();
m = 2;
q = .5;
xn = 0.f;
yn[0] = cos(q) + 1.;
yn[1] = -sin(q) + 1.;
xk = 1.;
h0 = .5;
yt[0] = cos(q * (xk * 2. - 1.)) + 1.;
yt[1] = sin(q * (xk * 2. - 1.)) + 1.;
k = 11;
imax = 13;
for (iniap0 = 1; iniap0 <= 2; ++iniap0) {
iniapr = iniap0;
h__ = h0;
de53d_c((U_fp)f_c, (U_fp)ftreat_c, &m, &xn, yn, &xk, &k, &iniapr, &imax,
&h__, y, rab);
dely[0] = yt[0] - y[0];
dely[1] = yt[1] - y[1];
/* Операторы вывода на печать: h0, k, iniapr, imax, y, yt, dely */
. . . . . . . . . . . . . . . . . . . . . . . . . .
return 0;
} /* main */
/*..........................................................................*/
int f_c(double *x, double *y, double *z__, int *m)
{
/* Builtin functions */
double cos(double), exp(double), sin(double);
/* Local variables */
static double q, r__;
/* Parameter adjustments */
--z__;
--y;
/* Function Body */
q = .5;
r__ = *x * 2. - 1.;
z__[1] = q * -2. * (y[2] - 1.) + (1. - exp(1. - y[1] + cos(q * r__))) / (*
x + 1.);
z__[2] = q * 2. * (y[1] - 1.) + (1. - exp(1. - y[2] + sin(q * r__))) / (*
x + 1.);
return 0;
} /* f_c */
/*..........................................................................*/
int ftreat_c(int *s, double *xi, double *xe, double *y, double *ay,
double *ady, int *m, int *kp1, int *kp2)
{
/* System generated locals */
int ay_dim1, ay_offset, ady_dim1, ady_offset, y_dim1, i__1, i__2,
i__3;
/* Local variables */
static int j, l;
/* Parameter adjustments */
y_dim1 = *m;
--y;
ady_dim1 = *m;
ady_offset = 1 + ady_dim1;
ady -= ady_offset;
ay_dim1 = *m;
ay_offset = 1 + ay_dim1;
ay -= ay_offset;
/* Операторы вывода на печать: xi, xe, y, ay, ady */
. . . . . . . . . . . . . . . . . . . . . . . . . .
return 0;
} /* ftreat_c */
Результаты:
---------------------------------------------------------
COEFFICIENTS AY AND ADY ON 1 SEGMENT:
xi= 0.00000000000000e+000 xe=5.00000000000000e-001
SOLUTION Y AT THE END OF THE SEGMENT [XI,XE], i.e. y[xe]
Y=
2.00000000000000e+000 1.00000000000000e+000
Chebyshev coefficients for 1 component
-------------------------------------------------------
for y for y'
0 3.90766440054603e+000 4.87106693080399e-001
1 6.13690356801086e-002 -2.40340620085586e-001
2 -1.50605601386582e-002 -3.84559236046984e-003
3 -1.60442087410513e-004 6.28342132945954e-004
4 1.96510520421465e-005 5.01773738245725e-006
5 1.25508828074456e-007 -4.91532402733373e-007
6 -1.02440704881216e-008 -2.61574052097324e-009
7 -4.67226947263059e-011 1.82980696465326e-010
8 2.85969964497412e-012 7.30383699891089e-013
9 1.01480150204011e-014 -4.00808130174984e-014
10 -5.00565809234601e-016 -2.73381577792220e-016
11 -3.10660883854795e-018 -3.55482787303685e-017
12 -3.70294570108005e-019
Chebyshev coefficients for 2 component
for y for y'
0 1.51289330691960e+000 1.90766440054603e+000
1 2.40340620085586e-001 6.13690356801088e-002
2 3.84559236046995e-003 -1.50605601386584e-002
3 -6.28342132945851e-004 -1.60442087410362e-004
4 -5.01773738245519e-006 1.96510520420718e-005
5 4.91532402808844e-007 1.25508828203872e-007
6 2.61574065406908e-009 -1.02440702819581e-008
7 -1.82980898893858e-010 -4.67231914440323e-011
8 -7.30192778241261e-013 2.86005609793178e-012
9 3.97340507167311e-014 9.14636340840561e-015
10 1.12470474692436e-016 -7.95553672851973e-016
11 -9.04038264604515e-018 1.48725433010699e-016
12 1.54922326052812e-018
COEFFICIENTS AY AND ADY ON 2 SEGMENT:
xi= 5.00000000000000e-001 xe=1.00000000000000e+000
SOLUTION Y AT THE END OF THE SEGMENT [XI,XE], i.e. y[xe]
Y=
1.87758256189037e+000 1.47942553860420e+000
Chebyshev coefficients for 1 component
for y for y'
0 3.90766440054603e+000 -4.87106693080399e-001
1 -6.13690356801086e-002 -2.40340620085586e-001
2 -1.50605601386582e-002 3.84559236046996e-003
3 1.60442087410518e-004 6.28342132945810e-004
4 1.96510520421438e-005 -5.01773738247782e-006
5 -1.25508828079410e-007 -4.91532402791218e-007
6 -1.02440704953297e-008 2.61574069858536e-009
7 4.67226941925040e-011 1.82980984607816e-010
8 2.85969623678284e-012 -7.30176194864496e-013
9 -1.01429241017564e-014 -3.95745462861571e-014
10 -4.94180644706201e-016 1.14340461964586e-016
11 1.29932343141575e-018 -4.00947096609962e-017
12 -4.17653225635377e-019
Chebyshev coefficients for 2 component
-------------------------------------------------------
for y for y'
0 2.48710669308040e+000 1.90766440054603e+000
1 2.40340620085586e-001 -6.13690356801084e-002
2 -3.84559236046994e-003 -1.50605601386584e-002
3 -6.28342132945852e-004 1.60442087410646e-004
4 5.01773738245700e-006 1.96510520420564e-005
5 4.91532402806672e-007 -1.25508827977917e-007
6 -2.61574062872712e-009 -1.02440702104685e-008
7 -1.82980898204373e-010 4.67222009846902e-011
8 7.30209069543573e-013 2.86008897636266e-012
9 3.97349893233514e-014 -1.11794660985140e-014
10 -1.40950094071378e-016 -8.30254918635087e-016
11 -9.43471498448963e-018 9.65414271962561e-017
12 1.00563986662767e-018
-------------------------------------------------------
h0=5.00000000000000e-001 k= 11 iniapr= 1 imax= 13
y[0]= 1.87758256189037e+000 y[1]=1.47942553860420e+000
yt[0]=1.87758256189037e+000 yt[1]= 1.47942553860420e+000
dely
0.00000000000000e+000 0.00000000000000e+000
******************************************************
COEFFICIENTS AY AND ADY ON 1 SEGMENT:
xi= 0.00000000000000e+000 xe=5.00000000000000e-001
SOLUTION Y AT THE END OF THE SEGMENT [XI,XE], i.e. y[xe]
Y=
2.00000000000000e+000 1.00000000000000e+000
Chebyshev coefficients for 1 component
-------------------------------------------------------
for y for y'
0 3.90766440054603e+000 4.87106693080399e-001
1 6.13690356801086e-002 -2.40340620085586e-001
2 -1.50605601386582e-002 -3.84559236046984e-003
3 -1.60442087410513e-004 6.28342132945954e-004
4 1.96510520421465e-005 5.01773738245725e-006
5 1.25508828074456e-007 -4.91532402733373e-007
6 -1.02440704881216e-008 -2.61574052097324e-009
7 -4.67226947263059e-011 1.82980696465326e-010
8 2.85969964497412e-012 7.30383699891089e-013
9 1.01480150204011e-014 -4.00808130174984e-014
10 -5.00565809234601e-016 -2.73381577792220e-016
11 -3.10660883854795e-018 -3.55482787303685e-017
12 -3.70294570108005e-019
Chebyshev coefficients for 2 component
-------------------------------------------------------
for y for y'
0 1.51289330691960e+000 1.90766440054603e+000
1 2.40340620085586e-001 6.13690356801088e-002
2 3.84559236046995e-003 -1.50605601386584e-002
3 -6.28342132945851e-004 -1.60442087410362e-004
4 -5.01773738245519e-006 1.96510520420718e-005
5 4.91532402808844e-007 1.25508828203872e-007
6 2.61574065406908e-009 -1.02440702819581e-008
7 -1.82980898893858e-010 -4.67231914440323e-011
8 -7.30192778241261e-013 2.86005609793178e-012
9 3.97340507167311e-014 9.14636340840561e-015
10 1.12470474692436e-016 -7.95553672851973e-016
11 -9.04038264604515e-018 1.48725433010699e-016
12 1.54922326052812e-018
COEFFICIENTS AY AND ADY ON 2 SEGMENT:
xi= 5.00000000000000e-001 xe=1.00000000000000e+000
SOLUTION Y AT THE END OF THE SEGMENT [XI,XE], i.e. y[xe]
Y=
1.87758256189037e+000 1.47942553860420e+000
Chebyshev coefficients for 1 component
-------------------------------------------------------
for y for y'
0 3.90766440054603e+000 -4.87106693080399e-001
1 -6.13690356801086e-002 -2.40340620085586e-001
2 -1.50605601386582e-002 3.84559236046991e-003
3 1.60442087410517e-004 6.28342132945817e-004
4 1.96510520421444e-005 -5.01773738250524e-006
5 -1.25508828079679e-007 -4.91532402802669e-007
6 -1.02440704950977e-008 2.61574068191242e-009
7 4.67226936695002e-011 1.82980962021866e-010
8 2.85969642735029e-012 -7.30163579597710e-013
9 -1.01431737831314e-014 -3.96093285524830e-014
10 -4.94137214012519e-016 1.44932787748089e-016
11 1.64696349713737e-018 -7.83514314814895e-017
12 -8.16160744598849e-019
Chebyshev coefficients for 2 component
-------------------------------------------------------
for y for y'
0 2.48710669308040e+000 1.90766440054603e+000
1 2.40340620085586e-001 -6.13690356801084e-002
2 -3.84559236046994e-003 -1.50605601386584e-002
3 -6.28342132945851e-004 1.60442087410608e-004
4 5.01773738245640e-006 1.96510520420275e-005
5 4.91532402806431e-007 -1.25508827997223e-007
6 -2.61574062872712e-009 -1.02440702297741e-008
7 -1.82980898031941e-010 4.67221816791152e-011
8 7.30208767893965e-013 2.86006001461213e-012
9 3.97348552098014e-014 -1.11794660985140e-014
10 -1.41191413758051e-016 -8.49560493568907e-016
11 -9.65409651782849e-018 1.15847002130076e-016
12 1.20673960552163e-018
-------------------------------------------------------
h0=5.00000000000000e-001 k= 11 iniapr= 2 imax= 13
y[0]= 1.87758256189037e+000 y[1]=1.47942553860420e+000
yt[0]=1.87758256189037e+000 yt[1]= 1.47942553860420e+000
dely
0.00000000000000e+000 0.00000000000000e+000
******************************************************
2) Решается задача Коши из примера 1 также с разбиением интервала интегрирования [0, 1] на два элементарных сегмента длиной h = 0,5. В этом примере компоненты решения y1(x) и y2(x) рассматриваются как координаты, а компоненты производной решения y'1(x) и y'2(x) отождествляются с компонентами скорости. Для данного примера функция обработки результатов ftreat_c устроена следующим образом. Коэффициенты Чебышёва для каждой координаты запоминаются (сохраняются) в отдельном массиве. Всего таких массивов четыре (по два массива на каждом сегменте). Коэффициенты Чебышёва для каждой компоненты скорости также запоминаются в отдельных массивах и таких массивов тоже четыре (по два массива на каждом сегменте). Эти массивы помещаются во внешнюю структуру first_ для первого сегмента и во внешнюю структуру second_ для второго сегмента. В этих же структурах размещаются значения координат в конце каждого сегмента (интервала) и границы каждого сегмента. Поскольку структуры first_ и second_ внешние, то на выходе из функции de53d_c вычисленные в ней коэффициенты Чебышёва будут доступны также в главной функции. Приводятся фрагмент вызывающей функции, функция f_c вычисления значений правой части системы, функция ftreat_c обработки результатов. Далее представлены значения параметров h, k, iniapr, imax, приближенное решение y, вычисленное в конце интервала xk = xf, точное решение yt (в точке xf) и абсолютная погрешность приближенного значения y. Затем даются те же результаты, что и в примере 1, а именно: границы элементарных сегментов, значения координат в конце каждого сегмента, коэффициенты Чебышёва для координат и компонент скорости на каждом сегменте при значении параметра iniapr = 1.
/* Common Block Declarations */
struct {
double ay1[13], ay2[13], av1[12], av2[12], ya[2], ina[2];
} first_;
struct {
double by1[13], by2[13], bv1[12], bv2[12], yb[2], inb[2];
} second_;
/*.............................................................................*/
int main(void)
{
/* System generated locals */
int i__1, i__2;
/* Builtin functions */
double cos(double), sin(double);
/* Local variables */
extern int f_c();
static double h__;
static int j, k, m;
static double q, y[2], h0, xk, xn, yn[2], yt[2];
static int kp1, kp2;
static double rab[446];
extern int de53d_c( );
static double dely[2];
static int imax, iniapr;
extern int ftreat_c();
/* AY1 - CHEBYSHEV COEFFICIENTS FOR 1 COORDINATE AT 1 INTERVAL */
/* AY2 - CHEBYSHEV COEFFICIENTS FOR 2 COORDINATE AT 1 INTERVAL */
/* AV1 - CHEBYSHEV COEFFICIENTS FOR 1 VELOCITY COMPONENT AT 1 INTERVAL */
/* AV2 - CHEBYSHEV COEFFICIENTS FOR 2 VELOCITY COMPONENT AT 1 INTERVAL */
/* BY1 - CHEBYSHEV COEFFICIENTS FOR 1 COORDINATE AT 2 INTERVAL */
/* BY2 - CHEBYSHEV COEFFICIENTS FOR 2 COORDINATE AT 2 INTERVAL */
/* BV1 - CHEBYSHEV COEFFICIENTS FOR 1 VELOCITY COMPONENT AT 2 INTERVAL */
/* BV2 - CHEBYSHEV COEFFICIENTS FOR 2 VELOCITY COMPONENT AT 2 INTERVAL */
/* YA - COORDINATES AT THE END OF 1 INTERVAL */
/* YB - COORDINATES AT THE END OF 2 INTERVAL */
/* INA - BOUNDARIES OF 1 INTERVAL */
/* INB - BOUNDARIES OF 2 INTERVAL */
m = 2;
q = .5;
xn = 0.f;
yn[0] = cos(q) + 1.;
yn[1] = -sin(q) + 1.;
xk = 1.;
h0 = .5;
yt[0] = cos(q * (xk * 2. - 1.)) + 1.;
yt[1] = sin(q * (xk * 2. - 1.)) + 1.;
k = 11;
imax = 13;
iniapr = 1;
h__ = h0;
de53d_c((U_fp)f_c, (U_fp)ftreat_c, &m, &xn, yn, &xk, &k, &iniapr, &imax,
&h__, y, rab);
dely[0] = yt[0] - y[0];
dely[1] = yt[1] - y[1];
/* Операторы вывода на печать: h0, k, iniapr, imax, y, yt, dely, */
/* ina,ya,ay1,av1,ay2,av2,inb.yb,by1,bv1,by2,bv2 */
. . . . . . . . . . . . . . . . . . . . . . . . . .
return 0;
} /* main */
/*.............................................................................*/
int f_c(double *x, double *y, double *z__, int *m)
{
/* Builtin functions */
double cos(double), exp(double), sin(double);
/* Local variables */
static double q, r__;
/* Parameter adjustments */
--z__;
--y;
/* Function Body */
q = .5;
r__ = *x * 2. - 1.;
z__[1] = q * -2. * (y[2] - 1.) + (1. - exp(1. - y[1] + cos(q * r__))) /
(*x + 1.);
z__[2] = q * 2. * (y[1] - 1.) + (1. - exp(1. - y[2] + sin(q * r__))) /
(*x + 1.);
return 0;
} /* f_c */
/*.............................................................................*/
int ftreat_c(int *s, double *xi, double *xe, double *y,
double *ay, double *ady, int *m, int *kp1, int *kp2)
{
/* System generated locals */
int ay_dim1, ay_offset, ady_dim1, ady_offset, i__1;
/* Local variables */
static int i__;
/* AY1 - CHEBYSHEV COEFFICIENTS FOR 1 COORDINATE AT 1 INTERVAL */
/* AY2 - CHEBYSHEV COEFFICIENTS FOR 2 COORDINATE AT 1 INTERVAL */
/* AV1 - CHEBYSHEV COEFFICIENTS FOR 1 VELOCITY COMPONENT AT 1 INTERVAL */
/* AV2 - CHEBYSHEV COEFFICIENTS FOR 2 VELOCITY COMPONENT AT 1 INTERVAL */
/* BY1 - CHEBYSHEV COEFFICIENTS FOR 1 COORDINATE AT 2 INTERVAL */
/* BY2 - CHEBYSHEV COEFFICIENTS FOR 2 COORDINATE AT 2 INTERVAL */
/* BV1 - CHEBYSHEV COEFFICIENTS FOR 1 VELOCITY COMPONENT AT 2 INTERVAL */
/* BV2 - CHEBYSHEV COEFFICIENTS FOR 2 VELOCITY COMPONENT AT 2 INTERVAL */
/* YA - COORDINATES AT THE END OF 1 INTERVAL */
/* YB - COORDINATES AT THE END OF 2 INTERVAL */
/* INA - BOUNDARIES OF 1 INTERVAL */
/* INB - BOUNDARIES OF 2 INTERVAL */
/* Parameter adjustments */
--y;
ady_dim1 = *m;
ady_offset = 1 + ady_dim1;
ady -= ady_offset;
ay_dim1 = *m;
ay_offset = 1 + ay_dim1;
ay -= ay_offset;
/* Function Body */
switch (*s) {
case 1: goto L10;
case 2: goto L30;
}
/* --------SAVING CHEBYSHEV COEFFICIENTS FOR COORDINATES AND */
/* -------- VELOCITY COMPONENTS AT THE FIRST INTERVAL */
L10:
i__1 = *kp2;
for (i__ = 1; i__ <= i__1; ++i__) {
first_1.ay1[i__ - 1] = ay[i__ * ay_dim1 + 1];
/* L15: */
first_1.ay2[i__ - 1] = ay[i__ * ay_dim1 + 2];
}
i__1 = *kp1;
for (i__ = 1; i__ <= i__1; ++i__) {
first_1.av1[i__ - 1] = ady[i__ * ady_dim1 + 1];
/* L20: */
first_1.av2[i__ - 1] = ady[i__ * ady_dim1 + 2];
}
/* --------SAVING COORDINATES AT THE END OF 1 INTERVAL */
first_1.ya[0] = y[1];
first_1.ya[1] = y[2];
/* --------SAVING BOUNDARIES OF 1 INTERVAL */
first_1.ina[0] = *xi;
first_1.ina[1] = *xe;
return 0;
/* --------SAVING CHEBYSHEV COEFFICIENTS FOR COORDINATES AND */
/* -------- VELOCITY COMPONENTS AT THE SECOND INTERVAL */
L30:
i__1 = *kp2;
for (i__ = 1; i__ <= i__1; ++i__) {
second_1.by1[i__ - 1] = ay[i__ * ay_dim1 + 1];
/* L35: */
second_1.by2[i__ - 1] = ay[i__ * ay_dim1 + 2];
}
i__1 = *kp1;
for (i__ = 1; i__ <= i__1; ++i__) {
second_1.bv1[i__ - 1] = ady[i__ * ady_dim1 + 1];
/* L40: */
second_1.bv2[i__ - 1] = ady[i__ * ady_dim1 + 2];
}
/* --------SAVING COORDINATES AT THE END OF 2 INTERVAL */
second_1.yb[0] = y[1];
second_1.yb[1] = y[2];
/* --------SAVING BOUNDARIES OF 2 INTERVAL */
second_1.inb[0] = *xi;
second_1.inb[1] = *xe;
return 0;
} /* ftreat_c */
Результаты:
---------------------------------------------------------
h0=5.00000000000000e-001 k= 11 iniapr= 1 imax= 13
y[0]= 1.87758256189037e+000 y[1]=1.47942553860420e+000
yt[0]=1.87758256189037e+000 yt[1]= 1.47942553860420e+000
dely
0.00000000000000e+000 0.00000000000000e+000
******************************************************
CHEBYSHEV COEFFICIENTS AT 1 SEGMENT
ina[0]= 0.00000000000000e+000 ina[1]=5.00000000000000e-001
COORDINATES AT THE END OF 1 SEGMENT
ya[0]= 2.00000000000000e+000 ya[1]=1.00000000000000e+000
Coefficients for 1 coordinate and 1 velocity component
-------------------------------------------------------
coordinate velocity
0 3.90766440054603e+000 4.87106693080399e-001
1 6.13690356801086e-002 -2.40340620085586e-001
2 -1.50605601386582e-002 -3.84559236046984e-003
3 -1.60442087410513e-004 6.28342132945954e-004
4 1.96510520421465e-005 5.01773738245725e-006
5 1.25508828074456e-007 -4.91532402733373e-007
6 -1.02440704881216e-008 -2.61574052097324e-009
7 -4.67226947263059e-011 1.82980696465326e-010
8 2.85969964497412e-012 7.30383699891089e-013
9 1.01480150204011e-014 -4.00808130174984e-014
10 -5.00565809234601e-016 -2.73381577792220e-016
11 -3.10660883854795e-018 -3.55482787303685e-017
12 -3.70294570108005e-019
Coefficients for 2 coordinate and 2 velocity component
-------------------------------------------------------
coordinate velocity
0 1.51289330691960e+000 1.90766440054603e+000
1 2.40340620085586e-001 6.13690356801088e-002
2 3.84559236046995e-003 -1.50605601386584e-002
3 -6.28342132945851e-004 -1.60442087410362e-004
4 -5.01773738245519e-006 1.96510520420718e-005
5 4.91532402808844e-007 1.25508828203872e-007
6 2.61574065406908e-009 -1.02440702819581e-008
7 -1.82980898893858e-010 -4.67231914440323e-011
8 -7.30192778241261e-013 2.86005609793178e-012
9 3.97340507167311e-014 9.14636340840561e-015
10 1.12470474692436e-016 -7.95553672851973e-016
11 -9.04038264604515e-018 1.48725433010699e-016
12 1.54922326052812e-018
******************************************************
CHEBYSHEV COEFFICIENTS AT 2 SEGMENT
inb[0]= 5.00000000000000e-001 inb[1]=1.00000000000000e+000
COORDINATES AT THE END OF 2 SEGMENT:
yb[0]= 1.87758256189037e+000 yb[1]=1.47942553860420e+000
Coefficients for 1 coordinate and 1 velocity component
-------------------------------------------------------
coordinate velocity
0 3.90766440054603e+000 -4.87106693080399e-001
1 -6.13690356801086e-002 -2.40340620085586e-001
2 -1.50605601386582e-002 3.84559236046996e-003
3 1.60442087410518e-004 6.28342132945810e-004
4 1.96510520421438e-005 -5.01773738247782e-006
5 -1.25508828079410e-007 -4.91532402791218e-007
6 -1.02440704953297e-008 2.61574069858536e-009
7 4.67226941925040e-011 1.82980984607816e-010
8 2.85969623678284e-012 -7.30176194864496e-013
9 -1.01429241017564e-014 -3.95745462861571e-014
10 -4.94180644706201e-016 1.14340461964586e-016
11 1.29932343141575e-018 -4.00947096609962e-017
12 -4.17653225635377e-019
Coefficients for 2 coordinate and 2 velocity component
-------------------------------------------------------
coordinate velocity
0 2.48710669308040e+000 1.90766440054603e+000
1 2.40340620085586e-001 -6.13690356801084e-002
2 -3.84559236046994e-003 -1.50605601386584e-002
3 -6.28342132945852e-004 1.60442087410646e-004
4 5.01773738245700e-006 1.96510520420564e-005
5 4.91532402806672e-007 -1.25508827977917e-007
6 -2.61574062872712e-009 -1.02440702104685e-008
7 -1.82980898204373e-010 4.67222009846902e-011
8 7.30209069543573e-013 2.86008897636266e-012
9 3.97349893233514e-014 -1.11794660985140e-014
10 -1.40950094071378e-016 -8.30254918635087e-016
11 -9.43471498448963e-018 9.65414271962561e-017
12 1.00563986662767e-018
******************************************************
3) Решается задача Коши из примера 1 также с разбиением интервала интегрирования [0, 1] на два элементарных
сегмента длиной h = 0,5. Система уравнений в примере 1 описывает в трехмерном пространстве переменных
x, y1, y2 движение точки по цилиндрической поверхности. Ось кругового цилиндра параллельна
оси x и проходит через точку (0, 1, 1). Фазовая траектория этой системы на плоскости
y1, y2 -- окружность с центром в точке (1, 1) радиуса 1. При изменении x
точка (y1(x), y2(x)) совершает движение по дуге окружности. В примере 3
функция обработки результатов выполняет следующие действия.
а) По вычисленным функцией de53d_c коэффициентам Чебышёва ay находится решение
системы y1(x), y2(x) для нескольких равноотстоящих значений аргумента x, принадлежащих
элементарному сегменту [xi, xe] (число этих значений задается переменной ni).
б) Определяется расстояние от точки (y1(x), y2(x)) до точки (1, 1); все расстояния
запоминаются в массиве r.
в) Вычисляется возможное отклонение этого расстояния от постоянного значения (равного 1). Это отклонение может
быть следствием приближенного интегрирования системы уравнений из-за замены точного решения задачи частичной суммой
ряда с приближенными коэффициентами. Все отклонения запоминаются в массиве deltr.
Массивы r и deltr помещаются во внешнюю структуру deltr_ вместе с переменной ni, представляющей число точек x,
взятых на элементарном сегменте. Поскольку структура deltr_ внешняя, то на выходе из функции de53d_c вычисленные в функции ftreat_c массивы r и deltr
будут доступны также и в главной функции. Приводятся фрагмент вызывающей функции, функции f_c вычисления
значений правой части системы, функция ftreat_c обработки результатов. Далее представлены значения параметров
h, k, iniapr, imax, приближенное решение y, вычисленное в конце интервала xk = xf, точное
решение yt (в точке xf) и абсолютная погрешность приближенного значения y. Затем даются расстояния от
точек (y1(x), y2(x)) до точки (1, 1) на фазовой плоскости и отклонения этих расстояний
от постоянного значения, равного 1, для каждого элементарного сегмента. Эти данные приводятся для двух значений параметра
iniapr (1 и 2).
/* Common Block Declarations */
struct {
double r__[32] /* was [16][2] */, deltr[32] /* was [16][2] */;
int ni;
} deltr_;
/*.............................................................................*/
int main(void)
{
/* System generated locals */
int i__1;
/* Builtin functions */
double cos(double), sin(double);
/* Local variables */
extern int f_c();
static double h__;
static int j, k, m;
static double q;
static int s;
static double y[2], h0, xk, xn, yn[2], yt[2], rab[446];
extern int de53d_c( );
static double dely[2];
static int imax, iniap0, iniapr;
extern int ftreat_c();
deltr_1.ni = 16;
m = 2;
q = .5;
xn = 0.f;
yn[0] = cos(q) + 1.;
yn[1] = -sin(q) + 1.;
xk = 1.;
h0 = .5;
yt[0] = cos(q * (xk * 2. - 1.)) + 1.;
yt[1] = sin(q * (xk * 2. - 1.)) + 1.;
k = 11;
imax = 13;
for (iniap0 = 1; iniap0 <= 2; ++iniap0) {
iniapr = iniap0;
h__ = h0;
de53d_c((U_fp)f_c, (U_fp)ftreat_c, &m, &xn, yn, &xk, &k, &iniapr, &imax,
&h__, y, rab);
dely[0] = yt[0] - y[0];
dely[1] = yt[1] - y[1];
/* Операторы вывода на печать: h0, k, iniapr, imax, y, yt, dely, deltr */
. . . . . . . . . . . . . . . . . . . . . . . . . .
return 0;
} /* main */
/*.......................................................................*/
int f_c(double *x, double *y, double *z__, int *m)
{
/* Builtin functions */
double cos(double), exp(double), sin(double);
/* Local variables */
static double q, r__;
/* Parameter adjustments */
--z__;
--y;
/* Function Body */
q = .5;
r__ = *x * 2. - 1.;
z__[1] = q * -2. * (y[2] - 1.) + (1. - exp(1. - y[1] + cos(q * r__))) / (*
x + 1.);
z__[2] = q * 2. * (y[1] - 1.) + (1. - exp(1. - y[2] + sin(q * r__))) / (*
x + 1.);
return 0;
} /* f_c */
/*.............................................................................*/
int ftreat_c(int *s, double *xi, double *xe, double *y,
double *ay, double *ady, int *m, int *kp1, int *kp2)
{
/* System generated locals */
int ay_dim1, ay_offset, ady_dim1, ady_offset, i__1;
double d__1, d__2;
/* Builtin functions */
double sqrt(double);
/* Local variables */
static int j;
static double x, bk[4] /* was [2][2] */, hi, yi[2];
extern int de70dc_c( );
/* Parameter adjustments */
ady_dim1 = *m;
ady_offset = 1 + ady_dim1;
ady -= ady_offset;
ay_dim1 = *m;
ay_offset = 1 + ay_dim1;
ay -= ay_offset;
/* Function Body */
hi = (*xe - *xi) / (double) deltr_1.ni;
x = *xi;
i__1 = deltr_1.ni;
for (j = 1; j <= i__1; ++j) {
x += hi;
/* ------CALCULATION OF SOLUTION YI AT POINT X (BY ALGORITHM CLENSHAW) */
de70dc_c(m, kp1, &x, &ay[ay_offset], yi, bk);
/* Computing 2nd power */
d__1 = yi[0] - 1.;
/* Computing 2nd power */
d__2 = yi[1] - 1.;
deltr_1.r__[j + (*s << 4) - 17] = sqrt(d__1 * d__1 + d__2 * d__2);
deltr_1.deltr[j + (*s << 4) - 17] = 1. - deltr_1.r__[j + (*s << 4) -
17];
/* L10: */
}
return 0;
} /* ftreat_c */
Результаты:
---------------------------------------------------------
h0=5.00000000000000e-001 k= 11 iniapr= 1 imax= 13
y[0]= 1.87758256189037e+000 y[1]=1.47942553860420e+000
yt[0]=1.87758256189037e+000 yt[1]= 1.47942553860420e+000
dely
0.00000000000000e+000 0.00000000000000e+000
-------------------------------------------------------
Phase motion for 1 segment
-------------------------------------------------------
distance error
1 1.00000000000000e+000 0.00000000000000e+000
2 1.00000000000000e+000 0.00000000000000e+000
3 1.00000000000000e+000 0.00000000000000e+000
4 1.00000000000000e+000 0.00000000000000e+000
5 1.00000000000000e+000 0.00000000000000e+000
6 1.00000000000000e+000 0.00000000000000e+000
7 1.00000000000000e+000 0.00000000000000e+000
8 1.00000000000000e+000 0.00000000000000e+000
9 1.00000000000000e+000 0.00000000000000e+000
10 1.00000000000000e+000 0.00000000000000e+000
11 1.00000000000000e+000 1.11022302462516e-016
12 1.00000000000000e+000 0.00000000000000e+000
13 1.00000000000000e+000 0.00000000000000e+000
14 1.00000000000000e+000 0.00000000000000e+000
15 1.00000000000000e+000 0.00000000000000e+000
16 1.00000000000000e+000 0.00000000000000e+000
Phase motion for 2 segment
-------------------------------------------------------
distance error
1 1.00000000000000e+000 0.00000000000000e+000
2 1.00000000000000e+000 0.00000000000000e+000
3 1.00000000000000e+000 0.00000000000000e+000
4 1.00000000000000e+000 0.00000000000000e+000
5 1.00000000000000e+000 0.00000000000000e+000
6 1.00000000000000e+000 0.00000000000000e+000
7 1.00000000000000e+000 0.00000000000000e+000
8 1.00000000000000e+000 0.00000000000000e+000
9 1.00000000000000e+000 0.00000000000000e+000
10 1.00000000000000e+000 0.00000000000000e+000
11 1.00000000000000e+000 0.00000000000000e+000
12 1.00000000000000e+000 0.00000000000000e+000
13 1.00000000000000e+000 1.11022302462516e-016
14 1.00000000000000e+000 0.00000000000000e+000
15 1.00000000000000e+000 0.00000000000000e+000
16 1.00000000000000e+000 0.00000000000000e+000
******************************************************
h0=5.00000000000000e-001 k= 11 iniapr= 2 imax= 13
y[0]= 1.87758256189037e+000 y[1]=1.47942553860420e+000
yt[0]=1.87758256189037e+000 yt[1]= 1.47942553860420e+000
dely
0.00000000000000e+000 0.00000000000000e+000
Phase motion for 1 segment
-------------------------------------------------------
distance error
1 1.00000000000000e+000 0.00000000000000e+000
2 1.00000000000000e+000 0.00000000000000e+000
3 1.00000000000000e+000 0.00000000000000e+000
4 1.00000000000000e+000 0.00000000000000e+000
5 1.00000000000000e+000 0.00000000000000e+000
6 1.00000000000000e+000 0.00000000000000e+000
7 1.00000000000000e+000 0.00000000000000e+000
8 1.00000000000000e+000 0.00000000000000e+000
9 1.00000000000000e+000 0.00000000000000e+000
10 1.00000000000000e+000 0.00000000000000e+000
11 1.00000000000000e+000 1.11022302462516e-016
12 1.00000000000000e+000 0.00000000000000e+000
13 1.00000000000000e+000 0.00000000000000e+000
14 1.00000000000000e+000 0.00000000000000e+000
15 1.00000000000000e+000 0.00000000000000e+000
16 1.00000000000000e+000 0.00000000000000e+000
Phase motion for 2 segment
-------------------------------------------------------
distance error
1 1.00000000000000e+000 0.00000000000000e+000
2 1.00000000000000e+000 0.00000000000000e+000
3 1.00000000000000e+000 0.00000000000000e+000
4 1.00000000000000e+000 0.00000000000000e+000
5 1.00000000000000e+000 0.00000000000000e+000
6 1.00000000000000e+000 0.00000000000000e+000
7 1.00000000000000e+000 0.00000000000000e+000
8 1.00000000000000e+000 0.00000000000000e+000
9 1.00000000000000e+000 0.00000000000000e+000
10 1.00000000000000e+000 0.00000000000000e+000
11 1.00000000000000e+000 0.00000000000000e+000
12 1.00000000000000e+000 0.00000000000000e+000
13 1.00000000000000e+000 1.11022302462516e-016
14 1.00000000000000e+000 0.00000000000000e+000
15 1.00000000000000e+000 0.00000000000000e+000
16 1.00000000000000e+000 0.00000000000000e+000
******************************************************