Подпрограмма: de51r_c
Назначение
Вычисление решения двухточечной краевой задачи для
нелинейного дифференциального уравнения второго порядка.
Математическое описание
Решается двухточечная краевая задача для нелинейного
дифференциального уравнения второго порядка
(1) y '' + f(x, y, y ') y ' + g(x, y, y ') y = r(x, y, y ') , xN≤ x ≤ xK
с заданными на концах отрезка интегрирования
[ xN, xK ] граничными условиями вида
(2) aN y '(xN) + bN y(xN) = CN ,
(3) aK y '(xK) + bK y(xK) = CK .
Решение вычисляется на заданной на отрезке
[ xN, xK ]
равномерной сетке с заданной абсолютной погрешностью EPS.
Используется метод линеаризации исходного уравнения (1) с
последующим применением метода конечных разностей для решения
линейного уравнения.
L.Fox, Proc. Roy. Soc. A 190, 31 - 59, 1947.
Использование
int de51r_c (S_fp f, integer *nx, real *xn, real *xk, real *cn,
real *ck, real *eps, integer *imax, integer *ider, real *y,
integer *ik, real *delty, real *rab, integer *ierr)
Параметры
|
f -
|
имя подпрограммы вычисления значений коэффициентов
f (x, y, y'), g (x, y, y'),
r (x, y, y') уравнения
(1). Первый оператор подпрограммы должен иметь вид:
|
|
|
int f (float *cf, float *cg, float *cr, float *x, float *y, int *i__,
float *rab)
Здесь:
|
|
x, y -
|
представляют первые два аргумента
функций f (x, y, y'),
g (x, y, y') и
r (x, y, y');
|
|
rab -
|
одномерный массив длины 11 * nx
(nx - число узлов сетки);
|
|
i -
|
указывает компоненту rab(i) массива
rab, которая представляет третий аpгумент функций
f (x, y, y'),
g (x, y, y'), r (x, y, y');
|
|
cf, cg, cr -
|
представляют вычисленные подпрограммой f значения
функций f (x, y, y'),
g (x, y, y') и r (x, y, y');
|
|
|
параметры cf, cg, cr, x, y и rab имеют тип
real, а i - integer;
|
|
nx -
|
число узлов равномерной сетки на отрезке
[xn, xk], (включая xn
и xk), в которых
требуется вычислить решение задачи;
nx ≥ 3 (тип: целый);
|
|
xn, xk -
|
концы отрезка интегрирования, в которых заданы
граничные условия (2) и (3), соответственно;
xn < xk (тип: вещественный);
|
|
cn, ck -
|
одномерные вещественные массивы длины 3, в
которых помещаются коэффициенты граничных условий
(2) и (3), соответственно, т.е.
cn (1) = an,
cn (2) = bn,
cn (3) = cn,
ck (1) = ak,
ck (2) = bk,
ck (3) = ck. Коэффициенты
должны удовлетворять условию
a2n + b2n ≠ 0
и a2k + b2k ≠ 0;
|
|
eps -
|
абсолютная погрешность, с которой требуется
вычислить решение краевой задачи; eps > 0 (тип:
вещественный);
|
|
imax -
|
указывает максимальное допустимое число
линеаризаций исходного уравнения (1) и число
итераций, используемых при решении получающегося
после линеаризации уравнения: число линеаризаций
pавно абсолютному значению | imax |
параметра imax, а число итераций pавно imax, если
imax > 0 и нулю, если imax < 0;
imax ≠ 0. Если
imax < 0, то
в результате работы подпрограммы значение imax
заменяется на противоположное; в этом случае imax
зажается переменной или элементом массива (тип: целый);
|
|
ider -
|
указывает на зависимость коэффициентов
f (x, y, y'), g (x, y, y'),
r (x, y, y') уравнения (1)
от производной y': если хотя бы один из этих
коэффициентов зависит от y', то ider > 0,
иначе ider = 0 (тип: целый);
|
|
y -
|
одномерный вещественный массив длины nx,
который должен содержать при обращении к
подпрограмме значения начального приближения решения
краевой задачи (1), (2), (3) во всех узлах
сетки. В результате работы подпрограммы в этом
массиве получаются значения конечного
приближения к решению (эти значения помещаются в y в
любом случае, даже и тогда, когда они не могут
быть вычислены с требуемой точностью);
|
|
ik -
|
целая переменная, содержащая число линеаризаций
уравнения (1), выполненных в действительности
подпрограммой для достижения требуемой точности eps
решения. В случае, когда данная точность
решения нелинейного уравнения не достигается за imax
линеаризаций или когда решение соответствующего
линейного уравнения не может быть вычислено с
точностью eps за imax итераций, то
значение параметра ik в результате работы
подпрограммы полагается равным - 1;
|
|
delty -
|
одномерный вещественный массив длины nx,
значение которого в результате работы подпрограммы
полагается равным разности между начальным и
конечным приближениями решения того линейного
уравнения, котоpое получается после выполнения
самой последней линеаризации исходного
уpавнения (1), причем эта разность вычисляется во
всех узлах сетки. В случае, когда imax < 0, все
элементы массива delty в результате работы
подпрограммы полагаются равными 0;
|
|
rab -
|
одномерный вещественный рабочий массив длины 11 * nx;
|
|
ierr -
|
целая переменная, служащая для сообщения об
ошибках, обнаруженных в процессе работы подпрограммы; при этом:
|
|
ierr=65 -
|
если в граничном условии (2) коэффициенты
an = bn = 0;
|
|
ierr=66 -
|
если в граничном условии (3) коэффициенты
ak = bk = 0;
|
|
ierr=67 -
|
если imax = 0;
|
|
ierr=68 -
|
если ider < 0;
|
|
ierr = 1 -
|
если заданная точность решения
нелинейного уравнения не достигается за
| imax |
линеаризаций или если решение соответствующего линейного
уpавнения не может быть вычислено с
точностью eps за imax итераций;
|
|
ierr = 3 -
|
если заданная точность решения
нелинейного уравнения не достигается за
| imax |
линеаризаций или если решение соответствующего линейного
уpавнения не может быть вычислено с
точностью eps за imax итераций, при этом
при аппроксимации исходного уравнения
линейным уравнением использовались
производные, вычисленные с меньшей точностью, чем это необходимо;
|
|
|
В каждом из выше указаных случаев решение задачи прекращается;
|
|
ierr = 2 -
|
если при аппроксимации исходного уpавнения линейным уравнением
использовались производные, вычисленные с
меньшей точностью, чем это
необходимо. При ierr = 1, 2, 3 полученное
приближенное решение запоминается в
массиве y; при желании интегрирование можно повторить обращением к
подпрограмме с новыми значениями параметров nx, y
или imax.
|
Версии
|
de51d_c -
|
вычисление решения двухточечной краевой задачи
для нелинейного дифференциального уравнения
второго порядка с повышенной точностью. При
этом параметры xn, xk, cn, ck, eps, y, delty,
rab должны иметь тип double.
|
Вызываемые подпрограммы
|
de50r_c -
|
вычисление решения двухточечной краевой задачи
для линейного дифференциального уравнения
второго порядка методом конечных разностей.
|
|
de50d_c -
|
вычисление решения двухточечной краевой задачи
для линейного дифференциального уравнения
второго порядка методом конечных разностей с повышенной точностью.
|
|
id10r_c -
|
подпрограмма вычисления первых и вторых
производных с помощью центральных разностей.
|
|
id10d_c -
|
подпрограмма вычисления первых и вторых
производных с помощью центральных разностей с двойной точностью.
|
|
|
Подпрограммы de50r_c и id10r_c вызываются подпрограммой
de51r_c, а подпрограммы de50d_c и id10d_c - подпрограммой de51d_c.
|
|
utde10_c -
|
подпрограмма выдачи диагностических сообщений
при работе подпрограммы de51r_c.
|
|
utde11_c -
|
подпрограмма выдачи диагностических сообщений
при работе подпрограммы de51d_c.
|
Замечания по использованию
|
|
Число узлов сетки nx, задаваемое пользователем при
обращении к подпрограмме, должно быть достаточно большим
для того, чтобы обеспечить сходимость
конечно - разностных аппроксимаций, используемых
для решения линейных уравнений, и последовательных
линеаризаций исходного уравнения (1).
В общем случае заданная точность eps не гарантируется.
Пусть ηi ,
i = 1, ..., nx, представляет разность между
текущим приближением к решению и приближением,
полученным в результате предыдущей линеаризации уравнения (1).
Тогда текущее приближение принимается за окончательное
решение, если для всех
i (1 ≤ i ≤ nx)
выполняется условие
| ηi | ≤ eps .
Значения параметров nx, xn, xk, cn, ck, eps, imax, в
результате работы подпрограммы сохраняются.
Kpоме того, de50r_c и de50d_c используют рабочие
подпрограммы id10rp_c и id10dp_c, соответственно.
|
Пример использования
Решается краевая задача:
y '' + ( xy' ) y' + e -y y = ( xy' + 2e-y - 2 ) cos x - ( 1 + 2xy' - e-y ) sin x
2y - y' = 0 при x = 0
y + y' = 1 при x = π / 2
(точное решение задачи y = sin x + 2cos x).
int main(void)
{
/* System generated locals */
int i__1;
/* Builtin functions */
double cos(double), sin(double);
/* Local variables */
extern int de51r_c(U_fp, int *, float *, float *, float *, float *,
float *, int *, int *, float *, int *, float *,
float *, int *);
static int ider, imax, ierr;
extern int f_c();
static float h__;
static int i__, n;
static float x, y[150], delty[150], ck[3], cn[3];
static int ik;
static float xk;
static int nx;
static float xn, yr[200], yx[150], rab[1650], eps;
nx = 150;
xn = 0.f;
xk = 3.14159265359f;
xk /= 2.f;
cn[0] = 2.f;
cn[1] = -1.f;
cn[2] = 0.f;
ck[0] = 1.f;
ck[1] = 1.f;
ck[2] = -1.f;
imax = 4;
eps = 1e-5f;
ider = 1;
i__1 = nx;
for (i__ = 1; i__ <= i__1; ++i__) {
yx[i__ - 1] = 1.5f;
/* L1: */
}
de51r_c((U_fp)f_c, &nx, &xn, &xk, cn, ck, &eps, &imax, &ider, yx, &ik,
delty, rab, &ierr);
printf("\n %5i \n", ik);
printf("\n %5i \n", ierr);
n = nx - 1;
h__ = (xk - xn) / n;
x = xn - h__;
i__1 = nx;
for (i__ = 1; i__ <= i__1; ++i__) {
x += h__;
y[i__ - 1] = (float)sin(x) + (float)cos(x) * 2.f;
yr[i__ - 1] = yx[i__ - 1] - y[i__ - 1];
printf("\n %16.7e %16.7e %16.7e %16.7e \n",
x, yx[i__-1], y[i__-1], yr[i__-1]);
/* L2: */
}
return 0;
} /* main */
int f_c(float *fc, float *gc, float *rc, float *x, float *y, int *i__,
float *rab)
{
/* Builtin functions */
double exp(double), cos(double), sin(double);
/* Parameter adjustments */
--rab;
/* Function Body */
*fc = *x * rab[*i__];
*gc = (float)exp(-(*y));
*rc = (*x * rab[*i__] - 2.f +
(float)exp(-(*y)) * 2.f) * (float)cos(*x) -
(*x * 2.f * rab[*i__] + 1.f -
(float)exp(-(*y))) * (float)sin(*x);
return 0;
} /* f_c */
Результаты: (приводятся только для первых 10 узлов)
yx(1) = 1.999999086
yx(2) = 2.010430009
yx(3) = 2.020637497
yx(4) = 2.030620414
yx(5) = 2.040377652
yx(6) = 2.049908127
yx(7) = 2.059210778
yx(8) = 2.068284572
yx(9) = 2.077128500
yx(10) = 2.085741581
ik = 3
ierr = 2