Текст подпрограммы и версий
de27r_c.zip , de27d_c.zip
Тексты тестовых примеров
tde27r_c.zip , tde27d_c.zip

Подпрограмма:  de27r_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, которая является концом интервала интегрирования. Интегрирование ведется с постоянным шагом, значение которого задается пользователем при обращении к подпрограмме.

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

Беленький В.З. Стандартная программа для интегрирования системы дифференциальных уравнений методом Адамса. Сб. "Вычислительные методы и программирование", вып. 3. Изд-во МГУ, 1965

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

    int de27r_c (real *f, integer *m, real *xn, real *yn, real *xk,
                 real *h, real *y, real *delty, real *df, real *rf, real *yp,
                 real *ry, real *rfn, integer *ierr)

Параметры

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

Версии

de27d_c - вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Адамса с постоянным шагом, при этом все вычисления над действительными числами ведутся с удвоенным числом значащих цифр. B этом случае параметры xn, yn, xk, h, y, delty, df, rf, yp, ry, rfn и параметры x, y и dy в подпрограмме f должны иметь тип double.

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

de26r_c - выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений первого порядка методом Адамса без контроля точности; вызывается при работе подпрограммы de27r_c.
de26d_c - выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений первого порядка методом Адамса без контроля точности, при этом все вычисления над действительными числами ведутся с удвоенным числом значащих цифр; вызывается при работе подпрограммы de27d_c.
utde12_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de27r_c.
utde13_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de27d_c.

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

 

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

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

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

           y1'  =   y2
           y2'  =  -y1

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

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

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

    /* Local variables */
    extern int de27r_c(U_fp, int *, float *, float *, float *, float *,
                       float *, float *, float *, float *, float *, float *,
                       float *, int *);
    static int ierr;
    extern int f_c();
    static float h__;
    static int i__, m;
    static float y[2], delty[2], df[10], rf[2], xk, xn, yn[2], yp[2], ry[10],
               rfn[2];

    m = 2;
    xn = 2.3561944875000003f;
    yn[0] = (float)sqrt(2.f) / 2.f;
    yn[1] = -yn[0];
    xk = 3.14159265f;
    for (i__ = 1; i__ <= 2; ++i__) {
        h__ = .01f;
        de27r_c((U_fp)f_c, &m, &xn, yn, &xk, &h__, y, delty, df, rf, yp, ry,
                rfn, &ierr);

    printf("\n %15.7e \n", xk);
    printf("\n %15.7e %15.7e \n", y[0], y[1]);
        xn = xk;
        xk = xn * .75f;
        yn[0] = 0.f;
        yn[1] = -1.f;
/* l30: */
    }
    return 0;
} /* main */

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

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


Результаты:

    после первого обращения к подпрограмме -
          xk                               y(1)                         y(2)
          3.141592654 + 00      -2.773237299-10      -1.000000000 + 00
 
    после второго обращения к подпрограмме -
          xk                               y(1)                         y(2)
          2.356194490 + 00       7.071067811-01      -7.071067815-01