|
Текст подпрограммы и версий de32r_c.zip , de32d_c.zip |
Тексты тестовых примеров tde32r_c.zip , tde32d_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. Количество таких узлов и максимальный порядок получаемых разностей определяются в подпрограмме автоматически и соответствуют порядку точности используемого метода Адамса, который задается пользователем при обращении к подпрограмме.
Указанная процедура вычисления необходимых для счета начальных значений, или "фронта Адамса", называется разгоном. Используемый в подпрограмме метод разгона является итерационным способом, опирающимся на экстраполяционную формулу Адамса, записанную в разностной форме.
B качестве критерия точности при разгоне принята оценка первого отброшенного члена экстраполяционной формулы, при этом все компоненты решения проверяются на точность по меpе погрешности, который заключается в следующем: если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, называемой границей перехода, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.
Беленький В.З. Стандартная программа для интегрирования системы дифференциальных уравнений методом Адамса. Сб. "Вычислительные методы и программирование", вып.3., Изд-во МГУ, 1965.
int de32r_c (S_fp f, integer *m, integer *iorder, real *xn,
real *yn, real *hmin, real *eps, real *p, real *h, real *df,
real *x, real *y, logical *bul, 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 (тип: вещественный); |
| hmin - | минимальное значение абсолютной величины шага, который разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| eps - | допустимая меpа погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
| p - | граница перехода, используемая при оценке погрешности решения (тип: вещественный); |
| 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; |
| bul - | логическая переменная, значение которой при обращении к подпрограмме полагается равным TRUE_, если заданный в h шаг выводит в конец интервала интегрирования,т.е. узел x (iorder) совпадает с концом интервала интегрирования, и FALSE_ в противном случае. B результате работы подпрограммы bul pавно FALSE_, если вместо исходного шага интегрирования при разгоне был использован меньший шаг; в противном случае значение переменной bul не меняется; |
|
rs - rf, r | одномерные вещественные рабочие массивы длины m; |
| ierr - | целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом: |
| ierr= 1 - | когда неправильно задан параметр iorder,т.е. iorder > 10; в этом случае разгон выполняется для значения iorder = 10; |
| ierr=65 - | когда на разгоне не может быть достигнута требуемая точность eps; в этом случае разгон можно начать сначала обращение к подпрограмме с новыми значениями параметров h, hmin и iorder. |
Версии
| de32d_c - | построение начальных значений с повышенной точностью при интегрировании системы обыкновенных дифференциальных уравнений первого порядка методом Адамса с контролем точности. При этом параметры xn, yn, hmin, eps, p, h, df, x, y, rs, rf, r и параметры x, y и dy в подпрограмме f должны иметь тип double. |
Вызываемые подпрограммы
| utde12_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de32r_c. |
| utde13_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de32d_c. |
| Kpоме того, подпрограммы de32r_c и de32d_c используют рабочие подпрограммы de28rs_c и de28ds_c, соответственно. |
Замечания по использованию
|
B общем случае требуемая точность при разгоне не гарантируется. При работе подпрограммы и ее версии значения параметров m, iorder, xn, yn, hmin, eps, p сохраняются. Подпрограмма и ее версия могут использоваться как для разгона, так и непосредственно для численного решения системы уравнений на всем интервале интегрирования. |
y1' = y2
y2' = -y1
y1 (3/4 π) = √2 /2
y2 (3/4 π) = -√2 /2
int main(void)
{
/* Builtin functions */
double sqrt(double);
/* Local variables */
extern int de32r_c(U_fp, int *, int *, float *, float *, float *,
float *, float *, float *, float *, float *, float *,
logical *, float *, float *, float *, int *);
static float hmin;
static int ierr;
extern int f_c();
static float h__;
static int i__, m;
static float p, r__[2], x[5], y[10] /* was [2][5] */,
df[10] /* was [2][5] */,
rf[2], rs[2], xn, yn[2];
static int iorder;
static logical bul;
static float eps;
#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]
bul = FALSE_;
iorder = 5;
p = 100.f;
eps = 1e-4f;
hmin = 1e-10f;
m = 2;
h__ = .01f;
xn = 2.3561944901925003f;
yn[0] = (float)sqrt(2.f) / 2.f;
yn[1] = -yn[0];
de32r_c((U_fp)f_c, &m, &iorder, &xn, yn, &hmin, &eps, &p, &h__, df, x, y,
&bul, rs, rf, r__, &ierr);
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));
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