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

Назначение

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

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

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

          Y' = F (X, Y) ,

          Y = ( y1, ..., yM ) ,
          F = ( f1 (X, y1,..., yM),..., fM (X, y1,..., yM) )     методом Гира.

Метод Гира для нежесткой системы является многошаговым предсказывающе - исправляющим методом Адамса, записанным в форме Нордсика, при этом предсказание и исправление имеют один и тот же порядок.

B случае, когда система уравнений является жесткой, интегрирование осуществляется специальным методом основанным на методе типа Адамса и использующим якобиан ( ∂F/∂Y ) системы, который вычисляется подпрограммой по формулам численного дифференцирования.

По заданному значению YX решения в узле X_n вычисляется значение этого решения в узле X_n + H. Считается, что решение в узле X_n + H вычислено с требуемой точностью EPS, если выполняется следующее соотношение:

                M   
             (  ∑ ( eI / YPM(I) )2 )1/2  ≤  EPS  ,
                I=1 

где e_I - вычисляемая подпрограммой оценка абсолютной погрешности приближенного значения I - й компоненты решения в узле X_n + H, YPM (I) - задается при обращении к подпрограмме и, по существу, позволяет использовать разный тип погрешности для каждой компоненты. Если, например, YPM(I) = 1, то для I - ой компоненты YX(I) решения будет использоваться абсолютная погрешность; если YPM (I) = | YX (I) | ≠ 0 , где YX(I) - значение I - й компоненты YX(I) в узле X_n, то для этой компоненты берется отношение абсолютной погрешности приближенного ее значения в узле X_n + H к абсолютной величине приближенного значения этой компоненты в предыдущем узле Х_n. Значение H может быть меньше или pавно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.

Gear C.W., The automatic integration of ordinary differential equations. Communications of the ACM, 14, 3 (March 1971), 176-179.

Gear C.W., Numerical initial value problems in ordinary differential equations, Prentice - Hall, Englewood Cliffs, N.J., 1971.

Gear C.W., The automatic integration of stiff ordinary differential equations. Information Processing 68, A.J.H.

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

    int de21r_c (S_fp f, integer *m, integer *istifj,
                 integer *iorder, integer *jstart, real *hmin, real *hmax, real *eps,
                 real *yx, real *x, real *h, logical *bul, real *ypm, real *delty,
                 real *rab, real *yp, integer *ierr)

   Параметры 

   Версии 

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

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

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

         y1'  =   y3
         y2'  =   y4
         y3'  =  - ( 1 - 0.002 x ) y1 / ( y12  +  y22 )3/2
         y4'  =  - ( 1 - 0.002 x ) y2 / ( y12  +  y22 )3/2  ,

        y1(0)  =  0.09411764706 ,     y2(0)  =  0 ,     y3(0)  =  0 ,     y4(0)  =  4.5 

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

int main(void)
{
    /* System generated locals */
    int i__1;

    /* Local variables */
    extern int de21r_c(U_fp, int *, int *, int *, int *, float *, float *,
                       float *, float *, float *, float *, logical *,
                       float *, float *, float *, float *, int *);
    static float hmin, hmax;
    static int ierr;
    extern int f_c();
    static float h__;
    static int i__, m;
    static float x, y[4], delty[4], yp[32];
    static int iorder, istifj, jstart;
    static float rab[68];
    static logical bul;
    static float eps, ypm[4];

    bul = FALSE_;
    m = 4;
    x = 0.f;
    h__ = .01f;
    y[0] = .09411764706f;
    y[1] = 0.f;
    y[2] = 0.f;
    y[3] = 4.5f;
    iorder = 4;
    jstart = 0;
    hmax = .1f;
    hmin = 1e-18f;
    eps = 1e-5f;
    istifj = 0;
    i__1 = m;
    for (i__ = 1; i__ <= i__1; ++i__) {
/* L10: */
        ypm[i__ - 1] = 1.f;
    }
    de21r_c((U_fp)f_c, &m, &istifj, &iorder, &jstart, &hmin, &hmax, &eps, y,
            &x, &h__, &bul, ypm, delty, rab, yp, &ierr);

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

int f_c(float *x, float *y, float *dy, int *m)
{
    /* Builtin functions */
    double sqrt(double);

    /* Local variables */
    static float r__, r2;

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

    /* Function Body */
    r2 = y[1] * y[1] + y[2] * y[2];
    r2 *= (float)sqrt(r2);
    r__ = -(1.f - *x * .002f) / r2;
    dy[1] = y[3];
    dy[2] = y[4];
    dy[3] = r__ * y[1];
    dy[4] = r__ * y[2];
    return 0;
} /* f_c */


Результаты:

          x  =   4.057708603-05
          h  =   4.057708603-05

          y(1)  =   9.411746119-02
          y(2)  =   1.825965265-04
          y(3)  =  -4.580764462-03
          y(4)  =   4.499991113

          ierr  =   0