|
Текст подпрограммы и версий id10r_c.zip , id10d_c.zip |
Тексты тестовых примеров tid10r_c.zip , tid10d_c.zip |
Вычисление первой или второй производной вещественной таблично - заданной функции на равномерной сетке с использованием центральных разностей.
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 заданных значений функций, а именно, значения функции y0, y1, y2, ..., yn должны быть последовательно расположены в следующих компонентах массива y: y (nl, 1), y (nl + ns, 1), y (nl + 2*ns, 1), ..., y (nu, 1); расположение в массиве y центральных разностей определяется ниже при описании параметра lo (тип: целый); |
| no - | заданный максимальный порядок разностей, которые используются в конечно - разностных аппроксимациях производных (тип: целый); |
| lo - | определяет порядок расположения заданных центральных разностей в массиве y (тип: целый); при этом если: |
| lo = 0 - | то δ 2j - 1yi - 1/2 должна быть расположена в элементе массива y (nl + i* ns - [ns/2], 2*j), а δ 2jyi - в элементе массива y (nl + i*ns, 2*j + 1); |
| lo = 1 - | то δ 2j - 1yi - 1/2 должна быть расположена в элементе массива y (nl + i* ns - [ns/2], j + 1), а δ 2jyi - в элементе массива y (nl + i*ns, j + 1); |
| lo = 2 - | то δ 2j - 1yi - 1/2 должна быть расположена в элементе массива y(nl + i* ns - [ns/2], j + 1), а δ 2jyi - в элементе массива 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