|
Текст подпрограммы и версий de34r_c.zip , de34d_c.zip , de18r_c.zip , de18d_c.zip |
Тексты тестовых примеров tde34r_c.zip , tde34d_c.zip , tde18r_c.zip , tde18d_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 '(XN) = DYN, DYN = ( y10', ..., yM0' ) ,
отыскиваются необходимые при численном интегрировании этой системы многошаговым методом Штермера значение решения Y в первом (после начальной точки XN) узле интегрирования ХN + Н, его конечная разность первого порядка и конечные разности правой части системы до четвертого порядка включительно в том же узле XN + H.
Четвертый порядок разностей соответствует шестому порядку точности используемого метода Штермера.
Указанная процедура вычисления необходимых для счета начальных значений называется разгоном.
Используемый в подпрограмме метод разгона является итерационным способом, опирающимся на пару формул, одна из которых представляет значение решения в узле XN + H через начальные условия (т.е. через его значение и значение его производной в начальной точке XN), а также через конечные разности правой части системы в точке XN, а другая является экстраполяционной формулой, применяемой для предсказания приближенного значения решения в используемом для интегрирования методе Штермера.B качестве критерия точности при разгоне принята оценка первого отброшенного члена экстраполяционной формулы Адамса - Штермера, при этом все компоненты решения проверяются на точность по меpе погрешности, который заключается в следующем: если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, называемой границей перехода, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.
Березин И.С., Жидков Н.П. Методы вычислений, т.2. Физматгиз, М., 1960.
int de34r_c (S_fp f, integer *m, integer *iorder, real *xn,
real *yn, real *dyn, real *hmin, real *eps, real *p, real *h,
real *x, real *yx, real *dy, real *df, real *rfn, real *rf,
real *r, integer *ierr)
Параметры
| f - | имя подпрограммы вычисления значений правой части системы. Первый оператор подпрограммы должен иметь вид: |
|
int f_c(float *x, float *y, float *dy, int *m) Здесь: x, y - значения независимой и зависимой переменных, соответственно; вычисленное значение правой части должно быть помещено в dy. В случае системы уравнений, т.е. когда m ≠ 1, параметры y и dy представляют одномерные массивы длиной m (тип параметров x, y и dy: вещественный); |
| m - | количество уравнений в системе (тип: целый); |
| iorder - | порядок точности без единицы того метода Штермера, для которого выполняется разгон и который будет использоваться при интегрировании данной системы уравнений (тип: целый); |
|
xn, yn - dyn | начальные значения аргумента, решения и его производной; в случае системы уравнений (т.е. m ≠ 1) yn и dyn представляют одномерные массивы длиной m (тип: вещественный); |
| hmin - | минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| eps - | допустимая меpа погрешности, с которой тpебуется вычислить все компоненты решения (тип: вещественный); |
| p - | граница перехода, используемая при оценке погрешности решения (тип: вещественный); |
| h - | вещественная переменная, содержащая значение шага интегрирования. Если для этого значения шага точность при разгоне достигается, то именно он и реализуется на разгоне, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность; |
| x - | вещественная переменная, значение которой на выходе из подпрограммы представляет первый (после начальной точки xn) узел интегрирования xn + h, в котоpом вычислены необходимые для интегрирования данной системы уравнений начальные значения; |
| yx, dy - | одномерные вещественные массивы длиной m, в которых запоминаются значения решения и его первой разности, вычисленные в узле x, при этом погрешность решения и разности имеет (iorder + 1) - й порядок по h; |
| df - | двумерный вещественный массив размера m * iorder, в котоpом запоминаются значения правой части системы и ее разностей до (iorder - 1) - го порядка включительно, вычисленные в узле x и умноженные на коэффициент h2 / 12, при этом элемент df (i, 1) этого массива содержит значение правой части i - го уравнения системы, а df (i, j + 1) - ее j - ю разность, погрешность которой имеет (iorder + 1) - й порядок по h; |
|
rfn - rf, r | одномерные вещественные рабочие массивы длины m; |
| ierr - | целая переменная, значение которой в pезультате работы подпрограммы полагается равным 65, если на разгоне не может быть достигнута требуемая точность eps; в этом случае разгон можно начать сначала обращением к подпрограмме с новыми значениями паpаметpов h и hmin. |
Версии
| de34d_c - | построение начальных значений с повышенной точностью при интегрировании системы обыкновенных дифференциальных уравнений второго порядка методом Штермера с контролем точности. При этом параметры xn, yn, dyn, hmin, eps, p, h, x, yx, dy, df, rfn, rf, r и параметры x, y и dy в подпрограмме должны иметь тип double. |
| de18r_c - | назначение этой подпрограммы такое же, как и подпрограммы de34r_c, но схема организации вычислений в ней несколько иная, чем в de34r_c, а именно: значение решения вычисляется через значение некоторой промежуточной переменной, введение которой позволяет уменьшить вычислительную погрешность в приближенном значении решения. Все формальные параметры de18r_c имеют тот же смысл, что и соответствующие им параметры de34r_c с той только разницей, что в массив dy заносится разделенная разность решения, вычисленная с порядком iorder по h, а не конечная разность как в de34r_c, а в массив df - конечные разности правой части системы, умноженные на коэффициент h / 12 и вычисленные с порядком iorder по h. |
| de18d_c - | то же, что и de18r_c, но все вычисления ведутся с удвоенным числом значащих цифр; при этом параметры xn, yn, dyn, hmin, eps, p, h, x, yx, dy, df, rfn, rf, r и параметры x, y и dy в подпрограмме f должны иметь тип double. |
Вызываемые подпрограммы
| utde12_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de34r_c. |
| utde13_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de34d_c. |
| utde16_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de18r_c. |
| utde17_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de18d_c. |
| Kpоме того, подпрограммы de34r_c и de34d_c используют рабочие подпрограммы de28rs_c и de28ds_c, соответственно. |
Замечания по использованию
|
B общем случае требуемая точность при разгоне не гарантируется. Значение праметра h, задаваемое при обращении к подпрограмме, должно быть таким, чтобы узел xn + (iorder - 1) * h не выходил за конец интервала интегрирования. Значение параметра iorder на входе в подпрограмму полагается равным 5. B дальнейшем предполагается расширить допустимое множество значений этого параметра. При работе подпрограммы и ее версии значения параметров m, iorder, xn, yn, dyn, hmin, eps, p сохраняются. |
y1'' = - y1 ,
y2'' = - y2 ,
y1 (3/4 * π) = √2 / 2, y2 (3/4 * π) = - √2 / 2 ,
y1' (3/4 * π) = - √2 / 2, y2' (3/4 * π) = - √2 / 2 .
int main(void)
{
/* Builtin functions */
double sqrt(double);
/* Local variables */
extern int de34r_c(U_fp, int *, int *, float *, float *, float *,
float *, float *, float *, float *, float *, float *,
float *, float *, float *, float *, float *, int *);
static float hmin;
static int ierr;
extern int f_c();
static float h__;
static int i__, m;
static float p, x, y[2], delty[2], df[10] /* was [2][5] */,
rf[2], dy[2], xn, yn[2];
static int iorder;
static float rfn[2], eps, dyn[2];
#define df_ref(a_1,a_2) df[(a_2)*2 + a_1 - 3]
h__ = .01f;
iorder = 5;
p = 100.f;
eps = 1e-4f;
hmin = 1e-10f;
m = 2;
xn = 2.3561944901925003f;
yn[0] = (float)sqrt(2.f) / 2.f;
yn[1] = -yn[0];
dyn[0] = yn[1];
dyn[1] = -yn[0];
de34r_c((U_fp)f_c, &m, &iorder, &xn, yn, dyn, &hmin, &eps, &p, &h__, &x,
y, dy, df, rfn, rf, delty, &ierr);
printf("\n %16.7e \n", x);
printf("\n %16.7e %16.7e \n", y[0], y[1]);
printf("\n %16.7e %16.7e \n\n", dy[0], dy[1]);
for (i__ = 1; i__ <= 2; ++i__) {
printf("\n %16.7e %16.7e %16.7e ",
df_ref(i__, 1), df_ref(i__, 2), df_ref(i__, 3));
printf("\n %16.7e %16.7e \n",
df_ref(i__, 4), df_ref(i__, 5));
}
return 0;
} /* main */
int f_c(float *t, float *y, float *z, int *m)
{
/* Parameter adjustments */
--z__;
--y;
/* Function Body */
z__[1] = -y[1];
z__[2] = -y[2];
return 0;
} /* f_c */
Результаты:
ierr = 0
x = 2.366194490 + 00
y(1) = 7.000004762 - 01 y(2) = - 7.141423761 - 01
dy(1) = - 7.106305006 - 03 dy(2) = - 7.035594917 - 03
df_ref(1, 1) = - 5.83333730151 - 06 df_ref(1, 2) = 5.92192083634 - 08
df_ref(1, 3) = 5.89251453187 - 10 d_ref(1, 4) = - 5.86490578325 - 12
df_ref(1, 5) = - 5.76760861293 - 14
df_ref(2, 1) = 5.95118646749 - 06 df_ref(2, 2) = 5.86299576250 - 08
df_ref(2, 3) = - 5.89250072347 - 10 df_ref(2, 4) = - 5.92374760355 - 12
df_ref(2, 5) = 6.01324545713 - 14