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

Назначение

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

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

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

          Y' = F (X, Y) ,

          Y = ( y1, ..., yM ) ,
          F = ( f1 (X, y1,..., yM),..., fM (X, y1,..., yM) )
 с начальными условиями, заданными в точке XN:
          Y(XN) = YN ,     YN = ( y10,..., yM0 ) , 

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

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

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

                M   
             (  ∑ δI2 )1/2  ≤  EPS  ,
                I=1 

где δ_I - одна из погрешностей следующих типов: абсолютная, относительная или стандартная. При этом под относительной погрешностью приближенного значения I - й компоненты решения в узле  x_n подразумевается отношение абсолютной погрешности  e_I этого значения в узле  x_n к абсолютной величине значения I - й компоненты в предыдущем узле  x_n - 1, т.е.  e_I / | y_I^n - 1|, а под стандартной погрешностью - отношение  e_I / YPM (I), где

           YPM(I)  =  max ( | yI0 |, | yI1 |,..., | yIn-1 | ) 

Tип погрешности специфицируется пользователем при обращении к подпрограмме.

Gear C.W. The automatic integration of ordinary differential equations. Communicatuons 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 de23r_c (real *f, integer *m, real *xn, real *yn, real *xk,
                 real *hmin, real *hmax, real *eps, integer *istifj, integer *iorder, 
                 integer *iu, real *h, real *y, real *ypm, real *delty, real *rab, 
                 real *yp, integer *ierr)

   Параметры 

   Версии 

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

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

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

           y1'  =  -20 y1  +  y2
            y2'  =   19 y1 - 2 y2 ,     0 ≤ x ≤ 1

           y1(0)  =  2 ,      y2(0)  =  18
 
 Собственные значения якобиевой матрицы системы
   λ1  =  -21 ,    λ2  =  -1 .
 Поэтому эту систему можно рассматривать как жесткую.

int main(void)
{
    /* Local variables */
    extern int de23r_c(U_fp, int *, float *, float *, float *, float *,
                       float *, float *, int *, int *, int *, 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 x, y[2], delty[2];
    static int iu;
    static float xk, xn, yn[2], yp[16];
    static int iorder, istifj;
    static float rab[38], eps, ypm[2];

    iu = 3;
    eps = 1e-5f;
    iorder = 6;
    hmin = 1e-10f;
    m = 2;
    xn = 0.f;
    yn[0] = 2.f;
    yn[1] = 18.f;
    xk = 1.f;
    istifj = 1;
    hmax = .1f;
    h__ = .01f;
    de23r_c((U_fp)f_c, &m, &xn, yn, &xk, &hmin, &hmax, &eps, &istifj, &iorder,
            &iu, &h__, y, ypm, delty, rab, yp, &ierr);

    printf("\n %16.7e \n", x);
    printf("\n %16.7e %16.7e \n", y[0], y[1]);
    return 0;
} /* main */

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

    /* Function Body */
    dy[1] = y[1] * -20.f + y[2];
    dy[2] = y[1] * 19.f - y[2] * 2.f;
    return 0;
} /* f_c */


Результаты:

          ierr  =  0

          y(1)  =  3.678796482-01
          y(2)  =  6.989709436