Текст подпрограммы и версий
de09r_c.zip , de09d_c.zip 		Тексты тестовых примеров
tde09r_c.zip , tde09d_c.zip

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

Назначение

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

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

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

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

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

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

методом Инглэнда пятого порядка точности. Решение вычисляется в одной точке XK, которая является концом интервала интегрирования.

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

1. 

Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. - М.: Наука, 1987.

2.  England R. Error estimates for Runge - Kutta type solutions to systems of ordinary differential equations. The Computer Journal. 1969 - v.12, No.2.

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

    int de09r_c (S_fp f, integer *m, real *xn, real *yn, real *xk,
                 real *hmin, real *eps, real *p, real *h, real *y, real *r1,
                 real *r2, real *r3, integer *ierr)

   Параметры
f - имя подпрограммы вычисления значений правой части системы. Первый оператор подпрограммы имеет вид:
int f (float *x, float *y, float *dy, int *m).
Здесь: x, y - значение независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в dy. В случае системы уравнений, т.е. когда m ≠ 1 параметры y и dy представляют одномерные массивы длины m (тип параметров x, y и dy: вещественный);
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 (тип: вещественный);
r1 - двумерный вещественный рабочий массив размера m*5;
r2, r3 - одномерные вещественные рабочие массивы длины m;
ierr - целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью eps. В этом случае интегрирование системы прекращается; при желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров hmin и h 
   Версии
de09d_c - вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Инглэнда с удвоенным числом значащих цифр. При этом параметры xn, yn, xk, hmin, eps, p, h, y, r1, r2, r3 и параметры x, y и dy в подпрограмме f должны иметь тип double. 
   Вызываемые подпрограммы
utde20_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de09r_c;
utde21_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de09d_c 
   Замечания по использованию
 

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

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

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

               y1' = - y1 - 5y2 ,     y1(0) = 1
               y2' =   y1  + y2 ,      y2(0) = 1 ,    0 ≤ x ≤ 2

  Точное решение системы:
               y1  =  cos2x - 3 sin2x ,
               y2  =  cos2x +  sin2x

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

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

    /* Local variables */
    extern int de09r_c(U_fp, int *, float *, float *, float *, float *,
                       float *, float *, float *, float *, float *,
                       float *, float *, int *);
    static float hmin;
    static int ierr;
    static float h__;
    static int m;
    static float p, y[2];
    extern int fengl_c();
    static float r1[10] /* was [2][5] */, r2[2], r3[2], y1, y2, xk, xn,
                 yn[2], eps;

    m = 2;
    xn = 0.f;
    yn[0] = 1.f;
    yn[1] = 1.f;
    xk = 1.f;
    hmin = 1e-12f;
    eps = 1e-10f;
    p = 100.f;
    h__ = .01f;
    de09r_c((U_fp)fengl_c, &m, &xn, yn, &xk, &hmin, &eps, &p, &h__, y, r1,
            r2, r3, &ierr);
    y1 = (float)cos(2.f) - (float)sin(2.f) * 3.f;
    y2 = (float)cos(2.f) + (float)sin(2.f);

    printf("\n %15.7e %15.7e \n", y[0], y[1]);
    printf("\n %15.7e %15.7e \n", y1, y2);
    return 0;
} /* main */

int fengl_c(float *x, float *y, float *z, int *m)
{
    /* Parameter adjustments */
    --z__;
    --y;

    /* Function Body */
    z__[1] = -y[1] - y[2] * 5.f;
    z__[2] = y[1] + y[2];
    return 0;
} /* fengl_c */
      y1 = cos(2.) - 3.*sin(2.) 
      y2 = cos(2.) + sin(2.) 


 Результаты: 

      ierr = 0 
      y(1) = -3.144039116632 + 00       y(2) = 4.931505897994-01
        y1 = -3.144039117018                  y2 = 4.931505902741-01