|
Текст подпрограммы и версий de35r_c.zip , de35d_c.zip |
Тексты тестовых примеров tde35r_c.zip , tde35d_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, а другая является экстраполяционной формулой, применяемой для предсказания приближенного значения решения в используемом для интегрирования методе Штермера.
Березин И.С., Жидков Н.П. Методы вычислений, т.2. Физматгиз, М., 1960.
int de35r_c (S_fp f, integer *m, integer *iorder, real *xn,
real *yn, real *dyn, real *h, real *x, real *yx, real *dy, real *df,
real *rfn, real *rf, real *r)
Параметры
| f - | имя подпрограммы вычисления значений правой части системы. Первый оператор подпрограммы должен иметь вид: |
|
int f (float *x, float *y, float *dy, int *m). Здесь: x, y - значения независимой и зависимой переменных, соответственно; вычисленное значение правой части должно быть помещено в dy. B случае системы уравнений, т.е. когда m ≠ 1, параметры y и dy представляют одномерные массивы длиной m (тип параметров x, y и dy: вещественный); |
| m - | количество уравнений в системе (тип: целый); |
| iorder - | порядок точности без единицы того метода Штермера, для которого выполняется разгон и который будет использоваться при интегрировании данной системы уравнений (тип: целый); |
|
xn, yn - dyn | начальные значения аргумента, решения и его производной; в случае системы уравнений (т.е. m ≠ 1) yn и dyn представляют одномерные массивы длиной m (тип: вещественный); |
| 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. |
Версии
| de35d_c - | построение начальных значений при интегрировании системы обыкновенных дифференциальных уpавнений второго порядка методом Штермера без контроля точности, при этом все вычисления с действительными числами выполняются с удвоенным числом значащих цифр; в этом случае параметры xn, yn, dyn, h, x, yx, dy, df, n, rf, r и параметры x, y и dy в подпрограмме f должны иметь тип double. |
Вызываемые подпрограммы
| Подпрограммы de35r_c и de35d_c используют рабочие подпрограммы de28rs_c и de28ds_c, соответственно. |
Замечания по использованию
|
Значение параметра h, задаваемое при обращении к подпрограмме, должно быть таким, чтобы узел xn + (iorder - 1) * h не выходил за конец интервала интегрирования. Значение параметра iorder на входе в подпрограмму полагается равным 5. В дальнейшем предполагается расширить допустимое множество значений этого параметра. При работе подпрограммы и ее версии значения параметров m, iorder, xn, yn, dyn, h сохраняются. |
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 de35r_c(U_fp, int *, int *, float *, float *, float *,
float *, float *, float *, float *, float *, float *,
float *, float *);
extern int f_c();
static float h__;
static int i__, m;
static float x, y[2], delty[2], df[10] /* was [2][5] */,
rf[2], dy[2], xn, yn[2];
static int iorder;
static float rfn[2], dyn[2];
#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];
dyn[0] = yn[1];
dyn[1] = -yn[0];
de35r_c((U_fp)f_c, &m, &iorder, &xn, yn, dyn, &h__, &x, y, dy, df, rfn,
rf, delty);
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 */
Результаты:
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 df_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