Текст подпрограммы и версий
de23r_c.zip , de23d_c.zip , de25r_c.zip , de25d_c.zip
Тексты тестовых примеров
tde23r_c.zip , tde23d_c.zip , tde25r_c.zip , tde25d_c.zip

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

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

где δI - одна из погрешностей следующих типов: абсолютная, относительная или стандартная. При этом под относительной погрешностью приближенного значения I - й компоненты решения в узле  xn подразумевается отношение абсолютной погрешности  eI этого значения в узле  xn к абсолютной величине значения I - й компоненты в предыдущем узле  xn - 1, т.е.  eI / | yIn - 1|, а под стандартной погрешностью - отношение  eI / 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)

Параметры

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 (тип: вещественный);
hmin - минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений; это значение должно быть много меньше среднего ожидаемого шага интегрирования, задаваемого параметром h (тип: вещественный);
hmax - максимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
eps - допустимая погрешность, с которой требуется вычислить все компоненты решения; тип погрешности специфицируется с помощью параметpа iu (тип: вещественный);
istifj - целый указатель метода численного интегрирования:
istifj=0 - интегрирование системы ведется методом Адамса;
istifj=1 - интегрирование ведется специальным методом, предназначенным для жестких систем;
iorder - целая переменная, указывающая максимальный допустимый порядок метода; iorder должен быть не больше 7 для метода Адамса и не больше 6 для метода интегрирования жестких систем;
iu - целый указатель типа погрешности численного решения:
iu = 1 - для стандартной погрешности;
iu = 2 - для относительной погрешности;
iu = 3 - для абсолютной погрешности;
h - вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления итегрирования, т.е. положительным, если xk > xn, отрицательным, если xk < xn, или без такого учета в виде абсолютной величины;
y - искомое решение задачи Коши, вычисленное подпрограммой для значения аргумента xk; для системы уравнений (когда m ≠ 1) задается массивом длиной m; в случае совпадения значений параметров xn и xk значение y полагается равным начальному значению yn (тип: вещественный);
         ypm -
       delty  
одномерные вещественные рабочие массивы длиной m;
rab - одномерный вещественный рабочий массив; при интегрировании нежесткой системы уравнений rab имеет размер 17*m, при интегрировании жесткой системы - m*(m + 17);
yp - двумерный вещественный рабочий массив размеpа 8*m;
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом:
ierr= 1 - когда неправильно задан параметр iorder, а именно, когда iorder превосходит максимальный допустимый порядок метода; в этом случае интегрирование системы ведется методом Гира порядка не выше 7 для нежесткой системы, и не выше 6 для жесткой;
ierr=65 - когда решение системы не может быть вычислено с требуемой точностью eps при заданных начальном шаге h, его минимальном значении hmin и порядке метода iorder;
ierr=66 - когда приближенное значение решения не может быть вычислено, т.к. итерационный процесс его определения не сходится для шагов интегрирования h, больших заданного минимального значения hmin;
ierr=67 - когда требуемая точность eps вычисления приближенного решения меньше той, которая может быть достигнута для данной задачи при тех размерах шага интегрирования, начальное значение которого задано параметром h;
  при ierr = 65, 66, 67 интегрирование системы прекращается; при желании интегрирование можно повторить обращением к подпрограмме с новыми значениями параметров h, hmin и iorder;
ierr=68 - когда приближенное значение решения для жесткой системы не может быть вычислено с заданной точностью; для достижения тебуемой точности следует воспользоваться версиями подпрограммы de23d_c, de25r_c или de25d_c.

Версии

de23d_c - вычисление решения задачи Коши для нежесткой и жесткой систем обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Гира с повышенной точностью. При этом параметры xn, yn, xk, hmin, hmax, eps, h, y, ypm, delty, rab, yp и параметры x, y и dy в подпрограмме f должны иметь тип double.
de25r_c - вычисление решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Гира с автоматическим выбором шага. Первый оператор подпрограммы имеет вид:
int de25r_c(real *f, real *fj, integer *m, real *xn, real *yn, real *xk, real *hmin, real *hmax, real *eps, integer *iorder, integer *iu, real *h, real *y, real *ypm, real *delty, real *rab, real *yp, integer *ierr)
Здесь: fj - имя подпрограммы вычисления якобиана правой части системы; первый оператор этой подпрограммы имеет вид:
int fj (float *x, float *y, float *z, int *m)
Здесь: x, y - значения независимой и зависимой переменных, соответственно, причем y представляет одномерный массив длины m; вычисленное значение якобиана должно быть помещено в двумерный массив z размера m*m, при этом частная производная от правой части i - ого уравнения по j - ой переменной y (j) запоминается в элементе z (i, j) (тип параметров x, y и z: вещественный). Остальные параметры подпрограммы de25r_c имеют тот же смысл, что и одноименные параметры подпрограммы de23r_c.
de25d_c - вычисление решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Гира с повышенной точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме de25r_c; при этом параметры xn, yn, xk, hmin, hmax, eps, h, y, ypm, delty, rab, yp и параметры x, y, dy и z в подпрограммах f и fj должны иметь тип double.

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

de21r_c - выполнение одного шага численного интегрирования нежесткой и жесткой систем обыкновенных дифференциальных уравнений первого порядка методом Гира с контролем точности; вызывается при работе подпрограммы de23r_c.
de21d_c - выполнение одного шага численного интегрирования нежесткой и жесткой систем обыкновенных дифференциальных уравнений первого порядка методом Гира с повышенной точностью; вызывается при работе подпрограммы de23d_c.
de24r_c - выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом Гира с контролем точности; вызывается при работе подпрограммы de25r_c.
de24d_c - выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом Гира с повышенной точностью; вызывается при работе подпрограммы de25d_c.
utde12_c - подпрограмма выдачи диагностических сообщений при работе подпрограмм de21r_c, de23r_c, de24r_c, de25r_c.
utde13_c - подпрограмма выдачи диагностических сообщений при работе подпрограмм de21d_c, de23d_c, de24d_c, de25d_c.

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

 

B общем случае заданая точность eps не гарантируется.

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

Значение hmin должно быть много меньше среднего ожидаемого шага интегрирования, задаваемого параметром h при обращении к подпрограмме и ее версиям, т.к. интегрирование системы начинается с метода первого порядка.

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

           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