Документ взят из кэша поисковой машины. Адрес оригинального документа : http://num-anal.srcc.msu.ru/lib_na/cat/de_htm_c/de14r_c.htm
Дата изменения: Tue Apr 28 12:37:09 2015
Дата индексирования: Sun Apr 10 02:23:55 2016
Кодировка: Windows-1251
БЧА НИВЦ МГУ. DE14R_C. Задача Коши для системы уравнений первого порядка (без контроля точности)
Текст подпрограммы и версий
de14r_c.zip , de14d_c.zip 		Тексты тестовых примеров
tde14r_c.zip , tde14d_c.zip

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

Назначение

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

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

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

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

 с начальными условиями, заданными в точке XN :
          Y(XN) = YN ,     YN = ( y10, ... , yM0 ) ,

классическим методом Рунге - Кутта четвертого порядка точности с постоянным шагом. Решение вычисляется в одной точке ХК, которая является концом интервала интегрирования.

Березин И.С., Жидков Н.П. Методы вычислений, т.2. Физматгиз, М., 1960.

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

    int de14r_c (S_fp f, integer *m, real *xn, real *yn, real *xk,
                 real *h, real *y, real *a, real *b, real *c)

   Параметры
f - имя подпрограммы вычисления значений правой части дифференциального уравнения. Первый оператоp подпрограммы должен иметь вид:
int f (float *t, float *y, float *dy, int *m).
Здесь: x, y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в dy. B случае системы уравнений, т.е. когда m ≠ 1 , параметры y и dy представляют массивы длины m (тип параметров x, y и dy: вещественный);
m - количество уравнений в системе (тип: целый);
xn, yn - начальные значения аргумента и решения; в случае системы уравнений (т.е. m ≠ 1) yn представляет одномерный массив длины m (тип: вещественный);
xk - значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования); xk может быть больше, меньше или pавно xn (тип: вещественный);
h - вещественная переменная, содержащая значение шага интегрирования. Может задаваться с учетом направления интегрирования, т.е. положительным, если xk > xn, отрицательным, если  xk < xn , или без такого учета с любым знаком " + " или " - ";
y - искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента xk. Для системы уравнений (когда m ≠ 1) задается одномерным массивом длины m. B случае совпадения значений параметров xn и xk значение y полагается равным начальному значению yn (тип: вещественный);
a, b, c - одномерные вещественные рабочие массивы длиной m. 
   Версии
de14d_c - вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования классическим методом Рунге - Кутта 4 - го порядка точности с постоянным шагом, при этом все вычисления над действительными числами выполняются с удвоенным числом значащих цифр. B этом случае параметры xn, yn, xk, h, y, a, b, c и параметры x, y и dy в подпрограмме f должны иметь тип double. 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
 

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

При работе подпрограммы и ее версии значения параметров m, xn, yn, xk сохраняются. Значение параметра h сох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 de14r_c(U_fp, int *, float *, float *, float *, float *,
                       float *, float *, float *, float *);
    static float a[4], b[4], c__[4];
    extern int f_c();
    static float h__;
    static int m, i;
    static float y[4], xk, xn, yn[4];

    m = 4;
    xn = 0.f;
    xk = 8.f;
    yn[0] = 0.f;
    yn[1] = 0.f;
    yn[2] = 0.f;
    yn[3] = 0.f;
    h__ = .01f;
    de14r_c((U_fp)f_c, &m, &xn, yn, &xk, &h__, y, a, b, c__);

    printf("\n %16.7e \n\n", xk);
    for (i = 1; i <= 4; ++i) {
         printf("\n %16.7e ", y[i-1]);
    }
    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 */


Результаты:

          y(1)  =  9.23847083543 - 03
          y(2)  =  4.82097225727 - 03
          y(3)  =  1.66711330351 + 00
          y(4)  =  1.88434937053 - 02