БЧА НИВЦ МГУ. ID10R_C. Численное дифференцирование на основе разностных формул

Документ взят из кэша поисковой машины. Адрес оригинального документа : http://num-anal.srcc.msu.su/lib_na/cat/i_htm_c/id10r_c.htm
Дата изменения: Thu May 21 14:58:18 2015
Дата индексирования: Sun Apr 10 01:01:24 2016
Кодировка: Windows-1251

Текст подпрограммы и версий
id10r_c.zip , id10d_c.zip

Тексты тестовых примеров
tid10r_c.zip , tid10d_c.zip

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

Назначение

Вычисление первой или второй производной вещественной таблично - заданной функции на равномерной сетке с использованием центральных разностей.

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

id10r_c вычисляет либо первые, либо вторые производные таблично - заданной функции на равномерной сетке с шагом h с использованием центральных разностей по формулам Стирлинга.

И.С.Березин, Н.П.Жидков, Методы вычислений, т.1, M., 1962.

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

    int id10r_c (real *h__, real *y, integer *mu, integer *nl,
            integer *ns, integer *nu, integer *no, integer *lo, integer *nth, 
            real *t, real *a, real *e, integer *la, integer *ierr)

   Параметры

h -	заданный шаг равномерной сетки (тип: вещественный);
y -	вещественный двумерный массив, в котоpом задается таблица значений функции и центральных разностей (см. описание параметра lo); размерность массива y должна быть при lo = 0 pавна mu на (no + 1), при lo = 1 - mu на [(no + 3)/2], при lo = 2 - mu на [(no + 4)/2];
mu -	заданное число стpок в описании размерности двумерного массива y в программе, вызывающей данную программу (тип: целый);
nl, ns - nu	определяют расположение в массиве y заданных значений функций, а именно, значения функции y₀, y₁, y₂, ..., y_n должны быть последовательно расположены в следующих компонентах массива y: y (nl, 1), y (nl + ns, 1), y (nl + 2*ns, 1), ..., y (nu, 1); расположение в массиве y центральных разностей определяется ниже при описании параметра lo (тип: целый);
no -	заданный максимальный порядок разностей, которые используются в конечно - разностных аппроксимациях производных (тип: целый);
lo -	определяет порядок расположения заданных центральных разностей в массиве y (тип: целый); при этом если:

lo = 0 -	то δ ^2j - 1y_i - 1/2 должна быть расположена в элементе массива y (nl + i* ns - [ns/2], 2j), а δ ^2jy_i - в элементе массива y (nl + ins, 2*j + 1);
lo = 1 -	то δ ^2j - 1y_i - 1/2 должна быть расположена в элементе массива y (nl + i* ns - [ns/2], j + 1), а δ ^2jy_i - в элементе массива y (nl + i*ns, j + 1);
lo = 2 -	то δ ^2j - 1y_i - 1/2 должна быть расположена в элементе массива y(nl + i* ns - [ns/2], j + 1), а δ ^2jy_i - в элементе массива y (nl + i*ns, j + 2);

nht -	полагается равным 1, если надо вычислить первую производную, или 2, если - вторую;
t -	заданная точность вычисления производных (тип: вещественный);
a -	вещественный вектоp длины nu, в котоpом помещаются вычисленные значения первой или второй производной последовательно в следующем порядке a (nl), a (nl + ns), a (nl + 2*ns), ... ,a (nu);
e -	вещественный вектоp длины nu, в котоpом помещаются последние члены ряда, использованные в конечно - разностной аппроксимации, в следующем порядке e (nl), e (nl + ns), e (nl + 2*ns),...,e (nu);
la -	целая переменная, значение которой полагается равным наивысшему порядку разностей, фактически использованных в вычислениях;
ierr -	целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом:

ierr= 1 -	когда для достижения заданной точности требуются разности более высокого порядка, чем они заданы при обращении к подпрограмме; производные вычисляются, хотя требуемая точность не достигнута;
ierr= 2 -	когда разности порядка no pасходятся; производные вычисляются, хотя требуемая точность не достигнута;
ierr= 3 -	когда разности порядка la (la < no) расходятся; производные вычисляются, хотя требуемая точность не достигнута;
ierr=65 -	когда значение nth ≠ 1 или ≠ 2;
ierr=66 -	когда неправильно задан порядок разностей; порядок разностей должен быть по крайней меpе, не ниже порядка вычисляемой производной;
ierr=67 -	когда не заданы разности нечетного порядка и требуется вычислить первую производную;
ierr=68 -	когда не заданы разности четного порядка и требуется вычислить вторую производную.

   Версии

id10d_c -

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

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

utid10_c -	подпрограмма выдачи диагностических сообщений при работе подпрограммы id10r_c.
utid11_c -	подпрограмма выдачи диагностических сообщений при работе подпрограммы id10d_c.

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

	1.	Для id10d_c параметры h, y, t, a и e должны иметь тип double.
	2.	Значение параметра mu должно быть больше или pавно значению параметра nu.
	3.	Формула вычисления первой производной использует разности только нечетных порядков. Поэтому можно задать lo = 2 и ns = 1.
	4.	Формула вычисления второй производной использует разности только четных порядков. Поэтому можно задать lo = 1 и ns = 1.
	5.	Заданная точность считается достигнутой, если следующий член ряда конечно - разностной аппроксимации меньше или pавен t.
	6.	Конечно - разностные аппроксимации не могут быть использованы в конце таблицы задания функции. Поэтому следует увеличить таблицу значений функции и ее разностей перед обращением к подпрограмме для вычисления производных вблизи конца таблицы.
	7.	Для вычисления соответствующего числа центральных разностей рекомендуется использовать подпрограмму i i10r_c.

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

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

    /* Local variables */
    extern int id10r_c(float *, float *, int *, int *, int *, int *, int *,
                       int *, int *, float *, float *, float *, int *, int *),
               ii10r_c(int *, int *, int *, int *, int *, int *, float *);
    static int ierr;
    static float a[10], e[10], h__;
    static int i__, j;
    static float t, x, y[50] /* was [10][5] */;
    static int la, nl, lo, no, ns, mu, nu;
    static float xn;
    static int nth;

#define y_ref(a_1,a_2) y[(a_2)*10 + a_1 - 11]

    h__ = .25f;
    xn = 2.f;
    x = xn - h__;
    nth = 1;
    t = 1e-4f;
    for (i__ = 1; i__ <= 10; ++i__) {
        for (j = 1; j <= 5; ++j) {
/* l3: */
            y_ref(i__, j) = 0.f;
        }
    }
    for (i__ = 1; i__ <= 10; ++i__) {
        x += h__;
        y_ref(i__, 1) = (float)exp((float)(-x));
/* l1: */
    }
    mu = 10;
    nl = 1;
    ns = 1;
    nu = 10;
    lo = 2;
    no = 6;
    ii10r_c(&mu, &nl, &ns, &nu, &no, &lo, y);

    for (i__ = 0; i__ <= 48; i__+= 2) {
        printf("\n  %16.7e %16.7e \n",y[i__], y[i__+1]);
    }
    id10r_c(&h__, y, &mu, &nl, &ns, &nu, &no, &lo, &nth, &t, a, e, &la, &ierr);

    for (i__ = 1; i__ <= 10; ++i__) {
        printf("\n %16.7e %16.7e \n",a[i__-1], e[i__-1]);
    }
/* l2: */
    return 0;
} /* main */


Результат:

         x                   a                              e
       2.         -6.207939493-02      0.000000000+00
       2.25     -1.081979987-01      2.160001819-05
       2.50     -8.205740855-02      1.682211116-05
       2.75     -6.392797405-02     -8.498944868-06
       3.00     -4.978715625-02     -6.618985062-06
       3.25     -3.877427628-02     -5.154870644-06
       3.50     -3.019743673-02     -4.014617251-06
       3.75     -2.343970501-02      7.495577201-05
       4.00     -1.832136756-02      7.671346584-05
       4.25     -8.102809959-03      0.000000000+00