Подпрограмма:  de31r_c

Назначение

Вычисление решения задачи Коши для жестких линейных однородных систем обыкновенных дифференциальных уравнений первого порядка с переменными коэффициентами экспоненциальным методом.

Математическое описание

Решается задача Коши для жесткой линейной однородной системы М обыкновенных дифференциальных уравнений первого порядка с переменными коэффициентами

(1)          Y ' = A(x) Y 

здесь Y = ( y₁, ..., y_M ), А(x) - квадратная вещественная матрица порядка M, элементами которой являются трижды непрерывно дифференцируемые функции. Начальные значения заданы в точке XN :

          Y(XN) = YN ,     YN = ( y10,..., yM0 ) , 

Предполагается, что матрица A (x) имеет все собственные числа с отрицательными вещественными частями, среди которых встречаются большие по абсолютной величине. Решение вычисляется в одной точке XK (которая является концом интервала интегрирования).

Интервал интегрирования разбивается на NS равных частей длины T, на каждой из которых исходная система уравнений (1) аппроксимируется системой уравнений с постоянными коэффициентами

               ~         ~   ~         ~
(2)          Y '  =  A  Y ,      A = CONST. 

Решение системы (2) вычисляется с помощью матричной экспоненты

                         ~
(3)          EXP( A T )  =
                                               ~           N
                             =  [ EXP(  A  H ) ]      ,   H = T/N ,   N > 0 

Значения N и NS задаются пользователем при обращении к подпрограмме, причем значение N выбирается таким образом, чтобы выполнялось условие

(4)               max        || A(x) H ||  ≤  1
           XN ≤ x ≤ XK 

(в качестве || A(x) H || можно взять максимальную сумму модулей элементов матрицы A (x) H по стpокам).

С.Ф.Залеткин, O численном решении задачи Коши для обыкновенных линейных однородных дифференциальных уравнений на больших отрезках интегрирования, "Вычислительные методы и программирование", вып.26, Изд-во МГУ, 1977.

Использование

    int de31r_c (S_fp fa, integer *m, real *xn, real *yn, real *xk,
                 integer *n, integer *ns, real *y, real *a, real *e, real *r,
                 real *r1, real *r2, real *r3, integer *ierr)

   Параметры 

   Версии 

   Вызываемые подпрограммы 

   Замечания по использованию  

Пример использования

Использование подпрограммы иллюстрируется на примере:

           y1'  =  - (2 + x) y1 / (1 + x) + 20 x y2 ,
           y2'  =  -20 x y1 - (2 + x) y2 / (1 + x) ,

          y1(0)  =  0 ,     y2(0)  =  1 ,     0 ≤ x ≤ 6  . 

Приводятся подпрограмма вычисления матрицы системы и фрагмент вызывающей программы, а также результаты счета.

int main(void)
{
    /* Local variables */
    extern int de31r_c(U_fp, int *, float *, float *, float *, int *, int *,
                       float *, float *, float *, float *, float *, float *,
                       float *, int *);
    static int ierr;
    static float a[4] /* was [2][2] */, e[4] /* was [2][2] */;
    static int m, n;
    static float r__[2], y[2], r1[4] /* was [2][2] */,
                               r2[4] /* was [2][2] */, r3[2];
    extern int fa_c();
    static float xk, xn, yn[2];
    static int is1;

    m = 2;
    xn = 0.f;
    yn[0] = 0.f;
    yn[1] = 1.f;
    xk = 6.f;
    n = 128;
    is1 = 128;
    de31r_c((U_fp)fa_c, &m, &xn, yn, &xk, &n, &is1, y, a, e, r__, r1, r2, r3,
            &ierr);

    printf("\n %15.7e \n", xk);
    printf("\n %15.7e %15.7e \n", y[0], y[1]);
    return 0;
} /* main */

int fa_c(float *a, float *x, int *m)
{
    /* System generated locals */
    int a_dim1, a_offset;

#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]

    /* Parameter adjustments */
    a_dim1 = *m;
    a_offset = 1 + a_dim1 * 1;
    a -= a_offset;

    /* Function Body */
    a_ref(1, 1) = -(*x + 2.f) / (*x + 1.f);
    a_ref(*m, *m) = a_ref(1, 1);
    a_ref(*m, 1) = *x * -20.f;
    a_ref(1, *m) = -a_ref(*m, 1);
    return 0;
} /* fa_c */


Результаты:

          y(1)  =   0.3395896401 * 10-3
          y(2)  =  -0.1004661347 * 10-3

          ierr  =  0