Текст подпрограммы и версий
de11r_c.zip , de11d_c.zip
Тексты тестовых примеров
tde11r_c.zip , tde11d_c.zip

Подпрограмма:  de11r_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 ) ,
 методом Хойна. 

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

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

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

    int de11r_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 *rab, integer *ierr)

Параметры

f - имя подпрограммы вычисления значений правой части дифференциального уравнения. Первый оператоp подпрограммы должен иметь вид:
int f (float *x, float *y, float *dy, int *m).
Здесь: x, y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в dy. B случае системы уравнений, т.е. когда 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  
вещественные одномерные рабочие массивы длины m;
rab - вещественный двумерный рабочий массив размеpа m*2;
ierr - целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая-нибудь компонента решения не может быть вычислена с требуемой точностью eps. B этом случае интегрирование системы прекращается. При желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров hmin и h.

Версии

de11d_c - вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Хойна с повышенной точностью. При этом параметры xn, yn, xk, hmin, eps, p, h, y, yp, delty, rab и параметры x, y и dy в подпрограмме f должны иметь тип double.

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

utde10_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de11r_c.
utde11_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы de11d_c.

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

 

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

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

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

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

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

          y '1  =  y2 ,
          y '2  =  - y1 ,      3π/4 ≤ x ≤ π

          y1(3π/4)  =  √2 / 2 ,      y2(3π/4)  =  - √2 / 2 

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

int main(void)
{
    /* Builtin functions */
    double sqrt(double);

    /* Local variables */
    extern int de11r_c(U_fp, int *, 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;
    static float p, y[2], delty[2], xk, xn, yp[2],
                  rab[4] /* was [2][2] */, eps;

    y[0] = (float)sqrt(2.f) / 2;
    y[1] = -y[0];
    xk = 3.14159265359f;
    m = 2;
    xn = 2.3561944901925003f;
    h__ = .01f;
    eps = 3e-5f;
    p = 1e-6f;
    hmin = 1e-8f;
    de11r_c((U_fp)f_c, &m, &xn, y, &xk, &hmin, &eps, &p, &h__, y, yp, delty,
            rab, &ierr);

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

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

    /* Function Body */
    z__[1] = y[2];
    z__[2] = -y[1];
    return 0;
} /* f_c */


Результаты:

          y (1)  =  -0.5928428474 * 10-7
          y (2)  =  -0.9999983798

          ierr  =   0