|
Текст подпрограммы и версий de50r_c.zip , de50d_c.zip |
Тексты тестовых примеров tde50r_c.zip , tde50d_c.zip |
Вычисление решения двухточечной краевой задачи для линейного дифференциального уравнения второго порядка методом конечных разностей.
Решается двухточечная краевая задача для линейного дифференциального уравнения второго порядка
(1) y '' + f(x) y ' + g(x) y = r(x) , xN≤ x ≤ xK
c заданными на концах отрезка интегрирования [ xN, xК ] граничными условиями вида
(2) aN y '(xN) + bN y(xN) = cN , a2N + b2N ≠ 0 (3) aK y '(xK) + bK y(xK) = cK , a2K + b2K ≠ 0
- методом конечных разностей.
Решение вычисляется на заданной на отрезке [ xN, xK ] равномерной сетке с заданной абсолютной погрешностью EPS.
L.Fox, Proc. Roy. Soc. A 190, 31 - 59, 1947.
int de50r_c (integer *nx, real *xn, real *xk, real *cn,
real *ck, real *f, real *g, real *r, real *eps, integer *imax,
real *y, integer *ik, real *delty, real *rab, integer *ierr)
Параметры
| nx - | число узлов равномерной сетки на отрезке [ xn, xk ], (включая xn и xk), в которых требуется вычислить решение задачи. Предполагается, что узлы расположены в порядке возрастания: xn = x1 < x2< ... < xnx = 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; |
| f, g, r - |
одномерные вещественные масивы длины nx, в
которых перед обращением к подпрограмме должны
быть помещены значения коэффициентов уравнения
(1), вычисленные в узлах равномерной сетки, на
которой требуется получить решение краевой задачи, а именно:
f (1) = f ( x1), f (2) = f ( x2 ) ,..., f (nx) = f ( xnx );
g (1) = g ( x1), g (2) = g ( x2 ) ,..., g (nx) = g ( xnx );
r (1) = r ( x1), r (2) = r ( x2 ) ,..., r (nx) = r ( xnx );
|
| eps - | абсолютная погрешность, с которой требуется вычислить решение краевой задачи (тип: вещественный); |
| imax - | максимальное допустимое число разностных итераций, котоpое разрешается использовать в процессе работы подпрограммы. Если эти итерации не нужны, то параметр imax при обращении к подпрограмме задается равным 0; в этом случае значение параметра eps при работе подпрограммы не используется (тип: целый); |
| y - | одномерный вещественный массив длины nx, в котоpом в результате работы подпрограммы получаются приближеннные значения решения краевой задачи, вычисленные на сетке (эти значения помещаются в y в любом случае, даже и тогда, когда они не могут быть вычислены с требуемой точностью); |
| ik - | число разностных итераций, использованных в действительности подпрограммой для достижения требуемой точности eps решения. B случае, когда данная точность не может быть достигнута за imax итераций, значение параметра ik в результате работы подпрограммы полагается равным - 1. Если imax = 0, то и значение ik в результате выполнения программы тоже полагается равным 0 (тип: целый); |
| delty - | одномерный вещественный массив длины nx, значение которого в результате работы подпрограммы полагается равным разности между начальным приближением решения и значением его последней итерации, вычисленной во всех nx узлах сетки. В случае imax = 0 все элементы массива delty полагаются равными нулю; |
| rab - | одномерный вещественный рабочий массив длины 8 * nx; |
| ierr - | целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом: |
| ierr=65 - | если в граничном условии (2) коэффициенты an = bn = 0; |
| ierr=66 - | если в граничном условии (3) коэффициенты ak = bk = 0; |
| ierr = 1 - | если приближенное решение задачи, вычисленное за imax итераций, не достигает требуемой точности eps; |
| В каждом из этих случаев решение задачи прекращается. При ierr = 1 интегрирование можно повторить обращением к подпрограмме с новыми значениями параметров nx, f, g, r или imax. |
Версии
| de50d_c - | вычисление решения двухточечной краевой задачи для линейного дифференциального уравнения второго порядка методом конечных разностей с повышенной точностью. При этом параметры xn, xk, cn, ck, f, g, r, eps, y, delty, rab должны иметь тип double. |
Вызываемые подпрограммы
| as09r_c - | вычисление решения систем линейных алгебраических уравнений с ленточной матрицей. |
| as09d_c - | вычисление решения систем линейных алгебраических уравнений с ленточной матрицей с двойной точностью. |
| i i10r_c - | подпрограмма вычисления центральных разностей таблично заданной функции. |
| i i10d_c - | подпрограмма вычисления центральных разностей таблично заданной функции с двойной точностью. |
| Подпрограммы as09r_c и i i10r_c вызываются подпрограммой de50r_c, а подпрограммы as09d_c и i i10d_c - подпрограммой de50d_c. |
| utde10_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de50r_c. |
| utde11_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de50d_c. |
Замечания по использованию
|
При задании значений коэффициентов уравнения в узлах равномерной сетки последняя должна быть достаточно густой для того, чтобы обеспечить сходимость метода конечных разностей. В общем случае заданная точность eps приближенного pешения не гарантируется. Если для двух последовательных итераций yi ( j ) и yi ( j - 1 ) разность ei ( j ) = yi ( j ) - yi ( j - 1 ) удовлетворяет критерию | ei( j ) | ≤ eps для i = 1, ..., nx , то yi( j ) принимается в качестве приближенного решения краевой задачи. Если же после выполнения imax итераций точность eps не достигается, то, тем не менее, значение приближенного решения, полученное на последней итерации, помещается в массив y. B этом случае параметр ierr = 1, а ik = - 1. Значения параметров nx, xn, xk, cn, ck, f, g, r, eps, imax сохраняются. |
Решается краевая задача:
y '' + x2 y ' + x y = 10e- x*x ( 2x2 - 1 ) ( 2 - x ) , - 1≤ x ≤ 1
y ' - 2y = 0 при x = - 1
y ' + 2y = 0 при x = 1
(точное решение задачи y (x) = 10e- x*x )
int main(void)
{
/* Builtin functions */
double exp(double);
/* Local variables */
extern int de50r_c(int *, float *, float *, float *, float *, float *,
float *, float *, float *, int *, float *, int *,
float *, float *, int *);
static int imax, ierr;
static float f[200], g[200], h__;
static int i__, n;
static float r__[200], x, y[200], delty[200], ck[3], cn[3];
static int ik;
static float xk, xn;
static int nx;
static float yr[200], yx[200], rab[1600], eps;
int i__1;
float r__1, r__2;
xn = -1.f;
xk = 1.f;
nx = 200;
cn[0] = 1.f;
cn[1] = -2.f;
cn[2] = 0.f;
ck[0] = 1.f;
ck[1] = 2.f;
ck[2] = 0.f;
imax = 4;
eps = 1e-5f;
n = nx - 1;
h__ = (xk - xn) / n;
x = xn - h__;
i__1 = nx;
for (i__ = 1; i__ <= i__1; ++i__) {
x += h__;
/* Computing 2nd power */
r__1 = x;
f[i__ - 1] = r__1 * r__1;
g[i__ - 1] = x;
/* Computing 2nd power */
r__1 = x;
/* Computing 2nd power */
r__2 = x;
r__[i__ - 1] = (float)exp(-(r__1 * r__1)) * 10.f *
(r__2 * r__2 * 2.f - 1.f) * (2.f - x);
/* L1: */
}
de50r_c(&nx, &xn, &xk, cn, ck, f, g, r__, &eps, &imax, yx, &ik, delty,
rab, &ierr);
printf("\n %5i \n", ik);
x = xn - h__;
i__1 = nx;
for (i__ = 1; i__ <= i__1; ++i__) {
x += h__;
/* Computing 2nd power */
r__1 = x;
y[i__ - 1] = (float)exp(-(r__1 * r__1)) * 10.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 */
Результаты: (приводятся только для первых 10 узлов)
yx(1) = 3.678799304
yx(2) = 3.753113958
yx(3) = 3.828156408
yx(4) = 3.903910587
yx(5) = 3.980359671
yx(6) = 4.057486083
yx(7) = 4.135271487
yx(8) = 4.213696790
yx(9) = 4.292742144
yx(10) = 4.372386948
ik = 2
ierr = 0