БЧА НИВЦ МГУ. DE95R_С (DE91R_С, DE93R_С, DE97R_С). Задача Коши для жестких систем и систем с большой константой Липшица

Текст подпрограммы и версий
de95r_c.zip de95d_c.zip de91r_c.zip de91d_c.zip de93r_c.zip de93d_c.zip de97r_c.zip de97d_c.zip

Тексты тестовых примеров
tde95r_c.zip tde95d_c.zip tde91r_c.zip tde91d_c.zip tde93r_c.zip tde93d_c.zip tde97r_c.zip tde97d_c.zip

Подпрограмма: de95r_c (версии: de91r_c, de93r_c, de97r_c)

Назначение

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

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

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

                             Y '   =   F(X,Y)  , 
       Y   =   ( y₁,...,y_M )  ,   F   =   ( f₁( X, y₁,..., y_M ),..., f_M( X, y₁,..., y_M ) )

методом типа Розенброка четвертого порядка.

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

Артемьев С.С., Демидов Г.В. A - устойчивый метод типа Розенброка четвертого порядка точности решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений/_htm_c/ Некоторые проблемы вычислительной и прикладной математики. Новосибирск: "Наука", 1975.

Современные численные методы решения обыкновенных дифференциальных уравнений/ Под ред. Дж.Холла и Дж.Уатта. М.: "Мир", 1979.

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

    int de95r_c (real *f, real *fj, real *fx, integer *m, real *xn,
                 real *yn, real *xk, real *hmin, real *eps, real *p, real *h__,
                 real *y, real *r__, integer *ierr)

   Параметры

f -	имя подпрограммы вычисления значений правой части дифференциальных уравнений. Первый оператор подпрограммы должен иметь вид: int f (float x, float y, float dy, int m). Здесь x, y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в dy; в случае системы уравнений, т.е. когда m № 1, параметры y и dy представляют одномерные массивы длиной m (тип параметров x, y и dy: вещественный);
fj -	имя подпрограммы вычисления значений элементов матрицы Якоби ¶f / ¶y правой части системы. Первый оператор подпрограммы должен иметь вид: int fj (float x, float y, float z, int m). Здесь x, y - значения независимой и зависимой переменных, соответственно. В случае системы уравнений, т.е. когда m № 1, параметр y представляет собой одномерный массив длины m, а параметр z - двумерный массив размера m*m. Значения элементов матрицы Якоби ¶f / ¶y правой части системы должны быть помещены в массив z, при этом частная производная от правой части i - го уравнения по j - ой переменной y (j) запоминается в элементе z (i,j) (тип параметров x, y и z: вещественный);
fx -	имя подпрограммы вычисления значений частных производных по x правой части системы ¶f / ¶x. Первый оператор подпрограммы должен иметь вид: int fx (float x, float y, float z1, int m). Здесь: x, y - значения независимой и зависимой переменных, соответственно. В случае системы уравнений, т.е. когда m № 1, параметры y и z1 представляют собой одномерные массивы длиной m. Значения частных производных по x правой части системы ¶f / ¶x должны быть помещены в массив z1; при этом частная производная от правой части i - го уравнения запоминается в элементе z1 (i) (тип параметров x, y и z1: вещественный);
m -	количество уравнений в системе (тип: целый);
xn, yn -	начальные значения аргумента и решения; в случае системы уравнений (т.е. m № 1) yn представляет одномерный массив длины m (тип: вещественный);
xk -	значение аргумента, при котором требуется вычислить решение задачи Коши (конец интервала интегрирования); xk может быть больше, меньше или равно xn (тип: вещественный);
hmin -	минимальное значение абсолютной величины шага, которое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
eps -	допустимая мера погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный);
p -	граница перехода, используемая при оценке меры погрешности решения (тип: вещественный);
h -	вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если xn < xk, отрицательным, если xn > xk, или без всякого учета в виде абсолютной величины;
y -	искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента xk; для системы уравнений (когда m № 1) задается одномерным массивом длины m. В случае совпадения значений параметров xn и xk значение y полагается равным начальному значению yn (тип: вещественный);
r -	одномерный рабочий массив вещественного типа длины 3mm + 11*m + 1;
ierr -	целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью eps; в этом случае интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров h и hmin.

   Версии

de97r_c -	вычисление решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом типа Розенброка четвертого порядка. Отличие подпрограммы de97r_c от подпрограммы de95r_c состоит в том, что при обращении к подпрограмме de97r_c не требуется задавать подпрограмму вычисления матрицы Якоби системы и подпрограмму вычисления частных производных правой части системы по x. Все частные производные от правой части вычисляются в подпрограмме de97r_c с помощью разностных отношений. Первый оператор подпрограммы de97r_c имеет вид: int de97r_c (real f, integer m, real xn, real yn, real xk, real hmin, real eps, real p, real h__, real y, real r__, integer ierr) Список формальных параметров подпрограммы de97r_c отличается от списка параметров подпрограммы de95r_c отсутствием параметров fj и fx. Параметры подпрограммы de97r_c имеют тот же смысл, что и одноименные параметры подпрограммы de95r_c;
de91r_c -	вычисление решения задачи Коши для автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом типа Розенброка четвертого порядка. Данная подпрограмма предназначена для интегрирования систем уравнений, в правые части которых не входит независимая переменная x, т.е. систем вида y ' = f (y). Первый оператор подпрограммы de91r_c имеет вид: int de91r_c (real f, real fj, integer m, real xn, real yn, real xk, real hmin, real eps, real p, real h__, real y, real r__, integer *ierr) Список формальных параметров подпрограммы de91r_c отличается от списка параметров de95r_c отсутствием параметра fx. Параметры подпрограммы de91r_c имеют тот же смысл, что и одноименные параметры подпрограммы de95r_c, кроме параметра r. В подпрограмме de91r_c параметр r представляет одномерный вещественный рабочий массив длины 3m² + 8m + 1.
de93r_c -	вычисление решения задачи Коши для автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Розенброка четвертого порядка. Данная подпрограмма предназначена для интегрирования систем уравнений, в правые части которых не входит независимая переменная x, т.е. систем вида y ' = f (y). Отличие подпрограммы de93r_c от подпрограммы de95r_c состоит в том, что при обращении к подпрограмме de93r_c не требуется задавать подпрограмму вычисления матрицы Якоби системы ¶f / ¶y и подпрограмму вычисления частных производных ¶f / ¶x. Все частные производные от правой части вычисляются в подпрограмме de93r_c с помощью разностных отношений. Первый оператор подпрограммы de93r_c имеет вид: int de93r_c (real f, integer m, real xn, real yn, real xk, real hmin, real eps, real p, real h__, real y, real r__, integer ierr) Список формальных параметров подпрограммы de93r_c отличается от списка параметров подпрограммы de95r_c отсутствием параметров fj и fx. Параметры подпрограммы de93r_c имеют тот же смысл, что и одноименные параметры подпрограммы de95r_c, кроме параметра r. В подпрограмме de93r_c параметр r представляет одномерный вещественный рабочий массив длины 3m² + 8m + 1.
de95d_c -	вычисление решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом типа Розенброка четвертого порядка с повышенной точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме de95r_c; при этом параметры xn, yn, xk, hmin, eps, p, h, y, r и параметры x, y, dy, z, z1 в подпрограммах f, fj и fx должны иметь тип double;
de97d_c -	вычисление решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом типа Розенброка четвертого порядка с повышенной точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме de97r_c; при этом параметры xn, yn, xk, hmin, eps, p, h, y, r и параметры x, y, dy в подпрограмме f должны иметь тип double;
de91d_c -	вычисление решения задачи Коши для автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом типа Розенброка четвертого порядка с повышенной точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме de91r_c; при этом параметры xn, yn, xk, hmin, eps, p, h, y, r и параметры x, y, dy, z в подпрограммах f и fj должны иметь тип double;
de93d_c -	вычисление решения задачи Коши для автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интергирования методом типа Розенброка четвертого порядка с повышенной точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме de93r_c; при этом параметры xn, yn, xk, hmin, eps, p, h, y, r и параметры x, y, dy в подпрограмме f должны иметь тип double.

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

de94r_c - de94d_c	выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка с обыкновенной и повышенной точностью. Вызываются при работе подпрограмм de95r_c и de95d_c соответственно;
de96r_c - de96d_c	выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка с обыкновенной и повышенной точностью. Вызываются при работе подпрограмм de97r_c и de97d_c соответственно;
de90r_c - de90d_c	выполнение одного шага численного интегрирования автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка с обыкновенной и повышенной точностью. Вызываются при работе подпрограмм de91r_c и de91d_c соответственно;
de92r_c - de92d_c	выполнение одного шага численного интегрирования автономной жесткой системы обыкновенных дифференциальных уравнений первого порядка методом типа Розенброка с обыкновенной и повышенной точностью. Вызываются при работе подпрограмм de93r_c и de93d_c соответственно;
utde20_c - utde21_c	подпрограммы выдачи диагностических сообщений. Подпрограмма utde20_c вызывается при работе подпрограмм de91r_c, de93r_c, de95r_c, de97r_c; подпрограмма utde21_c вызывается при работе подпрограмм de91d_c, de93d_c, de95d_c, de97d_c.

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

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

При работе подпрограммы и ее версий значения параметров m, xn, yn, xk, hmin, eps, p сохраняются. При работе подпрограмм f, fj и fx значения параметров m, x, y не должны изменяться.

Если после работы подпрограммы нет необходимости иметь начальное значение решения yn, то параметры yn и y при обращении к ней можно совместить. При этом следует иметь в виду, что в случае аварийного выхода из подпрограммы, т.е. со значением ierr = 65, значение параметра yn будет испорчено.

В подпрограмме de95r_c и ее версиях используется одношаговый метод типа Розенброка четвертого порядка точности. На каждом шаге h интегрирования, выполняемом из текущего узла интегрирования x_n, вычисляются четыре значения правой части системы дифференциальных уравнений. Среди этих значений есть одно, которое должно вычисляться при значении независимой переменной x, равном x_n - h. Это значение аргумента находится слева от точки x_n при h > 0 и справа от точки x_n при h < 0.

В частности, при выполнении первого шага h из начальной точки x_n (напомним, что абсолютная величина первого шаг h берется равной абсолютной величине значения, заданного параметром h при обращении к подпрограмме, а знак первого шага определяется соотношением значений параметров xn и xk) значение x_n - h не будет принадлежать интервалу интегрирования, ограниченному значениями, заданными параметрами xn и xk при обращении к подпрограмме de95r_c (или ее версиями). Это следует учитывать при составлении подпрограммы f вычисления правой части системы. Если правая часть системы не определена для x, не принадлежащих интервалу интегрирования, то попытка вычислить правую часть для указанного значения аргумента x может привести к аварийному прерыванию.

Если при решении системы диффиренциальных уравнений не задаются подпрограммы fj вычисления матрицы Якоби и fx вычисления частных производных по x правой части системы (т.е. производится обращение к подпрограммам de93r_c, de97r_c, de93d_c, de97d_c), то все частные производные по y и по x в этих подпрограммах аппроксимируются центральными разностными отношениями. Замена точных значений производных разностными аппроксимациями может привести к росту погрешности приближенного решения системы дифференциальных уравнений по сравнению с тем, что было бы, если бы все частные производные вычислялись точно с помощью подпрограмм fj и fx. Даже если будет мала погрешность аппроксимации частных производных (например, если правая часть системы является линейной функцией своих аргументов, то погрешность аппроксимации частных производных равна нулю), может оказаться значительной вычислительная погрешность, особенно, когда вычисления выполняются с одинарной точностью. При этом вычислительная погрешность может даже превосходить верхний предел погрешности приближенного решения, заданный при обращении к подпрограмме параметром eps. В этом случае целесообразно использовать не подпрограммы de93r_c, de97r_c, а их версии, выполняющие вычисления с удвоенным числом значащих цифр, т.е. подпрограммы de93d_c, de97d_c.

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

 
   Решается задача Коши

      y₁'  =  -100 y₁ 
      y₂'  =  -100 y₁  - 2 y₂  + 20 e^{- 100 x} + 2 e^{- x} cos x 
      y₃'  =  -100 y₁  + 9998 y₂ - 9990 y₃ - 10 y₄  + 20 e^{- 100 x} + 2 e^{- x} cos x 
      y₄'  =  -100 y₁  + 9988 y₂  + 20 y₃ - 10010 y₄  + 20 e^{- 100 x} + 2 e^{- x} cos x  

      0 Ј x Ј 10 ,  y₁(0) = 10 ,  y₂(0) = 11 ,  y₃(0) = 111 ,  y₄(0) = 11   

   Точное решение задачи имеет вид:

      y₁  =  10 e^{-100 x} 
      y₂  =  10 e^{- 100 x} + e^{- x} cos x + e^{- x} sin x 
      y₃  =  y₂  + 100 e^{- 10000 x} cos 10x 
      y₄  =  y₃  + 100 e^{- 10000 x} sin 10x 

   Матрица Якоби правой части системы имеет вид:

        |   -100       0            0             0      
        |   -100      -2            0             0      
        |   -100       9998     -9990     -10     
        |   -100       9988      20          -10010   

   Частные производные по  x правой части системы имеют вид:

        |                            0                 
        |   -2000 e^{- 100 x} - 2 e^{- x} (cos x + sin x)
        |   -2000 e^{- 100 x} - 2 e^{- x} (cos x + sin x)
        |   -2000 e^{- 100 x} - 2 e^{- x} (cos x + sin x)

Ниже приводятся фрагмент вызывающей программы для de95r_c, подпрограммы f, fj, fx и результаты счета, полученные после нескольких обращений к подпрограмме de95r_c.

int main(void)
{
    /* Builtin functions */
    double exp(double), cos(double), sin(double);

    /* Local variables */
    extern int de95r_c(U_fp, U_fp, U_fp, int *, float *, float *, float *,
                       float *, float *, float *, float *, float *, float *,
                       int *);
    static float hmin;
    static int ierr;
    extern int f_c();
    static float h__;
    static int m;
    static float p, r__[93], y[4], y1, y2, y3, y4;
    extern int fj_c();
    static int ih;
    extern int fx_c();
    static float xk, xn, yn[4], eps;

    m = 4;
    xn = 0.f;
    yn[0] = 10.f;
    yn[1] = 11.f;
    yn[2] = 111.f;
    yn[3] = 111.f;
    hmin = 1e-10f;
    eps = .01f;
    p = 1e3f;
    xk = 10.f;
/*  B›—€C‹EH€e TO—HOѓo PEEH€џ C€CTEM›: */
    y1 = (float)exp((float)(xk * -100.f)) * 10.f;
    y2 = y1 + (float)exp((float)(-xk)) * ((float)cos(xk) + (float)sin(xk));
    y3 = y2 + (float)exp((float)(xk * -1e4f)) * 100.f *
         (float)cos((float)(xk * 10.f));
    y4 = y3 + (float)exp((float)(xk * -1e4f)) * 100.f *
         (float)sin((float)(xk * 10.f));

    printf("\n %16.7e %16.7e ", y1, y2);
    printf("\n %16.7e %16.7e \n", y3, y4);
    ih = 0;
l10:
    ++ih;
    h__ = .01f;
    de95r_c((U_fp)f_c, (U_fp)fj_c, (U_fp)fx_c, &m, &xn, yn, &xk, &hmin, &eps, &p,
            &h__, y, r__, &ierr);

    printf("\n %7.1e \n", eps);
    printf("\n %16.7e %16.7e ", y[0], y[1]);
    printf("\n %16.7e %16.7e \n", y[2], y[3]);
    printf("\n %16.7e \n", h__);
    eps *= .01f;
    if (ih < 3) {
        goto l10;
    }
    return 0;
} /* main */

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

    /* Local variables */
    static float t1;

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

    /* Function Body */
    dy[1] = y[1] * -100.f;
    t1 = (float)exp((float)(*x * -100.f)) * 20.f +
         (float)exp((float)(-(*x))) * 2.f * (float)cos(*x);
    dy[2] = dy[1] - y[2] * 2.f + t1;
    dy[3] = dy[1] + y[2] * 9998.f - y[3] * 9990.f - y[4] * 10.f + t1;
    dy[4] = dy[1] + y[2] * 9988.f + y[3] * 20.f - y[4] * 10010.f + t1;
    return 0;
} /* f_c */

int fj_c(float *x, float *y, float *df, int *m)
{
#define df_ref(a_1,a_2) df[(a_2)*4 + a_1]

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

    /* Function Body */
    df_ref(1, 1) = -100.f;
    df_ref(1, 2) = 0.f;
    df_ref(1, 3) = 0.f;
    df_ref(1, 4) = 0.f;
    df_ref(2, 1) = -100.f;
    df_ref(2, 2) = -2.f;
    df_ref(2, 3) = 0.f;
    df_ref(2, 4) = 0.f;
    df_ref(3, 1) = -100.f;
    df_ref(3, 2) = 9998.f;
    df_ref(3, 3) = -9990.f;
    df_ref(3, 4) = -10.f;
    df_ref(4, 1) = -100.f;
    df_ref(4, 2) = 9988.f;
    df_ref(4, 3) = 20.f;
    df_ref(4, 4) = -10010.f;
    return 0;
} /* fj_c */


int fx_c(float *x, float *y, float *dx, int *m)
{
    /* Builtin functions */
    double exp(double), cos(double), sin(double);

    /* Parameter adjustments */
    --dx;
    --y;

    /* Function Body */
    dx[1] = 0.f;
    dx[2] = (float)exp((float)(*x * -100.f)) * -2e3f -
            (float)exp((float)(-(*x))) * 2.f * ((float)cos(*x) +
            (float)sin(*x));
    dx[3] = dx[2];
    dx[4] = dx[2];
    return 0;
} /* fx_c */


 Результаты: 

              y1                                   y2 
      0.000000000000+00    -6.279230870976-05 
              y3                                   y4 
     -6.279230870976-05    -6.279230870976-05 

      после первого обращения к подпрограмме -
      EPS = 1.0-02 
              y(1)                                 y(2) 
      7.146873880164-13    -6.764660892966-05 
              y(3)                                 y(4) 
     -6.764660892122-05    -6.764660892966-05 
      H = 2.560000000001+00 
   
      после второго обращения к подпрограмме -
      eps = 1.0-04 
              y(1)                                 y(2) 
      6.621227792960-20    -6.330159900469-05 
              y(3)                                 y(4) 
     -6.330159900392-05    -6.330159900414-05 
      h = 2.560000000001+00 

      после третьего обращения к подпрограмме -
      eps = 1.0-06
              y(1)                                 y(2) 
      0.000000000000+00    -6.286382905407-05 
              y(3)                                 y(4) 
     -6.286382905418-05    -6.286382905407-05 
      h = 1.280000000001+00 
   
      после четвертого обращения к подпрограмме -
      eps = 1.0-08 
              y(1)                                 y(2) 
      0.000000000000+00    -6.279451115976-05 
              y(3)                                 y(4) 
     -6.279451115976-05    -6.279451116020-05 
      h = 3.200000000002-01