Текст подпрограммы и версий
de59r_c.zip  de59d_c.zip  de60r_c.zip  de60d_c.zip  de57r_c.zip  de57d_c.zip  de58r_c.zip  de58d_c.zip
Тексты тестовых примеров
tde59r_c.zip  tde59d_c.zip  tde60r_c.zip  tde60d_c.zip  tde57r_c.zip  tde57d_c.zip  tde58r_c.zip  tde58d_c.zip

Подпрограмма:  de59r_c (версии: de60r_c, de57r_c, de58r_c)

Назначение

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

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

Решается двухточечная краевая задача для линейного самосопряженного уравнения второго порядка

(1)      d (k(x) dy/dx) / dx - g(x) y = - r(x)   ,      xN < x < xK
с разрывными коэффициентами
(2)      k(x) ,  g(x) ,  r(x)   Q( 2 ) [xN, xK]

( Q ( 2 ) [xN, xK] - класс функций, кусочно - непрерывных на  [xN, xK] вместе со вторыми производными ), удовлетворяющими неpавенствам

(3)      k(x) ≥ c > 0   ,   g(x) ≥ 0   .

Hа концах отрезка  [xN, xK] задаются произвольные комбинации следующих краевых условий

при  x = xN
(4)               1)   y( xN ) = GN   ,
(5)               2)   k dy/dx - ( σN y - μN ) = 0   ,   σN ≥ 0   ,

при  x = xK
(6)               1)   y( xK ) = GK   ,
(7)               2)   - k dy/dx - ( σK y - μK ) = 0   ,   σK ≥ 0   ;

при  σN > 0 условие (5) (а при  σK > 0 условие (7)) является условием третьего рода, при  σN = 0 (или при  σK = 0) - условием второго рода.

Решение вычисляется на равномерной сетке узлов, заданной на  [xN, xK], с помощью однородной консервативной разностной схемы, имеющей на указанном классе коэффициентов (2) второй порядок точности, при этом предполагается, что точки разрыва коэффициентов  k (x), g (x) и  r (x) уравнения (1) совпадают с узлами этой сетки.

Самарский A.A., Теория разностных схем, "Hаука", M., 1977.

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

    int de59r_c (real *ck, real *cg, real *cr, integer *nx,
                 real *xn, real *xk, logical *gn1, real *gn, logical *gk1, real *gk,
                 real *y, doublereal *alp, doublereal *bet, integer *ierr)

Параметры

ck - одномерный вещественный массив длины  nx - 1, содержащий на входе в подпрограмму значения коэффициента  k (x) уравнения, а именно:  i - й элемент  ck (i) массива содержит значение функции  k (x) в точке
x = xn + (i - 1/2) * (xk - xn)/(nx - 1); функция  k (x) должна принимать положительные значения;
cg, cr - одномерные вещественные массивы длины  2*nx - 2, содержащие на входе в подпрограмму значения коэффициентов  g (x) и  r (x), соответственно, а именно:
cg (1) содержит предел справа функции  g (x) в начале  xn отрезка интегрирования,
(2*i) - й и  (2*i + 1) - й элементы  cg (2*i),  cg (2*i + 1) массива  cg содержат, соответственно, пределы слева и справа функции  g (x) в узле сетки
x = xn + i*(xk - xn)/(nx - 1)  (i = 1, 2, ..., nx - 2),
(2*nx - 2) - й элемент  cg (2*nx - 2) массива  cg содержит предел слева функции  g (x) в конце  xk отрезка интегрирования,
элементы массива  cr содержат аналогичные значения, вычисленные для функции  r (x);
nx - число узлов равномерной сетки, включая концы отрезка интегрирования, в которых требуется вычислить решение краевой задачи;  nx > 1 (тип: целый);
xn, xk - концы отрезка интегрирования, в которых задаются граничные условия (4) или (5) и (6) или (7), соответственно;  xk > xn (тип: вещественный);
gn1 - логическая переменная (или логическая константа), имеющая на входе в подпрограмму значение true, если в начале  xn отрезка интегрирования задано граничное условие первого рода (4), и false в противном случае, т.е. когда в  xn задано условие 2 - го или 3 - его рода (5);
gn - если  gn1 = true, то  gn представляет значение решения в начале  xn отрезка интегрирования;
если  gn1 = false, то  gn является именем одномерного массива длины 2, первый элемент  gn (1) которого содержит значение  σn, а второй элемент  gn (2) - значение  mn, входящие в граничное условие (5) (тип: вещественный);
gk1 - логическая переменная (или логическая константа), имеющая на входе в подпрограмму значение true, если в конце  xk отрезка интегрирования задано граничное условие первого рода (6), и false в противном случае, т.е. когда в  xk задано условие 2 - го (σk = 0) или 3 - го (σk > 0) рода (7);
gk - если  gk1 = true, то  gk представляет значение решения в конце отрезка интегрирования;
если  gk1 = false, то  gk является именем одномерного массива длины 2, первый элемент  gk (1) которого содержит значение  σk, а второй элемент  gk (2) - значение  μk, входящие в граничное условие (7) (тип: вещественный);
y - одномерный вещественный массив длины  nx, содержащий на выходе из подпрограммы решение краевой задачи, вычисленное на равномерной сетке; при этой  i - й элемент  y (i) этого массива содержит приближенное значение решения в узле
xn + (i - 1)*(xk - xn)/(nx - 1);
         alp -
         bet  
одномерные вещественные рабочие массивы двойной точности длины  nx - 1;
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом:
ierr=65 - если в граничном условии (5) коэффициент  σn < 0;
ierr=66 - если в граничном условии (7) коэффициент  σk < 0;
ierr=67 - если какое - нибудь значение коэффициента  k (x), заданного в массиве  ck, меньше или pавно 0;
ierr=68 - если какое - нибудь значение коэффициента  g (x), заданного в массиве  cg, меньше 0;
ierr=69 - если значение  nx ≤ 1;
ierr=70 - если  xk ≤ xn;
  во всех указанных случаях интегрирование уравнения прекращается.

Версии

de60r_c - вычисление решения двухточечной краевой задачи для обыкновенного линейного самосопряженного дифференциального уравнения второго порядка с разрывными коэффициентами на неравномерной сетке.
de57r_c - вычисление решения двухточечной краевой задачи для обыкновенного линейного самосопряженного дифференциального уравнения второго порядка с непрерывными коэффициентами на равномерной сетке.
de58r_c - вычисление решения двухточечной краевой задачи для обыкновенного линейного самосопряженного дифференциального уравнения второго порядка с непрерывными коэффициентами на неравномерной сетке.
de57d_c - то же, что и de57r_c, при этом все арифметические операции с вещественными числами выполняются с удвоенным числом значащих цифр.
de58d_c - то же, что и de58r_c, при этом все арифметические операции с вещественными числами выполняются с удвоенным числом значащих цифр.
de59d_c - то же, что и de59r_c, при этом все арифметические операции с вещественными числами выполняются с удвоенным числом значащих цифр.
de60d_c - то же, что и de60r_c, при этом все арифметические операции с вещественными числами выполняются с удвоенным числом значащих цифр.
 

Первые операторы этих подпрограмм имеют следующий вид:

int de60r_c (real *ck, real *cg, real *cr, integer *nx, real *x, logical *gn1, real *gn, logical *gk1, real *gk, real *y, doublereal *alp, doublereal *bet, integer *ierr)

int de57r_c (R_fp fk, R_fp fg, R_fp fr, integer *nx, real *xn, real *xk, logical *gn1, real *gn, logical *gk1, real *gk, real *y, doublereal *alp, doublereal *bet, integer *ierr)

int de58r_c(R_fp fk, R_fp fg, R_fp fr, integer *nx, real *x, logical *gn1, real *gn, logical *gk1, real *gk, real *y, doublereal *alp, doublereal *bet, integer *ierr)

Параметр  x в подпрограммах de60r_c и de58r_c - имя одномерного вещественного массива длины  nx, в последовательных элементах которого, начиная с первого, размещаются в возрастающем порядке узлы неравномерной сетки, причем начало и конец отрезка интегрирования, в которых заданы граничные условия (4) - (7), хранятся, соответственно, в первом и последнем элементах  x (1), x (nx) этого массива.

Параметры  fk, fg, fr в подпрограммах de57r_c и de58r_c - имена вещественных подпрограмм - функций

float fk (float *x),

float fg (float *x),

float fr (float *x),

вычисляющих значения коэффициентов уравнения  k (x), g (x) и  r (x), соответственно, при вещественном значении аргумента  x = x.

Остальные параметры указанных версий имеют такой же смысл, что и одноименные параметры подпрограммы de59r_c.

При этом значение параметра  ierr = 70 в подпрограммах de60r_c и de58r_c указывает на то, что узлы сетки расположены в массиве  x не в возрастающем порядке, значения  ierr = 67, 68 в п/п de58r_c и de57r_c - что значения коэффициентов  k (x) и  g (x) уравнения, вычисленные с помощью подпрограмм - функций  fk и  fg, не удовлетворяют требуемое неpавенство (3).

B подпрограмме de60r_c в массиве  ck размещаются значения коэффициента  k (x) уравнения, вычисленные в "потоковых" точках, т.е. в середине между узлами сетки, а в массивах  cg и  cr - пределы слева и справа функций  g (x) и  r (x) в узлах неравномерной сетки.

Kpоме того, каждая из подпрограмм de57r_c, de58r_c, de59r_c и de60r_c имеет версию de57d_c, de58d_c, de59d_c и de60d_c, соответственно, с теми же самыми эффектом и списком параметpов, выполняющую все операции над действительными числами с удвоенным числом значащих цифр.

При этом параметры  fk, fg, fr (в подпрограммах de57d_c, de58d_c), являются именами подпрограмм - функций
double fk (double *x)
double fg (double *x)
double fr (double *x).

Формальные параметры этих подпрограмм :
         ck, cg -
         cr  
(в подпрограммах de59d_c, de60d_c),
xn, xk - (в подпрограммах de57d_c, de59d_c),
x - (в подпрограммах de58d_c, de60d_c),
       gn, gk -
     y, alp  
         bet  
(во всех подпрограммах)
  должны иметь тип double.
Тип остальных параметров не изменяется.

В версиях подпрограмм de57r_c, de58r_c, de59r_c и de60r_c, формальные параметры alp и bet имеют тип double.

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

utde14_c - подпрограмма выдачи диагностических сообщений при работе подпрограмм de57r_c, de58r_c, de59r_c и de60r_c.
utde15_c - подпрограмма выдачи диагностических сообщений при работе подпрограмм de57d_c, de58d_c, de59d_c и de60d_c.

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

 

Значение параметра  nx, задаваемое при обращении к подпрограмме, должно быть достаточно большим, чтобы получить хорошее приближение численного решения к точному решению задачи.

При интегрировании дифференциального уравнения (1) с непрерывными коэффициентами с помощью подпрограмм de57r_c, de58r_c, de57d_c, de58d_c значения коэффициента  k (x) вычисляются в "потоковых" точках, т.е. в середине между узлами, а значения коэффициентов  g (x) и  r (x) - в узлах сетки последовательно слева направо (начало и конец  xn, xk отрезка интегрирования также считаются узлами сетки).

Поэтому эти продпрограммы можно использовать также и в случае, когда какой - нибудь из коэффициентов  k (x), g (x) и  r (x) (или, быть может, все) задан в виде массива значений в указанных точках.

Эти значения должны быть расположены в последовательных элементах массива, начиная с первого.

При этом: если в  xn задано граничное условие 2 - го или 3 - го рода, то значения коэффициентов  g (x) и  r (x) должны быть заданы, начиная с первого узла  xn сетки, а если 1 - го рода, то начиная со второго узла  xn + (xk - xn)/(nx - 1); если в конце промежутка  xk задано условие 1 - го рода, то значения коэффициентов  g (x) и  r (x) используются только до  (nx - 1) - го узла включительно.

B этом случае каждое обращение к соответствующей подпрограмме - функции вычисления коэффициента сводится к выборке значения очередного элемента массива.

Значения параметров  ck, cg, cr, nx, xk, x, gn1, gn, gk1, gk сохраняются.

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

       d (k(x) dy/dx) / dx = 0 ,  y(0) = 1 ,  y(1) = 0 ,
 
         k(x) = 1 ,    x ≤ 0.55 ,
         k(x) = 10 ,  x > 0.55 ;

Точное решение:

       y = 1 - x/0.595 ,  x ≤ 0.55 ,  (1 - x) / 5.95 ,  x > 0.55 .

int main(void)
{
    /* Local variables */
    extern int de59r_c(float *, float *, float *, int *, float *, float *,
                       logical *, float *, logical *, float *, float *,
                       double *, double *, int *);
    static int ierr, i__;
    static float y[21], cg[40], ck[20], gk, cr[40], gn, xk;
    static int nx;
    static float xn;
    static logical gk1, gn1;
    static double bet[20], alp[20];

    nx = 21;
    xn = 0.f;
    xk = 1.f;
    gn1 = TRUE_;
    gn = 1.f;
    gk1 = TRUE_;
    gk = 0.f;
    for (i__ = 1; i__ <= 40; ++i__) {
        cg[i__ - 1] = 0.f;
        cr[i__ - 1] = 0.f;
/* L200: */
    }
    for (i__ = 1; i__ <= 10; ++i__) {
/* L5: */
        ck[i__ - 1] = 1.f;
    }
    for (i__ = 11; i__ <= 20; ++i__) {
/* L15: */
        ck[i__ - 1] = 10.f;
    }
    de59r_c(ck, cg, cr, &nx, &xn, &xk, &gn1, &gn, &gk1, &gk, y, alp, bet,
            &ierr);

    printf("\n %5i \n", ierr);
    printf("\n %16.7e \n", xn);
    printf("\n %16.7e \n", xk);
    printf("\n %5i \n", nx);
    for (i__ = 0; i__ <= 18; i__ += 3) {
         printf("\n %16.7e %16.7e %16.7e \n",
                y[i__], y[i__ + 1], y[i__ + 2]);
    }
    return 0;
} /* main */


Результаты:

    ierr  =  0

    y(1)    =  1.00000000000 + 00       y(11)  =  9.09090909100 - 02
    y(2)    =  9.09090909089 - 01        y(12)  =  8.18181818191 - 02
    y(3)    =  8.18181818179 - 01        y(13)  =  7.27272727281 - 02
    y(4)    =  7.27272727271 - 01        y(14)  =  6.36363636372 - 02
    y(5)    =  6.36363636362 - 01        y(15)  =  5.45454545463 - 02
    y(6)    =  5.45454545453 - 01        y(16)  =  4.54545454552 - 02
    y(7)    =  4.54545454545 - 01        y(17)  =  3.63636363643 - 02
    y(8)    =  3.63636363636 - 01        y(18)  =  2.72727272732 - 02
    y(9)    =  2.72727272727 - 01        y(19)  =  1.81818181822 - 02
    y(10)  =  1.81818181818 - 01         y20)  =  9.09090909109 - 03

Абсолютная погрешность во всех узлах сетки pавна  0.00000000000