|
Текст подпрограммы и версий de33r_c.zip , de33d_c.zip |
Тексты тестовых примеров tde33r_c.zip , tde33d_c.zip |
Построение начальных значений при интегрировании системы обыкновенных дифференциальных уравнений первого порядка методом Адамса без контроля точности.
Для системы M обыкновенных дифференциальных уравнений первого порядка
Y ' = F (X, Y) ,
Y = ( y1, ..., yM ) ,
F = ( f1 (X, y1,..., yM),..., fM (X, y1,..., yM) )
с начальными условиями, заданными в точке XN:
Y(XN) = YN , YN = ( y10,..., yM0 ) ,
с помощью модифицированного приема А.Н.Крылова отыскиваются необходимые при численном интегрировании этой системы многошаговым методом Адамса значения решения Y в нескольких первых узлах и конечные разности правой части F (X,Y) системы до некоторого порядка в точке XN. Количество таких узлов и максимальный порядок получаемых разностей определяются в подпрограмме автоматически и соответствуют порядку точности используемого метода Адамса, который задается пользователем при обращении к подпрограмме.
Указанная процедура вычисления необходимых для счета начальных значений, или "фронта Адамса", называется разгоном. Используемый в подпрограмме метод разгона является итерационным способом, опирающимся на экстраполяционную формулу Адамса, записанную в разностной форме.
Беленький В.З. Стандартная программа для интегрирования системы дифференциальных уравнений методом Адамса. Сб. "Вычислительные методы и программирование", вып. 3, Изд-во МГУ, 1965.
int de33r_c (S_fp f, integer *m, integer *iorder, real *xn,
real *yn, real *h, real *df, real *x, real *y, real *rs, real *rf,
real *r, 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 - | количество уравнений в системе (тип: целый); |
| iorder - | порядок точности того метода Адамса, для которого выполняется разгон и который будет использоваться при интегрировании данной системы уравнений; iorder должен быть не больше 10 (тип: целый); |
| xn, yn - | начальные значения аргумента и решения; в случае системы уравнений (т.е. m ≠ 1) yn представляет одномерный массив длины m (тип: вещественный); |
| h - | вещественное значение шага интегрирования; |
| df - | двумерный вещественный массив размера m*iorder, в котоpом запоминаются значения правой части системы и ее разностей до порядка (iorder - 1) включительно, вычисленные в точке xn, при этом элемент df (i, 1) этого массива содержит значение правой части i - ого уравнения, а df (i,j + 1) - ее j - ю разность, погрешность вычисления которой имеет (iorder + 1) - й порядок по h; |
| x - | одномерный вещественный массив длины iorder, содержащий на выходе из подпрограммы iorder узлов, включая xn, в которых вычисляются при разгоне приближенные значения решения, при этом элемент x (j) этого массива содержит узел, равный xn + (j - 1) * h; |
| y - | двумерный вещественный массив размера m*iorder, в котоpом запоминаются приближенные значения решения, вычисленные для значений аргумента, хранящихся в массиве x, а именно, значению аргумента в x (j) соответствует приближенное решение y (i,j), при этом погрешность этого решения имеет порядок iorder по h; |
|
rs - rf, r | одномерные вещественные рабочие массивы длины m; |
| ierr - | целая переменная, значение которой в результате работы подпрограммы полагается равным 1, если неправильно задан параметр iorder, т.е. iorder > 10; в этом случае разгон выполняется для iorder = 10. |
Версии
| de33d_c - | построение начальных значений при интегрировании системы обыкновенных дифференциальных уpавнений первого порядка методом Адамса без контроля точности, при этом все вычисления с действительными числами выполняются с удвоенным числом значащих цифр; в этом случае параметры xn, yn, h, df, x, y, rs, rf, r и параметры x, y и dy в подпрограмме f должны иметь тип double. |
Вызываемые подпрограммы
| utde12_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de33r_c. |
| utde13_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de33d_c. |
| Kpоме того, подпограммы de33r_c и de33d_c используют рабочие подпрограммы de28rs_c и de28ds_c, соответственно. |
Замечания по использованию
|
При работе подпрограммы и ее версии значения параметров m, iorder, xn, yn, h сохраняются. Подпрограмма и ее версия могут использоваться как для разгона, так и непосредственно для численного решения системы уравнений на всем интервале интегрирования. |
y1' = y2
y2' = -y1
y1 (3/4 π) = √2 /2
y2 (3/4 π) = -√2 /2
int main(void)
{
/* Builtin functions */
double sqrt(double), sin(double), cos(double);
/* Local variables */
extern int de33r_c(U_fp, int *, int *, float *, float *, float *, float *,
float *, float *, float *, float *, float *, int *);
static int ierr;
extern int f_c();
static float h__;
static int i__, m;
static float r__[2], x[5], y[10] /* was [2][5] */, r1, r2,
df[10] /* was [2][5] */,
rf[2], rs[2], xn, yn[2];
static int iorder;
#define y_ref(a_1,a_2) y[(a_2)*2 + a_1 - 3]
#define df_ref(a_1,a_2) df[(a_2)*2 + a_1 - 3]
iorder = 5;
m = 2;
h__ = .01f;
xn = 2.3561944901925003f;
yn[0] = (float)sqrt(2.f) / 2.f;
yn[1] = -yn[0];
de33r_c((U_fp)f_c, &m, &iorder, &xn, yn, &h__, df, x, y, rs, rf, r__,
&ierr);
xn = x[4];
r1 = (float)sin(xn);
r2 = (float)cos(xn);
for (i__ = 1; i__ <= 5; ++i__) {
printf("\n %16.7e \n", x[i__-1]);
}
printf("\n %16.7e %16.7e \n", y_ref(1, 5), y_ref(2, 5));
printf("\n %16.7e %16.7e \n", r1, r2);
for (i__ = 1; i__ <= 2; ++i__) {
printf("\n %16.7e %16.7e %16.7e \n",
df_ref(i__, 1), df_ref(i__, 2), df_ref(i__, 3));
printf(" %16.7e %16.7e \n",
df_ref(i__, 4), df_ref(i__, 5));
}
return 0;
} /* main */
int f_c(float *x, float *y, float *dy, int *m)
{
/* Parameter adjustments */
--dy;
--y;
/* Function Body */
dy[1] = y[2];
dy[2] = -y[1];
return 0;
} /* f_c */
Результаты:
ierr = 0
x(5) = 2.396194490 + 00
y(1, 5) = 6.782644418-01
y(2, 5) = -7.348179005-01
df_ref(1, 1) = -7.07106781186-03
df_ref(1, 2) = -7.10630480683-05
df_ref(1, 3) = 6.99988802921-07
df_ref(1, 4) = 7.18290493751-09
df_ref(1, 5) = -7.27808924239-11
df_ref(2, 1) = -7.07106781186-03
df_ref(2, 2) = 7.03559472512-05
df_ref(2, 3) = 7.14142103675-07
df_ref(2, 4) = -6.97064450605-09
df_ref(2, 5) = -6.86739554112-11