Текст подпрограммы и версий
de13r_c.zip , de13d_c.zip 		Тексты тестовых примеров
tde13r_c.zip , tde13d_c.zip

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

Назначение

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

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

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

          Y ' = F (X, Y) ,
          Y = ( y1, ... , yM ) ,
          F = ( f1 (X, y1, ... , yM), ... , fM (X, y1, ... , yM) )
 с начальными условиями, заданными в точке XN :

          Y(XN) = YN ,     YN = ( y10, ... , yM0 ) ,

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

О.Б.Арушанян, Стандартная программа решения задачи Коши для систем обыкновенных дифференциальных уравнений методом Рунге - Кутта, вып. 31, под общей редакцией В.В.Воеводина, НИВЦ МГУ, 1968.

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

    int de13r_c (S_fp f, integer *m, real *xn, real *yn, real *xk,
                 real *hmin, real *eps, real *p, real *h, real *y, real *yp,
                 real *delty, real *yr, real *dy, integer *ierr)

   Параметры
f - имя подпрограммы вычисления значений правой части дифференциального уравнения. Первый оператоp подпрограммы должен иметь вид:
int f (float *t, float *y, float *dy, int *m).
Здесь: x, y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в dy. В случае системы уравнений, т.е. когда m ≠ 1 , параметры y и dy представляют массивы длины m (тип параметров x, y и dy: вещественный);
m - количество уравнений в системе (тип: целый);
xn, yn - начальные значения аргумента и решения. B случае системы уравнений (т.е. m ≠ 1) yn представляет одномерный массив длины m (тип: вещественный);
xk - значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования). xk может быть больше, меньше или pавно xn (тип: вещественный);
hmin - минимальное значение абсолютной величины шага, который разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
eps - допустимая меpа погрешности, с которой тpебуется вычислить все компоненты решения (тип: вещественный);
p - граница перехода, используемая при оценке погрешности решения (тип: вещественный);
h - вещественная переменная, содержащая начальное значение шага интегрирования. Может задаваться с учетом направления интегрирования, т.е. положительным, если xk > xn, отрицательным, если xk < xn, или без такого учета в виде абсолютной величины;
y - искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента xk. Для системы уравнений (когда m ≠ 1) задается одномерным массивом длины m. B случае совпадения значений параметров xn и xk значение y полагается равным начальному значению yn (тип: вещественный);
            yp -
       delty  
       yr, dy   		вещественные одномерные рабочие массивы длины m;
ierr - целая переменная, значение которой в pезультате работы подпрограммы полагается pавным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью eps. B этом случае интегрирование системы прекращается. При желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметpов hmin и h. 
   Версии
de13d_c - вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования классическим методом Рунге - Кутта четвертого порядка с повышенной точностью. При этом параметры xn, yn, xk, hmin, eps, p, h, y, yp, delty, yr, dy и параметры x, y и dy в подпрограмме f должны иметь тип double. 
   Вызываемые подпрограммы
utde10_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de13r_c.
utde11_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de13d_c. 
   Замечания по использованию
 

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

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

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

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

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

          y1'  =  0.2 ( y4 - y1 )
          y2'  =  y1 + 2 ( y2 - y2 y3 )
          y3'  =  y4 - ( y3 - y2 y3 )
          y4'  =  10 y1 - ( 61 - 0.13 x ) y4 + 0.13 x ,     0 ≤ x ≤ 8 ,
     
          y1 (0)  =  y2 (0)  =  y3 (0)  =  y4 (0)  =  0

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

int main(void)
{
    /* Local variables */
    extern int de13r_c(U_fp, int *, float *, float *, float *, 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, i;
    static float p, delty[4], xf, yf[4], xi, yi[4], yp[4], rab[4],
                 err, rab1[4];

    m = 4;
    yi[0] = 0.f;
    yi[1] = 0.f;
    yi[2] = 0.f;
    yi[3] = 0.f;
    xi = 0.f;
    xf = 8.f;
    h__ = .01f;
    hmin = 1e-16f;
    err = 1e-4f;
    p = 1e-8f;
    de13r_c((U_fp)f_c, &m, &xi, yi, &xf, &hmin, &err, &p, &h__, yf, yp, delty,
            rab, rab1, &ierr);

    printf("\n %16.7e \n\n", xf);
    for (i = 1; i <= 4; ++i) {
         printf("\n %16.7e ", yf[i-1]);
    }
    printf("\n\n %5i \n", ierr);
    return 0;
} /* main */

int f_c(float *t, float *y, float *z, int *m)
{
    static float r23, ct;

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

    /* Function Body */
    r23 = y[2] * y[3];
    ct = *t * .13f;
    z__[1] = (y[4] - y[1]) * .2f;
    z__[2] = y[1] + (y[2] - r23) * 2.f;
    z__[3] = y[4] - (y[3] - r23);
    z__[4] = y[1] * 10.f - (61.f - ct) * y[4] + ct;
    return 0;
} /* f_c */


Результаты:

          yf(1)  =  0.00923847083381
          yf(2)  =  0.00482097150315
          yf(3)  =  1.66711319966
          yf(4)  =  0.0188434937933

          ierr  =  0