Текст подпрограммы и версий
de29r_c.zip , de29d_c.zip 		Тексты тестовых примеров
tde29r_c.zip , tde29d_c.zip

Подпрограмма:  de29r_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, которая является концом интервала интегрирования. Для всех компонент решения осуществляется контроль точности по меpе погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P (называемой границей перехода), то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.

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

Беленький В.З. Стандартная программа для интегрирования системы дифференциальных уравнений методом Адамса. Сб. "Вычислительные методы и программирование", вып. 3, Изд-во МГУ, 1965

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

    int de29r_c (real *f, integer *m, real *xn, real *yn, real *xk,
                 real *hmin, real *hmax, real *eps, real *p, real *h, real *y,
                 real *delty, real *df, real *rf, real *yp, real *ry, real *rfn,
                 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 - начальные значения аргумента и решения; в случае системы уравнений (т.е. m ≠ 1) yn представляет массив длины m (тип: вещественный);
xk - значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования); xk может быть больше, меньше или pавно xn (тип: вещественный);
hmin - минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
hmax - максимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
eps - допустимая меpа погрешности, с которой тpебуется вычислить все компоненты решения (тип: вещественный);
p - граница перехода, используемая при оценке меры погрешности решения (тип: вещественный);
h - вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если x > xn, отрицательным, если xk < xn, или без такого учета в виде абсолютной величины;
y - искомое решение задачи Коши, вычисленное подпрограммой для значения аргумента xk; для системы уравнений (когда m ≠ 1) задается одномерным массивом длины m; в случае совпадения значений параметров xk и xn значение y полагается равным начальному значению yn (тип: вещественный);
         delty -
         rf, yp  
         rfn   		одномерные вещественные рабочие массивы длины m;
df, ry - двумерные вещественные рабочие массивы размеpа m*5;
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом:
           ierr=65 -
         ierr=66   		когда какая - нибудь компонента решения не может быть вычислена с требуемой точностью eps при заданных начальном шаге h и его минимальном значении hmin; при этом ierr = 66 указывает, что требуемая точность не может быть достигнута при разгоне, а ierr = 65 - после разгона; при желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров h и hmin. 
   Версии
de29d_c - вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Адамса с повышенной точностью. При этом параметры xn, yn, xk, hmin, hmax, eps, p, h, y, delty, df, rf, yp, ry, rfn и параметры x, y и dy в подпрограмме f должны иметь тип double. 
   Вызываемые подпрограммы
de28r_c - выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом Адамса с контролем точности.
de28d_c - выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом Адамса с повышенной точностью.
      utde12_c -       utde13_c   подпрограммы выдачи диагностических сообщений.
  Подпрограммы de28r_c и utde12_c вызываются при работе подпрограммы de29r_c, а подпрограммы de28d_c и utde13_c - при работе подпрограммы de29d_c. 
   Замечания по использованию
 

B общем случае заданая точность не гарантируется.

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

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

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

           y1'  =   y2
           y2'  =  -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 de29r_c(U_fp, int *, float *, float *, float *, float *,
                       float *, float *, float *, float *, float *, float *,
                       float *, float *, float *, float *, float *, int *);
    static float hmin, hmax;
    static int ierr;
    extern int f_c();
    static float h__;
    static int m;
    static float p, y[2], delty[2], df[10];
    static int ih;
    static float rf[2], xk, xn, yn[2], yp[2], ry[10], rfn[2], eps;

    m = 2;
    xn = 2.3561947499999998f;
    yn[0] = (float)sqrt(2.f) / 2.f;
    yn[1] = -yn[0];
    xk = 3.141593f;
    hmin = 1e-18f;
    eps = 1e-5f;
    h__ = .01f;
    ih = 0;
    p = 100.f;
    hmax = .1f;
    de29r_c((U_fp)f_c, &m, &xn, yn, &xk, &hmin, &hmax, &eps, &p, &h__, y,
            delty, df, rf, yp, ry, rfn, &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 *t, float *y, float *z__, int *m)
{
    /* Parameter adjustments */
    --z__;
    --y;

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


Результаты:

          ierr  =  0
          
          y(1)  =   3.010855175-07
          y(2)  =  -1.000001503+00