Текст подпрограммы и версий
de03r_c.zip , de03d_c.zip 		Тексты тестовых примеров
tde03r_c.zip , tde03d_c.zip

Подпрограмма:  de03r_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) ) ,

с начальными условиями, заданными в точке XN: 

                         Y(XN)  =  YN ,   YN  =  ( y1, ... , yM ) .

где 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) является достаточно малым. Решение вычисляется в одной точке ХК, которая является концом интервала интегрирования. Для интегрирования системы применяется метод Лоусона.

Метод Лоусона является одношаговым и заключается в следующем. Допустим, что искомое решение системы (1) уже вычислено в некоторой точке  x = xn  интервала интегрирования, т.е. известно  yn ≈ y(xn). Для отыскания решения Y(xn + 1) = Y(xn + H) в следующем узле  xn + 1 = xn + H выполняются такие действия. Исходная система уравнений с помощью замены искомой функции Y (x) на  xn ≤ x ≤ 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) .

Все компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине меньше некоторой наперед заданной константы 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 de03r_c (real *fu, integer *m, real *xn, real *yn,
                 real *xk, real *a, real *hmin, real *eps, real *p, real *h,
                 real *y, real *r, 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 - количество уравнений в системе (тип: целый);
xn, yn - начальные значения аргумента и решения; в случае системы уравнений (т.е. m ≠ 1) yn представляет одномерный массив длины m (тип: вещественный);
xk - значение аргумента, при котором требуется вычислить решение задачи Коши (конец интервала интегрирования); xk может быть больше, меньше или равно xn (тип: вещественный);
a - вещественный двумерный массив размера m*m, содержащий элементы матрицы A;
hmin - минимальное значение абсолютной величины шага, которое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
eps - допустимая мера погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный);
p - граница перехода, используемая при оценке меры погрешности решения (тип: вещественный);
h - вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если xn < xk, отрицательным, если xn > xk, или без всякого учета в виде абсолютной величины;
y - искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента xk; для системы уравнений (когда m ≠ 1) задается одномерным массивом длины m. В случае совпадения значений параметров xn и xk значение y полагается равным начальному значению yn (тип: вещественный);
r - вещественный одномерный рабочий массив длины 4*m*m+7*m+1;
ierr - целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью eps. В этом случае интегрирование системы прекращается; при желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров hmin и h 
   Версии
de03d_c - вычисление решения задачи Коши для квазилинейной системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица в конце интервала интегрирования методом Лоусона с удвоенным числом значащих цифр. При этом параметры xn, yn, xk, a, hmin, eps, p, h, y, r и параметры x, y, r в подпрограмме fu должны иметь тип double. 
   Вызываемые подпрограммы
      de02r_c -
      de02d_c   		выполнение одного шага численного интегрирования квазилинейной системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона;
      utde20_c -
      utde21_c   		подпрограмма выдачи диагностических сообщений; Подпрограмма de02r_c, utde20_c вызывается при работе подпрограммы de03r_c, а подпрограммы de02d_c, utde21_c - при работе de03d_c. 
   Замечания по использованию
 

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

В общем случае заданная точность не гарантируется. При работе подпрограммы значения параметров m, a, xn, yn, xk, hmin, eps, p сохраняются.

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

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

Подпрограммы de03r_c и de03d_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 

        0 ≤ x ≤ 5

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

        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 de03r_c(U_fp, int *, 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__[31], y[2];
    extern int fu_c();
    static float xk, xn, yn[2], eps;

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

    m = 2;
    xn = 0.f;
    yn[0] = 50.f;
    yn[1] = 50.f;
    hmin = 1e-10f;
    eps = 1e-4f;
    p = 100.f;
    xk = 5.f;
    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;
    de03r_c((U_fp)fu_c, &m, &xn, yn, &xk, a, &hmin, &eps, &p, &h__, y, r__,
            &ierr);

    printf("\n %16.7e %16.7e \n", y[0], y[1]);
    printf("\n %16.7e \n", h__);
    printf("\n %5i \n", ierr);
    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 */


 Результаты: 

             y(1)                                   y(2)                               h 
     -3.508894600484-04      -8.340807464992-02       1.250000000001-03 

      ierr = 0