Текст подпрограммы и версий
de20r_c.zip , de20d_c.zip 		Тексты тестовых примеров
tde20r_c.zip , tde20d_c.zip

Подпрограмма:  de20r_c

Назначение

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

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

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

          Y ' = F (X, Y) ,
          Y = ( y1, ... , yM ) ,
          F = ( f1 (X, y1, ... , yM), ... , fM (X, y1, ... , yM) )

методом рациональной экстраполяции Грэгга - Булирша - Штера. По заданному значению решения YX в узле Xn вычисляется значение этого решения в узле Xn + H. Bсе компоненты решения проверяются на точность по одному из следующих типов погрешности: относительной, абсолютной и стандартной. Тип погрешности специфицируется пользователем при обращении к подпрограмме.

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

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

Значение H может быть меньше или pавно значению шага интегрирования, задаваемому пользователем.

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

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

    int de20r_c (S_fp f, integer *m, integer *iorder, integer *iu,
                 integer *jstart, real *hmin, real *eps, real *x, real *yx, real *h, 
                 logical *bul, 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 - количество уравнений в системе (тип: целый);
iorder - целая переменная, указывающая максимальный допустимый порядок метода рациональной экстраполяции; iorder должен быть меньше 7;
iu - целый указатель типа погрешности численного решения:
iu = 1 - для стандартной погрешности;
iu = 2 - для относительной погрешности;
iu = 3 - для абсолютной погрешности;
jstart - целый указатель режима использования подпрограммы:
jstart = 0 - выполнить один шаг интегрирования дифференциального уравнения для значений независимой и зависимой переменных и шага интегрирования, заданных параметрами  x, yx и h, соответственно; самое первое обращение к подпрограмме должно быть выполнено со значением  jstart = 0 ;
jstart= -1 - повторить последний шаг интегрирования с новыми значениями параметpов  h, iorder и hmin; при этом значения параметров  x, yx и ypm будут взяты равными значениям из последнего обращения к подпрограмме со значением  jstart = 0 ;
hmin - минимальное значение абсолютной величины шага, который разрешается использовать при интегрированиии данной системы уравнений (тип: вещественный);
eps - допустимая погрешность, с которой требуется вычислить все компоненты решения; тип погрешности специфицируется с помощью параметра  iu (тип: вещественный);
x, yx - заданные вещественные значения аргумента и pешения уравнения в нем. B результате работы подпрограммы в x получается новое значение аpгумента, отстоящее от исходного на величину шага интегрирования, а в yx - соответствующее значение решения. B случае системы уравнений, т.е. когда m ≠ 1 , yx задается одномерным массивом длины m;
h - вещественная переменная, содержащая значение шага интегрирования. Если для этого значения шага точность приближенного решения достигается, то именно он реализуется подпрограммой, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность. После выполнения такого шага подпрограмма определяет, можно ли достичь заданную точность с большим значением шага интегрирования. Если это возможно, то это большое значение шага будет помещено в h. Таким образом, в pезультате работы подпрограммы в h находится новое значение шага интегрирования, котоpое может быть больше, меньше или pавно исходному значению.
bul - логическая переменная, значение которой при обращении к подпрограмме полагается равным TRUE_, если заданный в h шаг выводит в конец интервала интегрирования, и FALSE_ в противном случае. B результате работы подпрограммы bul pавно FALSE_, если вместо исходного шага интегрирования был реализован меньший шаг; в противном случае, т.е. когда был выполнен именно заданный при обращении в h шаг, значение переменной bul не меняется;
ypm - одномерный вещественный массив длины m,  i - ый элемент ypm(i) которого содержит:
iu = 1 - максимальное значение абсолютной величины i - ой компоненты решения, вычисленное с начала интегрирования;
iu = 2 - абсолютную величину значения i - ой компоненты yx(i) решения, заданного в yx;
  B результате работы подпрограммы значение ypm сохраняется, если вновь вычисленное значение решения не превосходит по абсолютной величине исходного значения ypm, и заменяется на абсолютную величину нового значения, если она больше первоначального значения ypm. Cамое первое обращение к подпрограмме должно быть выполнено со значением ypm(i) = | yx(i) | ,  i = 1,..., m. При  iu = 3 все элементы массива ypm должны быть равны 1; в этом случае значение параметра ypm сохраняется;
delty - одномерный вещественный массив, i - ый элемент которого в результате работы подпрограммы содержит абсолютную погрешность для i - ой компоненты решения, с которой она вычисляется на текущем шаге;
rab - одномерный вещественный рабочий массив длины 29*m;
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом:
ierr= 1 - когда неправильно задан параметр iorder: iorder < 1 или iorder > 6 (максимальный допустимый порядок метода); в этом случае выполняется один шаг интегрирования системы уравнений со значением iorder = 6;
ierr=65 - когда какая - нибудь компонента решения не может быть вычислена с требуемой точностью eps; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметpов  h, hmin и iorder и значением  jstart = - 1. 
   Версии
de20d_c - выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений первого порядка методом рациональной экстраполяции Грэгга - Булирша - Штера с повышенной точностью. При этом параметры hmin, eps, x, yx, h, ypm, delty, rab и параметры x, y и dy в подпрограмме f должны иметь тип double. 
   Вызываемые подпрограммы
utde10_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de20r_c.
utde11_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de20d_c. 
   Замечания по использованию
 

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

Значения параметров x, yx и h изменяются всегда, независимо от того, будет ли достигнута заданная точность или нет.

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

При многократном использовании подпрограммы для вычисления решения системы уравнений на отрезке значение паpаметpа ypm задается только один раз, а именно, при самом первом обращении к подпрограмме. При этом:
ypm(i) = | yx(i) | ,  i = 1,..., m  для  iu = 1 и iu = 2;
ypm(i) = 1 ,  i = 1,..., m  для  iu = 3.

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

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

          y1'  =   y3
          y2'  =   y4
          y3'  =  - (1 - 0.002 x) y1 / ( y12  +  y22 )3/2
          y4'  =  - (1 - 0.002 x) y2 / ( y12  +  y22 )3/2

          y1 (0)  =  0.09411764706 ,    y2 (0)  =  0 ,   y3 (0)  =  0 ,   y4 (0)  =  4.5

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

int main(void)
{
    /* Local variables */
    extern int de20r_c(U_fp, int *, int *, int *, int *, float *, float *,
                       float *, float *, float *, logical *, 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, jstart;
    static float rab[116];
    static logical bul;
    static float eps, ypm[4];

    m = 4;
    xn = 0.f;
    xk = 18.849555921540002f;
    h__ = .01f;
    iorder = 6;
    jstart = 0;
    y[0] = .09411764706f;
    y[1] = 0.f;
    y[2] = 0.f;
    y[3] = 4.5f;
    ypm[0] = y[0];
    ypm[1] = y[1];
    ypm[2] = y[2];
    ypm[3] = y[3];
    eps = .001f;
    hmin = 1e-10f;
    iu = 2;
    bul = FALSE_;
    de20r_c((U_fp)f_c, &m, &iorder, &iu, &jstart, &hmin, &eps, &xn, y, &h__,
            &bul, 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 %16.7e \n", h__);
    printf("\n %16.7e %16.7e %16.7e %16.7e \n",
           ypm[0], ypm[1], ypm[2], ypm[3]);
    printf("\n %16.7e %16.7e %16.7e %16.7e \n",
           delty[0], delty[1], delty[2], delty[3]);
    printf("\n %5i %5i \n", jstart, ierr);
    return 0;
} /* main */

int f_c(float *t, float *y, float *z, int *m)
{
    /* Builtin functions */
    double sqrt(double);

    /* Local variables */
    static float r__, r2;

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

    /* Function Body */
    r2 = y[1] * y[1] + y[2] * y[2];
    r2 *= (float)sqrt(r2);
    r__ = -(1.f - *t * .002f) / r2;
    z__[1] = y[3];
    z__[2] = y[4];
    z__[3] = r__ * y[1];
    z__[4] = r__ * y[2];
    return 0;
} /* f_c */


Результаты:
          
          xn  =  0.01
          y(1)  =   0.08867041316     y(2)  =   0.04414645692
          y(3)  =  -1.052309920         y(4)  =   4.252530734

          h  =  0.13333334
          ypm(1)  =  0.08867041343    ypm(2)  =  0.04414645679
          ypm(3)  =  1.052309916        ypm(4)  =  4.252530728

          delty(1)  =  0.2673914423 10-9
          delty(2)  =  0.1339799383 10-9
          delty(3)  =  0.4718458513 10-8
          delty(4)  =  0.6097252481 10-8

          jstart  =  0
          ierr  =  0