Текст подпрограммы и версий
de00r_c.zip , de00d_c.zip , de84r_c.zip , de84d_c.zip
Тексты тестовых примеров
tde00r_c.zip , tde00d_c.zip , tde84r_c.zip , tde84d_c.zip

Подпрограмма:  de00r_c (версия de84r_c)

Назначение

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

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

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

(1)                             Y '  =   F(X, Y) ,

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

Предполагается, что среди характеристических корней матрицы Якоби ∂F / ∂y функции F имеются большие по модулю корни. По заданному значению решения YX в узле  xn  вычисляется значение этого решения в узле  xn + H . Предполагается также, что система (1) на [ xn, xn + H ] может быть достаточно хорошо аппроксимирована линейной однородной системой с постоянными коэффициентами, т.е. системой вида

                                     Y '  =   AY  , 

где A - некоторая постоянная матрица. Вычисление производится по методу Лоусона. Метод Лоусона заключается в следующем. Исходная система уравнений с помощью замены искомой функции  y (x) на [ xn, xn + H ] по формуле

                             y(x)  =  exp [ ( x - xn ) A0 ] Z(x)  ,

где A0 - некоторая постоянная матрица, преобразуется в систему уравнений относительно новой неизвестной функции Z (X) :

(2) Z ' (x)  =  F1(x, z)  =  exp [ - ( x - xn ) A0 ] { F(x, exp [ ( x - xn ) A0 ] Z(x)) -
                                      - A0 exp [ ( x - xn ) A0 ] Z(x) } 

Матрица Якоби ∂ F1 / ∂Z системы (2) и матрица Якоби ∂F / ∂y системы (1) связаны между собой соотношением

         ∂ F1 / ∂Z  =  exp [ - (x - xn) A0] { ∂F / ∂y - A0 } exp [ ( x - xn ) A0 ] 

Указанное преобразование выполняется самой подпрограммой. В качестве матрицы A0 подпрограмма выбирает матрицу ∂F / ∂y , вычисленную либо в точке ( xn, yn ), либо в некоторой предыдущей точке ( xn - k, yn - k ). Если шаг H достаточно мал, то данное преобразование позволяет уменьшить характеристические корни матрицы ∂ F1 / ∂Z по сравнению с характеристическими корнями матрицы ∂F / ∂y . А это приводит к уменьшению константы Липшица системы (2) по сравнению с константой Липшица системы (1).

Для решения системы (2) применяются классический метод Рунге - Кутта четвертого порядка точности, причем одновременно с решением (2) производится обратное преобразование от функции Z (x) к функции Y (x).

Все компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Абсолютная погрешность приближенного решения оценивается по правилу Рунге.

J.Douglas Lowson, Generalized Runge - Kutta processes for stable systems with large lipshitz constants, SIAM Journal on Numerical Analisys, 1967, Vol.4, No 3.

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

    int de00r_c (real *f, S_fp fj, integer *m, integer *jstart,
                 real *hmin, real *eps, real *p, real *yx, real *x, real *h,
                 logical *bul, real *xp, real *yp, real *a, real *e1, real *e2,
                 real *e3, real *e4, real *r, real *r0, real *r1, real *r2,
                 real *r3, real *yd, integer *ierr)

Параметры

f - подпрограмма вычисления правой части системы в точке  x . Первый оператор подпрограммы должен иметь вид:
int f (float *x, float *y, float *dy, int *m).
Здесь x и y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в dy. В случае системы уравнений, т.е. m ≠ 1, параметры y и dy представляют одномерные массивы длиной m (тип параметров x, y и dy: вещественный);
fj - подпрограмма вычисления матрицы Якоби правой части системы. Первый оператор подпрограммы должен иметь вид:
int fj (float *x, float *y, float *df, int *m).
Здесь x и y - значения независимой и зависимой переменных, соответственно, причем y представляет одномерный массив длины m. df - двумерный массив размера m*m. Вычисленное значение матрицы Якоби ∂F / ∂y должно быть помещено в массив df, при этом частная производная от правой части i - ого уравнения F по j - ой переменной y (j) запоминается в dF (i, j), т.е. dF (i, j) = ∂Fi / ∂yj (тип параметров x, y и df: вещественный);
m - количество уравнений в системе (тип: целый);
jstart - целый указатель режима использования подпрограммы, имеющий следующие значения:
0 - первое обращение к подпрограмме должно быть выполнено с нулевым значением jstart; при этом значение параметра в качестве матрицы A0 принимается матрица ∂F / ∂y в точке [ x = xn, y = yn ]; вычисленное значение матрицы A0 будет использоваться в дальнейшем для значений jstart = 1, 2, 3;
1 - выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных, заданных параметрами yx и x соответственно; при этом само значение величины шага интегрирования H, на который предполагается продолжить решение системы, должно совпадать с тем значением шага, который был фактически выполнен подпрограммой при предыдущем обращении к ней, а матрица A0 остается той же, что и на предыдущем шаге;
2 - то же, что и для jstart = 1, но только с той разницей, что величина шага интегрирования, на который предполагается продолжить решение, в два раза больше того значения шага, который был фактически выполнен при предыдущем обращении к подпрограмме; матрица A0 остается той же, что и на предыдущем шаге;
3 - то же, что и для jstart = 1, но только с той разницей, что шаг интегрирования может принимать произвольное значение, матрица A0 остается той же, что и на предыдущем шаге.
-1 - повторить последний шаг интегрирования с новыми значениями параметров h и/или hmin; при этом значении jstart в качестве матрицы A0 принимается матрица ∂F / ∂y в точке [ x = xp, y = yp ] таким образом, до тех пор, пока jstart принимает значения 1, 2 или 3, сохраняется одно и то же значение матрицы A0 , равное значению матрицы Якоби, вычисленному при последнем обращении к подпрограмме с параметром jstart = 0 или jstart = - 1; если требуется изменить матрицу A0 , то следует обратиться к подпрограмме со значением jstart = 0 или - 1; при выходе из подпрограммы параметр jstart полагается равным 2, если выполнив заданный в h шаг, подпрограмма рекомендовала для интегрирования вдвое большеe его значение, и 1 в противном случае, т.е. в том случае, когда рекомендованное значение шага равно только что выполненному шагу; рекомендованное значение шага на выходе из подпрограммы запоминается в переменной h;
hmin - минимальное значение абсолютной величины шага, которое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
eps - допустимая мера погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный);
p - граница перехода, используемая при оценке меры погрешности решения (тип: вещественный);
yx, x - заданные вещественные значения решения и соответствующее ему значение аргумента; в результате работы подпрограммы в x получается новое значение аргумента, в yx - соответствующее значение решения; в случае системы уравнений т.е. когда m ≠ 1, yx задается одномерным массивом длины m;
h - вещественная переменная, содержащая значение шага интегрирования; если для этого значения шага точность приближенного решения достигается, то именно он и реализуется подпрограммой, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность eps; на выходе из подпрограммы h содержит рекомендуемое подпрограммой значение следующего шага, определяемое ею с целью достижения более экономного способа интегрирования;
bul - логическая переменная, значение которой при обращении к подпрограмме полагается равным TRUE_, если заданный в H шаг выводит в конец интервала интегрирования, и FALSE_ в противном случае; в результате работы подпрограммы bul равно FALSE_, если вместо исходного шага интегрирования был реализован меньший шаг; в противном случае, т.е. когда был выполнен именно заданный при обращении в h шаг, значение параметра bul не меняется;
xp, yp - вещественная рабочая переменная и одномерный рабочий массив длины m соответственно; значения параметров xp, yp на выходе из подпрограммы равны тем значениям, которые имели параметры x, yx при входе в нее (т.е. предыдущий узел и решение в нем);
     a, e1, e2 -
    e3, e4  
вещественные двумерные рабочие массивы размера m*m;
         r, r0 -
         r1, r2  
         r3, yd  
вещественные одномерные рабочие массивы длины m;
ierr - целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью eps; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров h, hmin и значением jstart = - 1.

Версии

de00d_c - выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона с удвоенным числом значащих цифр. При этом параметры hmin, eps, p, yx, x, h, xp, yp, a, e1, e2, e3, e4, r, r0, r1, r2, r3, yd и параметры x, y, dy в подпрограмме f и x, y, df в подпрограмме fj должны иметь тип double;
de84r_c - выполнение одного шага численного интегрирования устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона. В отличии от подпрограммы de00r_c первый оператор подпрограммы de84r_c не содержит в списке формальных параметров имени подпрограммы вычисления матрицы Якоби и имеет следующий вид:
int de84r_c (real *f, integer *m, integer *jstart, real *hmin, real *eps, real *p, real *yx, real *x, real *h, logical *bul, real *xp, real *yp, real *a, real *e1, real *e2, real *e3, real *e4, real *r, real *r0, real *r1, real *r2, real *r3, real *yd, integer *ierr).
При использовании этой подпрограммы пользователь не составляет подпрограмму вычисления матрицы Якоби. Все параметры подпрограммы de84r_c имеют тот же смысл, что и одноименные параметры подпрограммы de00r_c;
de84d_c - выполнение одного шага численного интегрирования устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона с удвоенным числом значащих цифр. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме de84r_c, при этом параметры hmin, eps, p, yx, x, h, xp, yp, a, e1, e2, e3, e4, r, r0, r1, r2, r3, yd и параметры x, y, dy в подпрограмме f должны иметь тип double.

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

ame2r_c - подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы de00r_c;
ame2d_c - подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы de00d_c;
utde20_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de00r_c;
utde21_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de00d_c;
  Кроме того, при работе подпрограмм de00r_c, de84r_c и de00d_c, de84d_c вызываются рабочие подпрограммы de00rs_c, de04rm_c, de00rd_c и de00ds_c, de04dm_c, de00dd_c соответственно.
  При работе подпрограмм de84r_c и de84d_c дополнительно вызываются также рабочие подпрограммы de84rj_c и de84dj_c соответственно.

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

 

Данная подпрограмма предназначена для интегрирования таких нелинейных систем, которые могут быть достаточно хорошо аппроксимированы линейной однородной системой с постоянными коэффициентами, т.е. системой вида

               Y '  =  AY ,
  где   A - некоторая постоянная матрица. 

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

При работе подпрограммы значения параметров m, hmin, eps, p сохраняются.

Между последовательными обращениями к подпрограмме со значениями параметра jstart = 1, 2, 3 пользователь не должен изменять содержимое массивов a, e1, e2.

При работе подпрограмм f и fj значения параметров m, x, y не должны изменяться.

При обращении к подпрограмме со значением jstart = - 1 в качестве исходных значений аргумента и решения принимаются значения параметров xp и yp соответственно, т.е. те значения, которые эти параметры получили после самого последнего обращения к подпрограмме с неотрицательным значением jstart.

После работы подпрограммы в массиве r1 содержится значение оценки абсолютной погрешности на шаге, вычисленной по правилу Рунге.

Данная подпрограмма и ее версия предназначены также для интегрирования жестких дифференциальных уравнений (1).

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

     y1'  =  (- 500 + 2 x) y1  +  (1 - x / 10 )y2  +  (1 + 1 / ( x+1 ) ) sin y1 ,
     y1(   = 49 
     y2'  =  (- 1000 - 10 x) y1  +  (1  + x / 10 ) y2  +  (2  + 1 / ( x+2 ) ) cos y2 ,
     y2(   = 50 

              1 ≤ x ≤ 6 

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

int main(void)
{
    /* Local variables */
    extern int de00r_c(U_fp, U_fp, int *, int *, float *, float *,
                       float *, float *, float *, float *, logical *,
                       float *, float *, float *, float *, float *, float *,
                       float *, float *, float *, float *, float *, float *,
                       float *, int *);
    static float hmin;
    static int ierr;
    static float a[4] /* was [2][2] */;
    extern int f_c();
    static float h__;
    static int m;
    static float p, r__[2], x, e1[4] /* was [2][2] */, e2[4] /* was [2][2] */,
                               e3[4] /* was [2][2] */, e4[4] /* was [2][2] */,
                               r0[2], r1[2], r2[2], r3[2];
    extern int fj_c();
    static int ih;
    static float yd[2], xp, yp[2], yx[2];
    static int jstart;
    static logical bul;
    static float eps;

    m = 2;
    x = 1.f;
    yx[0] = 49.f;
    yx[1] = 50.f;
    hmin = 1e-10f;
    eps = 1e-8f;
    p = 100.f;
    jstart = 0;
    h__ = .01f;
    ih = 0;
l100:
    ++ih;
    de00r_c((U_fp)f_c, (U_fp)fj_c, &m, &jstart, &hmin, &eps, &p, yx, &x, &h__,
            &bul, &xp, yp, a, e1, e2, e3, e4, r__, r0, r1, r2, r3, yd,
            &ierr);

    printf("\n %16.7e %16.7e %16.7e \n", x, yx[0], yx[1]);
    printf("\n %16.7e %16.7e %16.7e \n\n", h__, r1[0], r1[1]);
    switch (ih) {
        case 1:  goto l101;
        case 2:  goto l102;
        case 3:  goto l103;
        case 4:  goto l104;
        case 5:  goto l105;
        case 6:  goto l106;
    }
l101:
    h__ = 1.f;
    jstart = 3;
    goto l100;
l102:
    eps = .01f;
    goto l100;
l103:
    jstart = -1;
    eps = 1e-8f;
    goto l100;
l104:
    eps = .01f;
    h__ = 1e-8f;
    jstart = 3;
    goto l100;
l105:
    h__ = -.001f;
    eps = 1e-8f;
    jstart = 0;
    goto l100;
l106:
    return 0;
} /* main */

int f_c(float *x, float *y, float *dy, int *m)
{
    /* Builtin functions */
    double sin(double), cos(double);

    /* Parameter adjustments */
    --dy;
    --y;

    /* Function Body */
    dy[1] = (*x * 2.f - 500.f) * y[1] + (1.f - *x * .1f) * y[2] +
            (1.f / (*x + 1.f) + 1.f) * (float)sin(y[1]);
    dy[2] = (-1e3f - *x * 10.f) * y[1] + (*x * .1f + 1.f) * y[2] +
            (1.f / (*x + 2.f) + 2.f) * (float)cos(y[2]);
    return 0;
} /* f_c */

int fj_c(float *x, float *y, float *df, int *m)
{
    /* Builtin functions */
    double cos(double), sin(double);

#define df_ref(a_1,a_2) df[(a_2)*2 + a_1]

    /* Parameter adjustments */
    df -= 3;
    --y;

    /* Function Body */
    df_ref(1, 1) = *x * 2.f - 500.f +
                   (1.f / (*x + 1.f) + 1.f) * (float)cos(y[1]);
    df_ref(1, 2) = 1.f - *x * .1f;
    df_ref(2, 1) = -1e3f - *x * 10.f;
    df_ref(2, 2) = *x * .1f + 1.f -
                   (1.f / (*x + 2.f) + 2.f) * (float)sin(y[2]);
    return 0;
} /* fj_c */


   Результаты: 

   После первого обращения к подпрограмме
             x                                    yx(1)                             yx(2) 
      1.000039062499+00       4.895766071531+00       4.808757932455+00 
             h__                                   r1(1)                             r1(2) 
      3.906250000002-05        9.597251282538-10       -7.233271996176-09 

   После второго обращения к подпрограмме
             x                                    yx(1)                             yx(2) 
      1.000069580077+00       4.733408100787+01       4.661901490984+01 
             h__                                   r1(1)                             r1(2) 
      3.051757812500-05        1.319373647371-10        6.946114202335-02 

   После третьего обращения к подпрограмме
             x                                    yx(1)                             yx(2) 
      1.000100097655+00       4.662140584405+01       4.517257692956+01 
             h__                                   r1(1)                             r1(2) 
      6.103515625000-05       -6.208817164108-11        1.830362084864-09 

   После четвертого обращения к подпрограмме
             x                                    yx(1)                             yx(2) 
      1.000100097655+00       4.662140584411+01       4.517257692973+01 
             h__                                   r1(1)                             r1(2) 
      3.051757812500-05       -5.820766091346-11        1.851003617047-09 

   После пятого обращения к подпрограмме
             x                                    yx(1)                             yx(2) 
      1.000100107654+00       4.662117408356+01       4.517210655950+01 
             h__                                   r1(1)                             r1(2) 
      2.000000000004-08        0.000000000000 + 00      0.000000000000 + 00 

   После шестого обращения к подпрограмме
             x                                    yx(1)                             yx(2) 
      1.000068857653+00       4.735108317091+01       4.665352269082+01 
             h__                                   r1(1)                             r1(2) 
     -3.125000000004-05        3.880510727560-11       -2.250696221985-09