Текст подпрограммы и версий
de22r_c.zip , de22d_c.zip 		Тексты тестовых примеров
tde22r_c.zip , tde22d_c.zip

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

Пусть Uin и Vin являются двумя последовательными приближениями к значению точного решения Yin в узле  xn. Тогда Vin принимается за искомое значение решения в узле  xn, если выполнен один из критериев сходимости:

1.  | Uin - Vin | ≤ EPS * | Zin | , где Zin = max{ Yi1,Yi2,...,Yin - 1} для  i = 1,...,M (стандартная погрешность)
2.  | Uin - Vin | ≤ EPS * | Vin | , для  i = 1,...,M (относительная погрешность)
3.  | Uin - Vin | ≤ EPS  , для  i = 1,...,M (абсолютная погрешность).

Bulirsch R., Stoer J., Numerical treatment of ordinary differential equations by extrapolation methods. Numerische Mathematik, 8(1) 1966, 1-13.

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

    int de22r_c (real *f, integer *m, real *xn, real *yn, real *xk,
                 real *hmin, real *eps, integer *iorder, integer *iu, real *h,
                 real *y, real *ypm, real *delty, real *rab, integer *ierr)

   Параметры
f - имя подпрограммы вычисления значений правой части дифференциального уравнения. Первый оператоp подпрограммы должен иметь вид:
int f (float *t, float *y, float *dy, int *m).
Здесь: x, y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в dy. B случае системы уравнений, т.е. когда m ≠ 1, параметры y и dy представляют массивы длины m (тип параметров x, y и dy: вещественный);
m - количество уравнений в системе (тип: целый);
xn, yn - начальные значения аргумента и решения. B случае системы уравнений (т.е. m ≠ 1) yn представляет массив длины m (тип: вещественный);
xk - значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования). xk может быть больше, меньше или pавно xn (тип: вещественный);
hmin - минимальное значение абсолютной величины шага, который разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
eps - допустимая погрешность, с которой требуется вычислить все компоненты решения; тип погрешности специфицируется с помощью параметpа iu (тип: вещественный);
iorder - целая переменная, указывающая максимальный допустимый порядок метода рацональной экстраполяции; iorder должен быть меньше 7;
iu - целый указатель типа погрешности численного решения:
iu = 1 - для стандартной погрешности;
iu = 2 - для относительной погрешности;
iu = 3 - для абсолютной погрешности;
h - вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если xk > xn, отрицательным, если xk < xn, или без такого учета в виде абсолютной величины;
y - искомое решение задачи Коши, вычисленное подпрограммой для значения аргумента xk. Для системы уравнений (когда m ≠ 1) y задается массивом длины m. B случае совпадения значений параметров xn и xk значение y полагается равным начальному значению yn (тип: вещественный);
         ypm -
       delty   		одномерные вещественные рабочие массивы длины m;
rab - одномерный вещественный рабочий массив длины 29*m;
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом:
ierr= 1 - когда неправильно задан параметр iorder, т.е.iorder < 1 или iorder > 6 максимальный допустимый порядок метода); в этом случае интегрирование данной системы уравнений ведется со значением iorder=6;
ierr=65 - когда какая - нибудь компонента pешения не может быть вычислена с требуемой точностью eps; в этом случае интегрирование прекращается. При желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров h, hmin и iorder. 
   Версии
de22d_c - вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом рациональной экстраполяции Грэгга - Булирша - Штера с повышенной точностью. При этом параметры xn, yn, xk, hmin, eps, h, y, ypm, delty, rab и параметры x, y и dy в подпрограмме f должны иметь тип double. 
   Вызываемые подпрограммы
de20r_c - выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений первого порядка методом рациональной экстраполяции Грэгга - Булирша - Штера; вызывается при работе подпрограммы de22r_c.
de20d_c - выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений первого порядка методом рацональной экстраполяции Грэгга - Булирша - Штера с повышенной точностью; вызывается при работе подпрограммы de22d_c.
utde10_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de22r_c.
utde11_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de22d_c. 
   Замечания по использованию
 

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

При работе подпрограммы значения параметров m, xn, yn, xk, hmin, eps, iorder, iu сохраняются. Если после работы подпрограммы нет необходимости иметь начальное значение решения yn, то параметры yn и y при обращении к ней можно совместить.

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

Использование подпрограммы иллюстрируется на примере: 

           y1'  =  0.2 ( y4 - y1 )
           y2'  =  y1  +  2 ( y2 - y2 y3 )
           y3'  =  y4  -  ( y3 - y2 y3 )
           y4'  =  10 y1 -  ( 61 - 0.13 x ) y4  +  0.13 x ,    0 ≤ x ≤ 8

          y1(0)  =  y2(0)  =  y3(0)  =  y4(0)  =  0 .

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

int main(void)
{
    /* Local variables */
    extern int de22r_c(U_fp, int *, float *, float *, float *, float *,
                       float *, int *, int *, float *, float *, float *,
                       float *, float *, int *);
    static float hmin;
    static int ierr;
    extern int f_c();
    static float h__;
    static int m;
    static float y[4], delty[4];
    static int iu;
    static float xk, xn;
    static int iorder;
    static float rab[116], eps, ypm[4];

    m = 4;
    y[0] = 0.f;
    y[1] = 0.f;
    y[2] = 0.f;
    y[3] = 0.f;
    xn = 0.f;
    xk = 8.f;
    h__ = .01f;
    hmin = 1e-6f;
    eps = 1e-5f;
    iorder = 4;
    iu = 2;
    de22r_c((U_fp)f_c, &m, &xn, y, &xk, &hmin, &eps, &iorder, &iu, &h__, y,
            ypm, delty, rab, &ierr);

    printf("\n %16.7e \n", xn);
    printf("\n %16.7e %16.7e %16.7e %16.7e \n", y[0], y[1], y[2], y[3]);
    printf("\n %5i \n", ierr);
    return 0;
} /* main */

int f_c(float *t, float *y, float *z, int *m)
{
    static float r23, ct;

    /* Parameter adjustments */
    --z__;
    --y;

    /* Function Body */
    r23 = y[2] * y[3];
    ct = *t * .13f;
    z__[1] = (y[4] - y[1]) * .2f;
    z__[2] = y[1] + (y[2] - r23) * 2.f;
    z__[3] = y[4] - (y[3] - r23);
    z__[4] = y[1] * 10.f - (61.f - ct) * y[4] + ct;
    return 0;
} /* f_c */


Результаты:

          y(1)  =  0.9238470846 * 10-2
          y(2)  =  0.4820971549 * 10-2
          y(3)  =  1.667112950
          y(4)  =  0.188434937 * 10-1

          ierr  =  0