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

Назначение

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

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

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

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

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

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

     Y ' (XN) = DYN   ,   DYN = ( y10', ..., yM0' )   ,   -

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

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

Решение и производная вычисляются в одной точке  XK, которая является концом интервала интегрирования.

Bсе компоненты решения проверяются на точность по меpе погрешности, т.е. если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой заданной константы  P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.

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

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

    int de47r_c (real *f, integer *m, real *xn, real *yn,
                 real *dyn, real *xk, real *hmin, real *hmax, real *eps, real *p,
                 real *h, real *y, real *dy, real *z, real *delty, real *df,
                 real *rfn, real *rf, real *yp, real *dyp, real *zp, integer *ierr)

   Параметры 

   Версии 

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

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

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

     y1'' = 0.5 ( y2 - y1 + y1' + y2' ) - 0.5   ,   y1 (0) = 1   ,   y1' (0) = 1.5   ,
     y2 = y1'' - 0.5 x   ,   y2 (0) = - 1   ,   y2 (0) = - 0.5   .

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

    /* Local variables */
    extern int de47r_c(U_fp, int *, float *, float *, float *, float *,
                       float *, float *, float *, float *, float *, float *,
                       float *, float *, float *, float *, float *, float *,
                       float *, float *, float *, int *);
    static float hmin, hmax;
    static int ierr;
    static float h__;
    static int m;
    static float p, y[2], z__[2], delty[2];
    extern int f2_c();
    static float y1, y2, df[10] /* was [2][5] */, rf[2], dy[2], xk, xn, yn[2],
                  yp[2], zp[2], dy1, dy2, rfn[2], eps, dyn[2], dyp[2];

    m = 2;
    xn = 0.f;
    yn[0] = 1.f;
    yn[1] = -1.f;
    dyn[0] = 1.5f;
    dyn[1] = -.5f;
    xk = 5.f;
    hmin = 1e-10f;
    hmax = .1f;
    eps = 1e-4f;
    p = 100.f;
    h__ = .01f;
    de47r_c((U_fp)f2_c, &m, &xn, yn, dyn, &xk, &hmin, &hmax, &eps, &p, &h__,
            y, dy, z__, delty, df, rfn, rf, yp, dyp, zp, &ierr);

    printf("\n %5i \n", ierr);
    printf("\n %16.7e \n", xk);
    printf("\n %16.7e %16.7e \n", y[0], y[1]);
    printf("\n %16.7e %16.7e \n", dy[0], dy[1]);
    y1 = (float)sin(xk) + (float)cos(xk) + xk * .5f;
    y2 = -y1 + xk;
    dy1 = (float)cos(xk) - (float)sin(xk) + .5f;
    dy2 = -dy1 + 1.f;

    printf("\n\n %16.7e \n", xk);
    printf("\n %16.7e %16.7e \n", y1, y2);
    printf("\n %16.7e %16.7e \n", dy1, dy2);
    return 0;
} /* main */

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

    /* Function Body */
    z__[1] = (y[2] + dy[1] - y[1] + dy[2]) * .5f - .5f;
    z__[2] = y[1] - *x * .5f;
    return 0;
} /* f2_c */


Результаты:

      ierr = 0

                      xk                              y(1)                             y(2)
       5.000000000000 + 00     1.824737966524 + 00    3.175262014483 + 00

                                                          y1                                y2
                                              1.824737910809 + 00      3.175262089189 + 00

                                                         dy(1)                            dy(2)
                                              1.742586640374 + 00    - 7.425866599151 - 01

                                                         dy1                              dy2
                                              1.742586460128 + 00    - 7.425864601282 - 01