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

Назначение

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

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

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

(1)                         Y ' (X)   =   A(X) * Y(X) + φ(X) ,

                                 Y   =   ( y1,..., yM ) , 
                                 A(X)   =   ( ai j(X) ) ,     i, j  =  1, ..., M ,
                                φ(X)    =   ( φ1(X), ... , φM(X) ) 

с начальными условиями, заданными в точке XN:

                         Y(XN)   =   YN ,   YN   =   ( y1 0,..., yM 0 ) . 

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

Метод Лоусона является одношаговым методом и заключается в следующем. Допустим, что искомое решение системы (1) уже вычислено в некоторой точке  x = x_n  интервала интегрирования, т.е. известно  y_n ≈ y (x_n). Для отыскания решения Y(x_n + 1) = Y(x_n + H) в следующем узле  x_n + 1 = x_n + H выполняются такие действия. Исходная система уравнений с помощью замены искомой функции Y (x) на  x_n ≤ x ≤  x_n + H  по формуле

                             Y(x)  =  exp [ ( x - xn ) A0 ] Z(x)  ,

где A⁰ - некоторая постоянная матрица, преобразуется в систему уравнений относительно новой неизвестной функции Z (X):

(2)  Z ' (x)  =  A1(x) Z(x)  +  φ1(x)  =  exp [ - ( x - xn ) A0 ] { A(x) - A0 }
                                   exp [ ( x - xn ) A0 ] Z(x) + exp [ ( - ( x - xn ) A0 ] φ(x)

           xn ≤ x ≤ xn + H  

Данное преобразование выполняется самой подпрограммой. В качестве матрицы A⁰ подпрограмма выбирает матрицу A⁰ = A (x_n + H /2). Если шаг H достаточно мал, то преобразование позволяет уменьшить характеристические корни матрицы A1 (x) по сравнению с характеристическими корнями исходной матрицы A (x). Это приводит к уменьшению константы Липшица системы (2) по сравнению с константой Липшица системы (1). Для решения системы (2) применяются формулы классического метода Рунге - Кутта четвертого порядка точности, причем одновременно с решением (2) производится обратное преобразование от функции Z (x) к функции Y (x).

Все компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Абсолютная погрешность приближенного решения оценивается по правилу Рунге.

J.Douglas Louson. Generalized Runge - Kutta processes for stable systems with large Lipshitz constants, SIAM Journal on Numerical Analisys. Vol 4, No.3, 1967.

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

    int de05r_c(real *fa, real *fi, integer *m, real *xn,
                real *yn, real *xk, real *hmin, real *eps, real *p, real *h,
                real *y, real *r, integer *ierr)

   Параметры 

   Версии 

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

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

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

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

       0 ≤ x ≤ 6 

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

int main(void)
{
    /* Local variables */
    extern int de05r_c(U_fp, U_fp, int *, float *, float *, float *,
                       float *, float *, float *, float *, float *,
                       float *, int *);
    static float hmin;
    static int ierr;
    static float h__;
    static int m;
    static float p, r__[37], y[2];
    extern int fa_c(), fi_c();
    static float xk, xn, yn[2], eps;

    m = 2;
    xn = 0.f;
    yn[0] = 2.f;
    yn[1] = 18.f;
    hmin = 1e-10f;
    eps = 1e-5f;
    p = 100.f;
    xk = 6.f;
    h__ = .01f;
    de05r_c((U_fp)fa_c, (U_fp)fi_c, &m, &xn, yn, &xk, &hmin, &eps, &p, &h__,
            y, r__, &ierr);

    printf("\n %16.7e %16.7e \n", y[0], y[1]);
    printf("\n %16.7e \n", h__);
    printf("\n %5i \n", ierr);
    return 0;
} /* main */

int fa_c(float *a, float *x, int *m)
{
#define a_ref(a_1,a_2) a[(a_2)*2 + a_1]

    /* Parameter adjustments */
    a -= 3;

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

#undef a_ref

int fi_c(float *r1, float *x, int *m)
{
    /* Parameter adjustments */
    --r1;

    /* Function Body */
    r1[1] = 0.f;
    r1[2] = 0.f;
    return 0;
} /* fi_c */

 
   Результаты: 

             y(1)                                y(2)                             h 
      5.911150077338-03    -2.487346326713-03     1.600000000001-01 

      ierr = 0