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

Назначение

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

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

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

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

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

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

отыскиваются необходимые при численном интегрировании этой системы методом Штермера значения решения  Y и производной  Y ' в первом (после начальной точки  XN) узле интегрирования   ХN + Н, первая разностная производная назад и конечные разности правой части системы до четвертого порядка включительно в том же узле  XN + H.

Четвертый порядок разностей соответствует пятому порядку точности используемого метода Штермера.

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

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

Первая формула представляет значение решения в узле  ХN + Н через начальные условия (т.е. через его значение и значение его производной в начальной точке  XN), а также через конечные расности правой части системы в точке  XN.

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

Третья формула является экстраполяционной формулой, применяемой для предсказания приближенного значения решения в используемом для интегрирования методе Штермера.

При этом приближенные значения вычисляются через значения его разностной производной назад.

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

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

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

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

    int de38r_c (S_fp f, integer *m, integer *iorder, real *xn,
                 real *yn, real *dyn, real *hmin, real *eps, real *p, real *h,
                 real *x, real *yx, real *dyx, real *z, real *df, real *rfn,
                 real *rf, real *r, integer *ierr)

   Параметры 

   Версии 

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

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

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

int main(void)
{
    /* Local variables */
    extern int de38r_c(U_fp, int *, int *, float *, float *, float *,
                       float *, float *, float *, float *, float *, float *,
                       float *, float *, float *, float *, float *, float *,
                       int *);
    static float hmin;
    static int ierr;
    static float h__;
    static int i__;
    static float p, r__[2], x, z__[2];
    extern int f2_c();
    static float df[10] /* was [2][5] */, rf[2], xn, yn[2], yx[2], rfn[2],
                 eps, dyn[2], dyx[2];

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

    yn[0] = 1.f;
    yn[1] = -1.f;
    dyn[0] = 1.5f;
    dyn[1] = -.5f;
    xn = 0.f;
    p = 100.f;
    h__ = .01f;
    hmin = 1e-11f;
    eps = 1e-4f;
    de38r_c((U_fp)f2_c, &c__2, &c__5, &xn, yn, dyn, &hmin, &eps, &p, &h__, &x,
            yx, dyx, z__, df, rfn, rf, r__, &ierr);

    printf("\n %5i \n", ierr);
    printf("\n %16.7e \n", x);
    printf("\n %16.7e %16.7e \n", yx[0], yx[1]);
    printf("\n %16.7e %16.7e \n", dyx[0], dyx[1]);
    printf("\n %16.7e %16.7e \n", z__[0], z__[1]);
    for (i__ = 1; i__ <= 2; ++i__) {
         printf("\n %16.7e %16.7e %16.7e ",
                df_ref(i__, 1), df_ref(i__, 2), df_ref(i__, 3));
         printf("\n %16.7e %16.7e \n",
                df_ref(i__, 4), df_ref(i__, 5));
    }
    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

          x = 1.000000000001 - 02

          yx(1) = 1.014949833751                yx(2) = - 1.004949833752

          dyx(1) = 1.489950167079              dyx(2) = - 4.899501670834 - 01

          z__(1) = 1.494983375078                    z__(2) = - 4.949833750843 - 01

          df_ref(1, 1) = - 8.416248614562 - 04    df_ref(1, 2) = - 8.291528123650 - 06
          df_ref(1, 3) =   8.333255330228 - 08     df_ref(1, 4) =   8.377316618179 - 10
          df_ref(1, 5) = - 8.498979298110 - 12

          df_ref(2, 1) =   8.416248614518 - 04    df_ref(2, 2) =   8.291528119209 - 06
          df_ref(2, 3) = - 8.333252932146 - 08    df_ref(2, 4) = - 8.377654125979 - 10
          df_ref(2, 5) =   8.512301974406 - 12

 Точные значения решения и производной такие:

    y1  =   1.014949833745   ,    y2  =  - 1.004949833748
    y1' =   1.489950167079   ,    y2' =  - 4.899501670798 - 01