Текст подпрограммы и версий
de41r_c.zip , de41d_c.zip 		Тексты тестовых примеров
tde41r_c.zip , tde41d_c.zip

Подпрограмма:  de41r_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 )   ,
      Y ' (XN) = DYN   ,    DYN = ( y10', ..., yM0' )   ,   -

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

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

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

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

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

    int de41r_c (real *f, integer *m, real *xn, real *yn,
                 real *dyn, real *xk, real *h, real *y, real *dy, real *delty,
                 real *df, real *rfn, real *rf, real *yp, real *dyp, integer *ierr)

   Параметры
f - имя подпрограммы вычисления значений правой части системы дифференциальных уpавнений. Первый оператор подпрограммы должен иметь вид:

int f (float *x, float *y, float *dy, int *m)

Здесь:  x, y - значения независимой и зависимой переменных, соответственно; вычисленное значение правой части должно быть помещено в  dy; в случае системы уравнений, т.е. когда  m ≠ 1, параметры  y и  dy представляют одномерные массивы длиной  m (тип параметров  x, y и  dy: вещественный);

m - количество уравнений в системе (тип: целый);
         xn, yn -
         dyn   		начальные значения аргумента, решения и его производной, соответственно; в случае системы уравнений (т.е.  m ≠ 1)  yn и  dyn представляют одномерные массивы длиной  m (тип: вещественный);
xk - значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования).  xk может быть больше, меньше или pавно  xn (тип: вещественный);
h - вещественная переменная, содержащая значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если  xk > xn, отрицательным, если  xk < xn, или без такого учета в виде абсолютной величины;
y - искомое решение задачи Коши, вычисленное подпрограммой для значения аргумента  xk. Для системы уравнений (когда  m ≠ 1)  y задается массивом длиной  m. b случае совпадения значений параметров  xn и  xk значение   y полагается равным начальному значению  yn (тип: вещественный);
            dy -
       delty  
       rfn, rf  
         yp, dyp   		вещественные одномерные рабочие массивы длиной  m;
df - двумерный вещественный рабочий массив размеpа  m * 5;
ierr - целая переменная, значение которой в pезультате работы подпрограммы полагается равным 65, если при заданном значении шага интегрирования не может быть достигнута сходимость численного решения к точному значению решения; в этом случае интегрирование системы можно повторить обращением к подпрограмме с новым (меньшим) значением параметра  h. 
   Версии
de41d_c - вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений второго порядка в конце интервала интегрирования методом Штермера с постоянным шагом, при этом все вычисления над действительными числами ведутся с удвоенным числом значащих цифр. B этом случае параметры  xn, yn, dyn, xk, h, y, dy, delty, df, rfn, rf, yp, dyp и параметры  x, y и  dy в подпрограмме  f должны иметь тип double. 
   Вызываемые подпрограммы
de40r_c - выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений второго порядка методом Штермера без контроля точности; вызывается при работе подпрограммы de41r_c.
de40d_c - выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений второго порядка методом Штермера без контроля точности, при этом все вычисления над действительными числами ведутся с удвоенным числом значащих цифр; вызывается при работе подпрограммы de41d_c.
utde12_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de41r_c.
utde13_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de41d_c. 
   Замечания по использованию
 

При работе подпрограммы и ее версии значения параметров   m, xn, yn, dyn, xk и  h сохраняются.

Если после работы подпрограммы нет необходимости иметь начальные значения решения  yn и его производной  dyn, то параметры  yn и  y, а в случае производной  -  параметры   dyn и  dy, при обращении к ней можно совместить.

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

     y1'' = - 6y1 - 7y2   ,
     y2'' = - 3y1 - 2y2 + 2x   ,
     y1 (0) = 0   ,   y2 (0) = 0   ,   y1' (0) = - 1 / 9   ,   y2' (0) = - 3   .

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

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

    /* Local variables */
    extern int de41r_c(U_fp, int *, float *, float *, float *, float *,
                       float *, float *, float *, float *, float *, float *,
                       float *, float *, float *, int *);
    static int ierr;
    extern int f_c();
    static float h__;
    static int k, m;
    static float y[2], delty[2], y1, y2, y3, df[10], rf[2], dy[2],
                  xk, xn, yn[2], yp[2], rfn[2], dyn[2], dyp[2];

    h__ = .1f;
    m = 2;
    xn = 0.f;
    xk = 5.f;
    yn[0] = 0.f;
    yn[1] = 0.f;
    dyn[0] = -.1111111111111111f;
    dyn[1] = -3.f;
    y3 = ((float)exp(5.f) - (float)exp(-5.f)) / 3.f;
    y1 = y3 - (float)sin(15.f) * 7.f / 9.f + 7.7777777777777777f;
    y2 = -y3 - (float)sin(15.f) / 3.f - 6.666666666666667f;
    for (k = 1; k <= 2; ++k) {
        de41r_c((U_fp)f_c, &m, &xn, yn, dyn, &xk, &h__, y, dy, delty, df, rfn,
                rf, yp, dyp, &ierr);

        printf("\n %16.7e \n", xk);
        printf("\n %16.7e %16.7e \n", y[0], y[1]);
        printf("\n %16.7e \n", h__);
        printf("\n %16.7e %16.7e \n", y1, y2);
        xk = -5.f;
        y1 = -y1;
        y2 = -y2;
        h__ = .1f;
/* L200: */
    }
    return 0;
} /* main */

int f_c(float *t, float *y, float *z, int *m)
{
    /* Parameter adjustments */
    --z__;
    --y;

    /* Function Body */
    z__[1] = y[1] * -6. - y[2] * 7.;
    z__[2] = y[1] * -3. - y[2] * 2. + *t * 2.;
    return 0;
} /* f_c */


Результаты:

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

                   xk                           y(1)                             y(2)
       5.000000000 + 00      5.674080692 + 01     - 5.635223804 + 01

                   h__                           y1                                 y2
       1.000000000 - 02      5.674080540 + 01     - 5.635223633 + 01

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

                   xk                           y(1)                             y(2)
     - 5.000000000 + 00    - 5.674080382 + 01       5.635223450 + 01

                   h__                         y1                                 y2
     - 1.000000000 - 02    - 5.674080540 + 01       5.635223633 + 01