Текст подпрограммы и версий
de02r_c.zip , de02d_c.zip
Тексты тестовых примеров
tde02r_c.zip , tde02d_c.zip

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

Назначение

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

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

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

(1)                         Y ' (X)   =   A * Y(X) + U(X, Y)  ,

                                   Y    =   ( y1,..., yM ) , 
                                   A    =   ( ai j ) ,     i, j  =  1, ..., M ,
                           U(X, Y)   =   ( U1(X, Y), ... , UM(X, Y) ) , 

где A - постоянная числовая матрица.

Предполагается, что среди характеристических корней матрицы A имеются большие по модулю корни, а константа Липшица для функции U (X,Y), т.е. константа L из условия Липшица

| Ui ( x, y1(1),..., yM(1) ) - Ui ( x, y1(2), ... , yM(2) ) |  ≤
                                                                                      M
                                                                             ≤  L  ∑   | yj(1) - yj(2)  |  , 
                                                                                     j=1 

независящая от  i,  x,  y(1),  y(2) , невелика. Также предполагается, что нелинейный член U (x, y) является достаточно малым. По заданному значению решения YX в узле  xn  вычисляется значение этого решения в узле  xn + H . Вычисление производится по методу Лоусона.

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

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

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

(2)   Z ' (x)  =  U1(x, z)  =  exp [ - ( x - xn ) A ]  U( x, exp [( x - xn ) A ] Z(x) ) 

Данное преобразование выполняется самой подпрограммой. Характеристические корни матрицы Якоби

   ∂U1 / ∂Z   =   exp [ - ( x - xn ) A ]  (∂U / ∂y)  exp [ ( x - xn ) A ] 

являются характеристическими корнями матрицы ∂U / ∂y и в силу малости константы Липшица функции U (x, y) невелики. Поэтому система (2) не жесткая и может быть решена традиционными методами численного интегрирования. Данная подпрограмма решает ее методом Рунге - Кутта  4 - ого порядка, при этом одновременно с решением (2) производится обратное преобразование от функции Z (x) к функции  y (x) .

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

Все компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине меньше некоторой наперед заданной константы 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 de02r_c (real *fu, integer *m, integer *jstart, real *a,
                 real *hmin, real *eps, real *p, real *yx, real *x, real *h,
                 logical *bul, real *xp, real *yp, real *e1, real *e2, real *e3,
                 real *e4, real *r, real *r0, real *r1, real *r2, real *r3,
                 real *yd, integer *ierr)

Параметры

fu - подпрограмма вычисления функции u (x, y) в правой части системы. Первый оператор подпрограммы должен иметь вид:
int fu (float *x, float *y, float *u, int *m).
Здесь y и u - одномерные массивы длины m. В массив u помещается значение функции u (x, y), вычисленное при значении аргументов x и y (тип параметров x, y, u: вещественный);
m - количество уравнений в системе (тип: целый);
jstart - целый указатель режима использования подпрограммы, имеющий следующие значения:
0 - первое обращение к подпрограмме должно быть выполнено с нулевым значением jstart;
+1 - выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных и шага интегрирования, заданных параметрами x, yx и h соответственно; при этом само значение величины шага интегрирования h, на который предполагается продолжить решение системы, должно совпадать с тем значением шага, который был фактически выполнен подпрограммой при предыдущем обращении к ней. Если же требуется выполнить шаг интегрирования, отличный от реально сделанного подпрограммой во время предыдущего обращения к ней, то для этого необходимо обратиться к подпрограмме со значением jstart = 0;
+2 - то же, что и для jstart = + 1, но только с той разницей, что величина шага интегрирования, на который предполагается продолжить решение, в два раза больше того значения шага, который был фактически выполнен при предыдущем обращении к подпрограмме;
-1 - повторить последний шаг интегрирования с новыми значениями параметров h и/или hmin;
  при выходе из подпрограммы параметр jstart полагается равным 2, если выполнив заданный в h шаг, подпрограмма рекомендовала использовать для интегрирования вдвое большее его значение, и 1 в противном случае, т.е. в том случае, когда рекомендованное значение шага равно только что выполненному шагу; рекомендованное значение шага на выходе из подпрограммы запоминается в переменной h;
a - вещественный двумерный массив размера m*m, содержащий элементы матрицы A;
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 при входе в нее (т.е. предыдущий узел и решение в нем);
       e1, e2 -
       e3, e4  
вещественные двумерные рабочие массивы размера m*m;
     r, r0, r1 -
     r2, r3, yd  
вещественные одномерные рабочие массивы длины m;
ierr - целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью eps; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров h, hmin и значением jstart = - 1.

Версии

de02d_c - выполнение одного шага численного интегрирования квазилинейной системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона с удвоенным числом значащих цифр. При этом параметры hmin, eps, p, a, yx, x, h, xp, yp, e1, e2, e3, e4, r, r0, r1, r2, r3, yd и параметры x, y, dy в подпрограмме fu должны иметь тип double.

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

ame2r_c - подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы de02r_c;
ame2d_c - подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы de02d_c;
utde20_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de02r_c;
utde21_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de02d_c;
  Кроме того, при работе подпрограмм de02r_c и de02d_c вызываются рабочие подпрограммы de02rs_c, de04rq_c и de02ds_c, de04dq_c соответственно.

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

 

Данная подпрограмма предназначена для интегрирования квазилинейных систем, имеющих малый нелинейный член U (x, y).

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

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

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

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

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

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

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

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

        y1'  =  - 500 y1  +  y2  +  sin ( y1 + y2 ) ( 1 + e - x/2 ) ,
        y1(0) = 50 
        y2'  =  - 1000 y1  +  y2  +  sin ( y1 - y2 ) ( 1 + e - x/2 ) ,
        y2(0) = 50 
     
        x ≥ 0  

Приводятся подпрограмма вычисления функций

        U1 (x, y)  =  sin ( y1 + y2 ) ( 1 + e - x/2 )  ,
        U2 (x, y)  =  sin ( y1 - y2 ) ( 1 + e - x/2 )

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

int main(void)
{
    /* Local variables */
    extern int de02r_c(U_fp, int *, int *, float *, float *, float *,
                       float *, float *, float *, float *, logical *,
                       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] */, 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];
    static int ih;
    extern int fu_c();
    static float yd[2], xp, yp[2], yx[2];
    static int jstart;
    static logical bul;
    static float eps;

#define a_ref(a_1,a_2) a[(a_2)*2 + a_1 - 3]

    m = 2;
    x = 0.f;
    yx[0] = 50.f;
    yx[1] = 50.f;
    hmin = 1e-10f;
    eps = 1e-8f;
    p = 100.f;
    jstart = 0;
    h__ = .01f;
    a_ref(1, 1) = -500.f;
    a_ref(1, 2) = 1.f;
    a_ref(2, 1) = -1e3f;
    a_ref(2, 2) = 1.f;
    ih = 0;
l100:
    ++ih;
    de02r_c((U_fp)fu_c, &m, &jstart, a, &hmin, &eps, &p, yx, &x, &h__, &bul,
            &xp, yp, 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;
    goto l100;
l102:
    eps = .1f;
    goto l100;
l103:
    jstart = -1;
    eps = 1e-10f;
    goto l100;
l104:
    eps = 1e-4f;
    h__ = -1e-5f;
    jstart = -1;
    goto l100;
l105:
    h__ = -1e-4f;
    eps = 1e-8f;
    goto l100;
l106:
    return 0;
} /* main */


int fu_c(float *x, float *y, float *u, int *m)
{
    /* Builtin functions */
    double exp(double), sin(double);

    /* Local variables */
    static float t;

    /* Parameter adjustments */
    --u;
    --y;

    /* Function Body */
    t = (float)exp((float)(-(*x) / 2.f)) + 1.f;
    u[1] = t * (float)sin((float)(y[1] + y[2]));
    u[2] = t * (float)sin((float)(y[1] - y[2]));
    return 0;
} /* fu_c */


 Результаты: 
 
 После первого обращения к подпрограмме
              x                                    yx(1)                                  yx(2) 
      1.953125000001 - 05       4.951502369251 + 01       4.902915763238 + 01 
              h__                                  r1(1)                                   r1(2) 
      1.953125000001 - 05       3.523503740624 - 09        -1.435788969199 - 10 

 После второго обращения к подпрограмме
              x                                    yx(1)                                  yx(2) 
      5.004882812498 - 05       4.876668035128 + 01       4.753104564862 + 01 
              h__                                  r1(1)                                   r1(2) 
      3.051757812500 - 05       5.044663945828 - 10         1.331015179553 - 09 

 После третьего обращения к подпрограмме
              x                                    yx(1)                                  yx(2) 
      8.056640624998 - 05       4.802966483223 + 01       4.605556652445 + 01 
              h__                                  r1(1)                                   r1(2) 
      6.103515625000 - 05      -2.398155629634 - 08        -9.972912569839 - 10 

 После четвертого обращения к подпрограмме
              x                                    yx(1)                                  yx(2) 
      5.767822265623 - 05       4.858137986262 + 01       4.716006872657 + 01 
              h__                                  r1(1)                                   r1(2) 
      7.629394531250 - 06      -4.268561800317 - 11        -3.880510727560 - 12 

 После пятого обращения к подпрограмме
              x                                    yx(1)                                  yx(2) 
      4.004882812492 - 05       4.901063239173 + 01       4.801943553990 + 01 
              h__                                  r1(1)                                   r1(2) 
     -2.000000000005 - 05       8.537123600634 - 11         7.761021455121 - 12 

 После шестого обращения к подпрограмме
              x                                    yx(1)                                  yx(2) 
     -5.995117187529 - 05       4.950228966336 + 01       4.900345166435 + 01 
              h__                                  r1(1)                                   r1(2) 
     -1.000000000002 - 04      -7.190586378176 - 09        -8.769954244295 - 10