Подпрограмма:  de01r_c (версия de85r_c)

Назначение

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

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

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

(1)                             Y '  =  F(X, Y) ,

          Y  =  ( y1,..., yM ) ,   F  =  ( f1( X, y1,..., yM ),..., fM( X, y1,..., yM ) ) 

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

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

Предполагается, что среди характеристических корней матрицы Якоби ∂F / ∂y функции F имеются большие по модулю корни. Решение вычисляется в одной точке XK, которая является концом интервала интегрирования. Предполагается также, что система (1) может быть достаточно хорошо аппроксимирована линейной однородной системой с кусочнопостоянными коэффициентами, т.е. на каждом из [ x_n, x_n + H ], на которые разбивается весь интервал интегрирования, она может быть аппроксимирована системой вида

                                     Y '  =   AY  , 

где A - постоянная матрица, вообще говоря, своя для каждого отрезка. Для интегрирования системы применяется метод Лоусона.

Метод Лоусона является одношаговым методом и заключается в следующем. Допустим, что искомое решение системы (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)  =  F1(x, z)  =  exp [ - ( x - xn ) A0 ] { F(x, exp [ ( x - xn ) A0 ] Z(x)) -
                                      - A0 exp [ ( x - xn ) A0 ] Z(x) } 

Матрица Якоби ∂ F1 / ∂Z системы (2) и матрица Якоби ∂F / ∂y системы (1) связаны между собой соотношением

         ∂ F1 / ∂Z  =  exp [ - (x - xn) A0] { ∂F / ∂y - A0 } exp [ ( x - xn ) A0 ] 

Указанное преобразование выполняется самой подпрограммой. В качестве матрицы A⁰ подпрограмма выбирает матрицу ∂F / ∂y, вычисленную либо в точке ( x_n, y_n ), либо в некоторой предыдущей точке ( x_n - k, y_n - k ). Если шаг H достаточно мал, то данное преобразование позволяет уменьшить характеристические корни матрицы ∂ F1 / ∂Z по сравнению с характеристическими корнями матрицы ∂F / ∂y. А это приводит к уменьшению константы Липшица системы (2) по сравнению с константой Липшица системы (1).

Для решения системы (2) применяются классический метод Рунге - Кутта четвертого порядка точности, причем одновременно с решением (2) производится обратное преобразование от функции Z (x) к функции Y (x).

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

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

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

    int de01r_c (real *f, real *fj, integer *m, real *xn, real *yn,
                 real *xk, real *hmin, real *eps, real *p, integer *n, real *h, 
                 real *y, real *r, integer *ierr)

   Параметры 

   Версии 

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

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

                     Y '  =  AY 

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

                                     
     y1'  =  (- 500 + 2 x) y1  +  (1 - x /10 ) y2  +  (1 + 1/(x +    ) sin y1 ,
     y1(   = 49 
     y2'  =  (- 1000 - 10 x) y1  +  (1 + x /10 ) y2  +  (2  + 1/(x + 2) ) cos y2 ,
     y2(   = 50 

              1 ≤ x ≤ 6  

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

int main(void)
{
    /* Local variables */
    extern int de01r_c(U_fp, U_fp, int *, float *, float *, float *,
                       float *, float *, float *, int *, float *, float *,
                       float *, int *);
    static float hmin;
    static int ierr;
    extern int f_c();
    static float h__;
    static int m, n;
    static float p, r__[35], y[2];
    extern int fj_c();
    static float xk, xn, yn[2], eps;

    m = 2;
    xn = 1.f;
    yn[0] = 49.f;
    yn[1] = 50.f;
    hmin = 1e-10f;
    eps = 1e-4f;
    p = 100.f;
    h__ = .01f;
    xk = 6.f;
    n = 2;
    de01r_c((U_fp)f_c, (U_fp)fj_c, &m, &xn, yn, &xk, &hmin, &eps, &p, &n,
            &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 f_c(float *x, float *y, float *dy, int *m)
{
    /* Builtin functions */
    double sin(double), cos(double);

    /* Parameter adjustments */
    --dy;
    --y;

    /* Function Body */
    dy[1] = (*x * 2.f - 500.f) * y[1] + (1.f - *x * .1f) * y[2] +
            (1.f / (*x + 1.f) + 1.f) * (float)sin(y[1]);
    dy[2] = (-1e3f - *x * 10.f) * y[1] + (*x * .1f + 1.f) * y[2] +
            (1.f / (*x + 2.f) + 2.f) * (float)cos(y[2]);
    return 0;
} /* f_c */

int fj_c(float *x, float *y, float *df, int *m)
{
    /* Builtin functions */
    double cos(double), sin(double);

#define df_ref(a_1,a_2) df[(a_2)*2 + a_1]

    /* Parameter adjustments */
    df -= 3;
    --y;

    /* Function Body */
    df_ref(1, 1) = *x * 2.f - 500.f +
                   (1.f / (*x + 1.f) + 1.f) * (float)cos(y[1]);
    df_ref(1, 2) = 1.f - *x * .1f;
    df_ref(2, 1) = -1e3f - *x * 10.f;
    df_ref(2, 2) = *x * .1f + 1.f -
                   (1.f / (*x + 2.f) + 2.f) * (float)sin(y[2]);
    return 0;
} /* fj_c */


Результаты: 

             y(1)                                    y(2)                              h 
     -4.175740488392-02       -5.087148371158+01      2.500000000001-03 
      ierr = 0 

Результаты решения этой же задачи, полученные с помощью подпрограммы de85r_c

     -4.175740472033-02       -5.087148363417+01      5.000000000003-03 
      ierr = 0

Результаты решения этой же задачи, полученные с помощью подпрограммы de25r_c

     -4.175740429849-02       -5.087148203468+01      1.756831014061-03 
      ierr = 0