|
Текст подпрограммы и версий (Си) de45d_c.zip |
Тексты тестовых примеров (Си) tde45d1_c.zip , tde45d2_c.zip , tde45d3_c.zip , tde45d4_c.zip |
Вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом рядов Чебышёва.
Решается задача Коши для системы M обыкновенных дифференциальных уравнений пepвoгo порядка
Y ' = F (X, Y) ,
Y = ( y1, ... , yM ) ,
F = ( f1 (X, y1, ... , yM), ... , fM (X, y1, ... , yM) )
с начальными условиями, заданными в точке XN :
Y(XN) = YN , YN = ( y10, ... , yM0 ) ,
методом рядов Чебышёва. Решение вычисляется в одной точке XK, которая является концом интервала интегрирования. Предполагается, что правая часть системы имеет непрерывные ограниченные частные производные по переменным X, Y. Тогда решение системы разлагается на промежутке интегрирования в равномерно сходящийся ряд по смещенным многочленам Чебышёва первого рода.
Интервал интегрирования разбивается на элементарные сегменты длины H: [xs, xs + H], x0 = XN, s = 0,1,...,[(XK - XN) / H] - 1 (квадратные скобки означают целую часть числа). Если длина интервала интегрирования не является целым кратным H, то последний элементарный сегмент считается нестандартным; в этом случае s принимает значения s = 0,1,...,[(XK - XN)/H]. На каждом элементарном сегменте решение исходной задачи Коши приближенно представляется в виде (K + 1) - й частичной суммы смещенного ряда Чебышёва
K+1
Y(xs + αH) ≈ ∑ ai*[Y] Ti*(α) , 0 ≤ α ≤ 1 ,
i=0
здесь ai*[Y] - коэффициенты ряда, Ti*(α) - смещенные многочлены Чебышёва первого рода на [0,1]. Значения H и K задаются пользователем при обращении к подпрограмме. Они выбираются таким образом, чтобы на каждом элементарном сегменте [xs, xs + H] ряд Чебышёва для решения задачи Коши (и его производной) являлся быстросходящимся рядом, а его сумма хорошо приближалась частичной суммой (K + 1) - го порядка (для производной - K - го порядка). Решение в конце каждого сегмента вычисляется по формуле
K+1
Y(xs + H) = Y(xs+1) ≈ ∑ ai*[Y] ,
i=0
В качестве решения в конце интервала интегрирования принимается значение решения в конце последнего элементарного сегмента.
Коэффициенты 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, 24 - 30, 2012.
int de45d_c (U_fp f, int *m, double *xn, double *yn,
double *xk, int *k, int *iniapr, int *imax,
double *h, double *y, double *rab)
Параметры
| f - |
имя функции вычисления значений правой
части дифференциального уравнения. Первый
оператор функции должен иметь вид: int f_c (double *x, double *y, double *dy, int *m) Здесь: x, y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в dy. В случае системы уравнений, т.е. когда m ≠ 1 , параметры y и dy представляют массивы длины m (тип параметров x, y и dy: с двойной точностью); |
| m - | количество уравнений в системе (тип: целый); |
| xn, yn - | начальные значения аргумента и решения; в случае системы уравнений (т.е. когда m ≠ 1) yn представляет массив длины m (тип: с двойной точностью); |
| 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 при обращении к ней можно совместить. При интегрировании системы уравнений с помощью функции de45d_c используется внешняя структура, элементы которой нельзя портить: struct { double wc1, wc2, wc3; double hd2, hd4, hold; int lasn; }com70d_; |
Вычисления во всех четырех примерах на Си проводились с 15-16 значащими цифрами.
1) y1' = x/y2 , y1(0) = 3 ,
y2' = -x/y1 , y(0) = 1/6, 0 ≤ x ≤ 3√2 .
Точное решение системы:
y1(x) = 3*ex*x , y2(x) = (1/6)*e-x*x .
Решение вычисляется в точке xk = 3√2. Приводятся функция вычисления значений правой части и фрагмент вызывающей функции, а также результаты счета, включая абсолютную dely и относительную relerr погрешности приближенного решения,точные значения решения yt, значения параметров h, k, imax, iniapr при двух способах определения начального приближения коэффициентов Чебышёва производной решения.
int main(void)
{
/* Builtin functions */
double sqrt(double), exp(double);
/* Local variables */
extern int f_c();
static double h__;
static int k, m;
static double y[2], h0, xk, xn, yn[2], yt[2], rab[1692];
extern int de45d_c( );
static double dely[2];
static int imax, iniapr;
static double relerr[2];
m = 2;
xn = 0.0;
yn[0] = 3.;
yn[1] = .16666666666666666;
xk = sqrt(2.) * 3.;
yt[0] = exp(xk * xk) * 3.;
yt[1] = exp(-xk * xk) * .16666666666666666;
h0 = .1;
k = 25;
imax = 31;
for (iniapr = 1; iniapr <= 2; ++iniapr) {
h__ = h0;
de45d_c((U_fp)f_c, &m, &xn, yn, &xk, &k, &iniapr, &imax, &h__, y, rab);
dely[0] = yt[0] - y[0];
dely[1] = yt[1] - y[1];
relerr[0] = dely[0] / y[0];
relerr[1] = dely[1] / y[1];
/* Операторы вывода на печать: h0, k, iniapr, imax, y, yt, dely, relerr */
. . . . . . . . . . . . . . . . . . . . . . . . . .
}
return 0;
} /* main */
/*..........................................................................*/
int f_c(double *x, double *y, double *z__, int *m)
{
/* Parameter adjustments */
--z__;
--y;
/* Function Body */
z__[1] = *x / y[2];
z__[2] = -(*x) / y[1];
return 0;
} /* f_c */
Результаты:
--------------
h0= 1.00000000000000e-001 k= 25 iniapr= 1 imax= 31
y[0]= 1.96979907411991e+008 y[1]= 2.53832995745212e-009
yt[0]= 1.96979907411992e+008 yt[1]= 2.53832995745210e-009
dely=
1.54972076416016e-006 -2.23338765389317e-023
relerr=
7.86740528270662e-015 -8.79864986557921e-015
----------------------------------------------------------------------------
h0= 1.00000000000000e-001 k= 25 iniapr= 2 imax= 31
y[0]= 1.96979907411991e+008 y[1]= 2.53832995745212e-009
yt[0]= 1.96979907411992e+008 yt[1]= 2.53832995745210e-009
dely=
1.51991844177246e-006 -1.98523347012727e-023
relerr=
7.71610902726995e-015 -7.82102210273708e-015
----------------------------------------------------------------------------
Видно, что при втором способе выбора начального приближения коэффициентов Чебышёва производной (т.е. при iniapr = 2) приближенное решение вычислено с лучшей точностью, чем при iniapr = 1, т.е. с меньшей абсолютной и относительной погрешностью.
2) Решается задача Коши из примера 1 с другими значениями параметров h, imax: h = 0,4, imax = 39. Так как значение h здесь увеличено по сравнению с примером 1, то и число итераций также возрастает. Приводятся результаты счета, включающие те же данные, что и в примере 1.
h0= 4.00000000000000e-001 k= 25 iniapr= 1 imax= 39 y[0]= 1.96979907411974e+008 y[1]= 2.53832995745239e-009 yt[0]= 1.96979907411992e+008 yt[1]= 2.53832995745210e-009 dely= 1.86860561370850e-005 -2.91581165924942e-022 relerr= 9.48627521587977e-014 -1.14871262133939e-013 ---------------------------------------------------------------------------- h0= 4.00000000000000e-001 k= 25 iniapr= 2 imax= 39 y[0]= 1.96979907411991e+008 y[1]= 2.53832995745211e-009 yt[0]= 1.96979907411992e+008 yt[1]= 2.53832995745210e-009 dely= 9.83476638793945e-007 -1.69572025573371e-023 relerr= 4.99277642940996e-015 -6.68045637942126e-015 ----------------------------------------------------------------------------
Здесь так же, как и в примере 1, при втором способе выбора начального приближения (т.е. при iniapr = 2) приближенное решение вычислено с лучшей точностью, чем при iniapr = 1, т.е. с меньшей абсолютной и относительной погрешностью.
3) Решается система обыкновенных дифференциальных уравнений из примера 1 на отрезке -3√2 ≤ x ≤ 0. Начальные условия задаются в точке t = -3√2:
y1(t) = 3*et*t , y2(t) = (1/6)*e-t*t .
Решение вычисляется в точке xk = 0. Приводятся фрагмент вызывающей функции и результаты счета, включающие те же данные, что и в примере 1.
int main(void)
{
/* Builtin functions */
double sqrt(double), exp(double);
/* Local variables */
extern int f_();
static double h__;
static int k, m;
static double y[2], h0, xk, xn, yn[2], yt[2], rab[1692];
extern int de45d_( );
static double dely[2];
static int imax, iniapr;
static double relerr[2];
m = 2;
xn = sqrt(2.) * -3.;
yn[0] = exp(xn * xn) * 3.;
yn[1] = exp(-xn * xn) * .16666666666666666;
yt[0] = 3.;
yt[1] = .16666666666666666;
xk = 0.0;
h0 = .1;
k = 25;
imax = 30;
for (iniapr = 1; iniapr <= 2; ++iniapr) {
h__ = h0;
de45d_c((U_fp)f_c, &m, &xn, yn, &xk, &k, &iniapr, &imax, &h__, y, rab);
dely[0] = yt[0] - y[0];
dely[1] = yt[1] - y[1];
relerr[0] = dely[0] / y[0];
relerr[1] = dely[1] / y[1];
/* Операторы вывода на печать: h0, k, iniapr, imax, y, yt, dely, relerr */
. . . . . . . . . . . . . . . . . . . . . . . . . .
}
return 0;
} /* main */
/*..........................................................................*/
int f_c(double *x, double *y, double *z__, int *m)
{
/* Parameter adjustments */
--z__;
--y;
/* Function Body */
z__[1] = *x / y[2];
z__[2] = -(*x) / y[1];
return 0;
} /* f_c */
Результаты:
--------------
h0= 1.00000000000000e-001 k= 25 iniapr= 1 imax= 30
y[0]= 3.00000000000001e+000 y[1]= 1.66666666666666e-001
yt[0]= 3.00000000000000e+000 yt[1]= 1.66666666666667e-001
dely=
-1.06581410364015e-014 3.05311331771918e-016
relerr=
-3.55271367880049e-015 1.83186799063151e-015
----------------------------------------------------------------------------
h0= 1.00000000000000e-001 k= 25 iniapr= 2 imax= 30
y[0]= 2.99999999999999e+000 y[1]= 1.66666666666667e-001
yt[0]= 3.00000000000000e+000 yt[1]= 1.66666666666667e-001
dely=
7.10542735760100e-015 -3.33066907387547e-016
relerr=
2.36847578586701e-015 -1.99840144432528e-015
----------------------------------------------------------------------------
4) Решается система обыкновенных дифференциальных уравнений из примера 1, но начальное условие задается в точке t = 3√2:
y1(t) = 3*et*t , y2(t) = (1/6)*e-t*t ,
а решение задачи Коши вычисляется в точке xk = 0. Интегрирование системы выполняется справа налево с отрицательным шагом h = -0.1. Приводятся фрагмент вызывающей функции и результаты счета, включающие такие же данные, как и в примере 1.
int main(void)
{
/* Builtin functions */
double sqrt(double), exp(double);
/* Local variables */
extern int f_c();
static double h__;
static int k, m;
static double y[2], h0, xk, xn, yn[2], yt[2], rab[1692];
extern int de45d_c( );
static double dely[2];
static int imax, iniapr;
static double relerr[2];
m = 2;
xn = sqrt(2.) * 3.;
yn[0] = exp(xn * xn) * 3.;
yn[1] = exp(-xn * xn) * .16666666666666666;
yt[0] = 3.;
yt[1] = .16666666666666666;
xk = 0.0;
h0 = -.1;
k = 25;
imax = 30;
for (iniapr = 1; iniapr <= 2; ++iniapr) {
h__ = h0;
de45d_c((U_fp)f_c, &m, &xn, yn, &xk, &k, &iniapr, &imax, &h__, y, rab);
dely[0] = yt[0] - y[0];
dely[1] = yt[1] - y[1];
relerr[0] = dely[0] / y[0];
relerr[1] = dely[1] / y[1];
/* Операторы вывода на печать: h0, k, iniapr, imax, y, yt, dely, relerr */
. . . . . . . . . . . . . . . . . . . . . . . . . .
}
return 0;
} /* main */
/*..........................................................................*/
int f_c(double *x, double *y, double *z__, int *m)
{
/* Parameter adjustments */
--z__;
--y;
/* Function Body */
z__[1] = *x / y[2];
z__[2] = -(*x) / y[1];
return 0;
} /* f_c */
Результаты:
--------------
h0= -1.00000000000000e-001 k= 25 iniapr= 1 imax= 30
y[0]= 3.00000000000001e+000 y[1]= 1.66666666666666e-001
yt[0]= 3.00000000000000e+000 yt[1]= 1.66666666666667e-001
dely=
-1.06581410364015e-014 3.05311331771918e-016
relerr=
-3.55271367880049e-015 1.83186799063151e-015
----------------------------------------------------------------------------
h0= -1.00000000000000e-001 k= 25 iniapr= 2 imax= 30
y[0]= 2.99999999999999e+000 y[1]= 1.66666666666667e-001
yt[0]= 3.00000000000000e+000 yt[1]= 1.66666666666667e-001
dely=
7.10542735760100e-015 -3.33066907387547e-016
relerr=
2.36847578586701e-015 -1.99840144432528e-015
----------------------------------------------------------------------------
Примеры 3 и 4 относятся к тому случаю, когда решение задачи Коши вычисляется при первом способе выбора начального приближения для коэффициентов Чебышёва производной (т.е. при iniapr = 1) с лучшей точностью, чем при iniapr = 2, т.е. с меньшей абсолютной и относительной погрешностью.