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

Назначение

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

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

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

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

многошаговым предсказывающе - исправляющим методом пятого порядка точности. Суть метода состоит в том, что сначала по явной формуле Адамса - Штермера строятся предсказанные значения приближенного решения, которые затем исправляются по неявной формуле Адамса - Штермера.

По значению решения  YX в узле  X_n вычисляется значение этого решения в узле  X_n + H. При этом, если  X_n является началом  XN интервала интегрирования, то для вычисления решения в   XN + H требуется еще и значение производной решения в  XN. Вместе с решением подпрограмма определяет также всю информацию, необходимую для вычисления приближенного решения в любом следующем узле.

Bсе компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой заданной константы  P (называемой границей перехода), то контроль точности для этой компоненты ведется по относительной погрешности, иначе  -  по абсолютной. Значение  H может быть меньше или pавно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.

Березин И.С., Жидков Н.П. Методы вычислений, т.2, Физматгиз, М., 1960.

Бахвалов H.C., Численные методы, т.1, "Hаука", 1975.

Жоголев E.A. Программа интегрирования систем обыкновенных дифференциальных уравнений второго порядка методом Штермера. Сб. "Вычислительные методы и программирование", вып. 1, Изд - во МГУ, 1962.

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

    int de42r_c (S_fp f, integer *m, integer *jstart, real *hmin,
                 real *hmax, real *eps, real *p, real *yx, real *x, real *dy,
                 real *h, logical *bul, real *delty, real *df, real *rfn, real *rf,
                 real *yp, real *dyp, integer *ierr)

   Параметры 

   Версии 

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

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

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

     y1'' = - y1   ,
     y2'' = - y2   ;

     y1 (3/4 π) = √2 / 2   ,      y2 (3/4 π) = - √2 / 2   ,
     y1' (3/4 π) = - √2 / 2   ,   y2' (3/4 π) = - √2 / 2   .

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

int main(void)
{
    /* Builtin functions */
    double sqrt(double), sin(double), cos(double);

    /* Local variables */
    extern int de42r_c(U_fp, int *, int *, float *, float *, float *,
                       float *, float *, float *, float *, float *,
                       logical *, float *, float *, float *, float *,
                       float *, float *, int *);
    static float hmin, hmax;
    static int ierr;
    extern int f_c();
    static float h__;
    static int m;
    static float p, x, y[2], delty[2], y1, y2, df[10];
    static int ih;
    static float rf[2], dy[2], yp[2];
    static int jstart;
    static logical bul;
    static float rfn[2], eps, dyp[2];

    h__ = -.01f;
    hmin = 1e-5f;
    m = 2;
    x = 2.3561947499999998f;
    y[0] = (float)sqrt(2.f) / 2.f;
    y[1] = -y[0];
    eps = 1e-4f;
    dy[0] = y[1];
    dy[1] = -y[0];
    ih = 0;
    bul = FALSE_;
    jstart = 0;
    p = 100.f;
    hmax = .1f;
l6:
    de42r_c((U_fp)f_c, &m, &jstart, &hmin, &hmax, &eps, &p, y, &x, dy, &h__,
            &bul, delty, df, rfn, rf, yp, dyp, &ierr);
    ++ih;
    y1 = (float)sin(x);
    y2 = (float)cos(x);

    printf("\n %16.7e %16.7e %16.7e \n", x, y[0], y[1]);
    printf("\n %16.7e %16.7e %16.7e \n\n", h__, y1, y2);
    if (ih == 1) {
        goto l6;
    }
    if (ih == 5) {
        goto l20;
    }
    jstart = -1;
    if (ih == 4) {
        goto l8;
    }
    if (ih == 3) {
        goto l7;
    }
    h__ = .005f;
    h__ = -h__;
    goto l6;
l7:
    h__ = .02f;
    goto l6;
l8:
    h__ = .01f;
    goto l6;
/* l53: */
    if (ierr != 0) {
        goto l20;
    }
l20:
    return 0;
} /* main */

int f_c(float *x, float *y, float *dy, int *m)
{
    /* Parameter adjustments */
    --dy;
    --y;

    /* Function Body */
    dy[1] = -y[1];
    dy[2] = -y[2];
    return 0;
} /* f_c */


Результаты:

после первого обращения к подпрограмме -

                   x                            y(1)                              y(2)
       2.346194490 + 00      7.141423761 - 01      - 7.000004762 - 01

                   h__                          y1                                 y2
     - 1.000000000 - 02      7.141423761 - 01      - 7.000004762 - 01

после второго обращения к подпрограмме -

                   x                            y(1)                              y(2)
       2.336194490 + 00      7.211065574 - 01      - 6.928241717 - 01

                   h__                          y1                                 y2
     - 1.000000000 - 02      7.211065574 - 01      - 6.928241717 - 01

после третьего обращения к подпрограмме -

                   x                           y(1)                               y(2)
       2.341194490 + 00      7.176334371 - 01      - 6.964210292 - 01

                   h__                         y1                                  y2
     - 5.000000000 - 03      7.176334371 - 01      - 6.964210292 - 01

после четвертого обращения к подпрограмме -

                   x                           y(1)                               y(2)
       2.366194490 + 00      7.000004762 - 01      - 7.141423761 - 01

                   h__                         y1                                  y2
       2.000000000 - 02      7.000004762 - 01      - 7.141423761 - 01

после пятого обращения к подпрограмме -

                   x                           y(1)                                y(2)
       2.356194490 + 00      7.071067812 - 01      - 7.071067812 - 01

                   h__                         y1                                   y2
       1.000000000 - 02      7.071067812 - 01      - 7.071067812 - 01