|
Текст подпрограммы и версий de38r_c.zip , de38d_c.zip |
Тексты тестовых примеров tde38r_c.zip , tde38d_c.zip |
Построение начальных значений при интегрировании системы обыкновенных дифференциальных уравнений второго порядка с правой частью, зависящей от производной, методом Штермера с контролем точности.
Для системы M обыкновенных дифференциальных уравнений второго порядка
Y '' = F( X, Y, Y ' ) ,
Y = ( y1, ..., yM ) ,
F = ( f1( X, y1, ..., yM, y1', ..., yM' ), ..., fM( X, y1, ..., yM, y1', ..., yM' ) ) ,
с начальными условиями, заданными в точке XN:
Y(XN) = YN, YN = ( y10, ..., yM0 ) ,
Y ' (XN) = DYN, DYN = ( y10', ..., yM0' ) ,
отыскиваются необходимые при численном интегрировании этой системы методом Штермера значения решения Y и производной Y ' в первом (после начальной точки XN) узле интегрирования ХN + Н, первая разностная производная назад и конечные разности правой части системы до четвертого порядка включительно в том же узле XN + H.
Четвертый порядок разностей соответствует пятому порядку точности используемого метода Штермера.
Указанная процедура вычисления необходимых для счета начальных значений называется разгоном.
Используемый в подпрограмме метод разгона является итерационным способом, опирающимся на три формулы.
Первая формула представляет значение решения в узле ХN + Н через начальные условия (т.е. через его значение и значение его производной в начальной точке XN), а также через конечные расности правой части системы в точке XN.
Вторая формула является предсказывающей формулой Адамса, по которой вычисляется производная решения.
Третья формула является экстраполяционной формулой, применяемой для предсказания приближенного значения решения в используемом для интегрирования методе Штермера.
При этом приближенные значения вычисляются через значения его разностной производной назад.
B качестве критерия точности при разгоне принята оценка первого отброшенного члена экстраполяционной формулы Штермера, при этом все компоненты решения проверяются на точность по мере погрешности, который заключается в следующем: если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.
Березин И.С., Жидков Н.П. Методы вычислений, т.2. Физматгиз, М., 1960.
Бахвалов H.C., Численные методы, т.1, "Hаука", 1975.
int de38r_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 *dyx, real *z, real *df, real *rfn,
real *rf, real *r, integer *ierr)
Параметры
| f - |
имя подпрограммы вычисления значений правой
части системы. Первый оператор подпрограммы
должен иметь вид:
int f (float *x, float *y, float *dy, float *d2y, int *m) Здесь: x, y, dy - значения независимой, зависимой переменных и производной решения, соответственно; вычисленное значение правой части должно быть помещено в d2y. В случае системы уравнений, т.е. когда m ≠ 1, параметры y, dy и d2y представляют одномерные массивы длины m (тип параметров x, y, dy и d2y: вещественный); |
| m - | количество уравнений в системе (тип: целый); |
| iorder - | порядок точности без единицы того метода Штермера, для которого выполняется разгон и который будет использоваться при интегрировании данной системы уравнений (тип: целый); |
|
xn, yn - dyn | начальные значения аргумента, решения и его производной; в случае системы уравнений (т.е. m ≠ 1) yn и dyn представляют одномерные массивы длиной m (тип: вещественный); |
| hmin - | минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| eps - | допустимая меpа погрешности, с которой тpебуется вычислить все компоненты решения (тип: вещественный); |
| p - | граница перехода, используемая при оценке меры погрешности решения (тип: вещественный); |
| h - | вещественная переменная, содержащая значение шага интегрирования. Если для этого значения шага точность при разгоне достигается, то именно он и реализуется на разгоне, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность; |
| x - | вещественная переменная, значение которой на выходе из подпрограммы представляет первый (после начальной точки xn) узел интегрирования xn + h, в котоpом вычислены необходимые для интегрирования данной системы уравнений начальные значения; |
|
yx - dyx, z | одномерные вещественные массивы длины m, в которых запоминаются значения решения, его производной и первой разностной производной назад, вычисленные в узле x; при этом погрешность решения имеет (iorder + 2) - й порядок по h, а погрешность dyx и z - (iorder + 1) - й порядок по h; |
| df - | двумерный вещественный массив размера m * iorder, в котоpом запоминаются значения правой части системы и ее разностей до (iorder - 1) - го порядка включительно, вычисленные в узле x и умноженные на коэффициент h / 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. |
Версии
| de38d_c - | построение начальных значений при интегрировании системы обыкновенных дифференциальных уpавнений второго порядка, с правой частью, зависящей от производной, методом Штермера с контpолем точности, при этом все промежуточные вычисления выполняются с удвоенным числом значащих цифр. В этом случае параметры xn, yn, dyn, hmin, eps, p, h, x, yx, dyx, z, df, rfn, rf, r и параметры x, y, dy и d2y в подпрограмме f должны иметь тип double. |
Вызываемые подпрограммы
| utde16_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de38r_c. |
| utde17_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de38d_c. |
Замечания по использованию
|
B общем случае требуемая точность при разгоне не гарантируется. Значение параметра h, задаваемое при обращении к подпрограмме, должно быть таким, чтобы узел xn + (iorder - 1) * h не выходил за конец интервала интегрирования. Значение параметра iorder на входе в подпрограмму полагается равным 5. В дальнейшем предполагается расширить допустимое множество значений этого параметра. При работе подпрограммы и ее версии значения параметров m, iorder, xn, yn, dyn, hmin, eps, p сохраняются. |
int main(void)
{
/* Local variables */
extern int de38r_c(U_fp, int *, int *, float *, float *, float *,
float *, float *, float *, float *, float *, float *,
float *, float *, float *, float *, float *, float *,
int *);
static float hmin;
static int ierr;
static float h__;
static int i__;
static float p, r__[2], x, z__[2];
extern int f2_c();
static float df[10] /* was [2][5] */, rf[2], xn, yn[2], yx[2], rfn[2],
eps, dyn[2], dyx[2];
#define df_ref(a_1,a_2) df[(a_2)*2 + a_1 - 3]
yn[0] = 1.f;
yn[1] = -1.f;
dyn[0] = 1.5f;
dyn[1] = -.5f;
xn = 0.f;
p = 100.f;
h__ = .01f;
hmin = 1e-11f;
eps = 1e-4f;
de38r_c((U_fp)f2_c, &c__2, &c__5, &xn, yn, dyn, &hmin, &eps, &p, &h__, &x,
yx, dyx, z__, df, rfn, rf, r__, &ierr);
printf("\n %5i \n", ierr);
printf("\n %16.7e \n", x);
printf("\n %16.7e %16.7e \n", yx[0], yx[1]);
printf("\n %16.7e %16.7e \n", dyx[0], dyx[1]);
printf("\n %16.7e %16.7e \n", z__[0], z__[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 f2_c(float *x, float *y, float *dy, float *z, int *m)
{
/* Parameter adjustments */
--z__;
--dy;
--y;
/* Function Body */
z__[1] = (y[2] + dy[1] - y[1] + dy[2]) * .5f - .5f;
z__[2] = y[1] - *x * .5f;
return 0;
} /* f2_c */
Результаты:
ierr = 0
x = 1.000000000001 - 02
yx(1) = 1.014949833751 yx(2) = - 1.004949833752
dyx(1) = 1.489950167079 dyx(2) = - 4.899501670834 - 01
z__(1) = 1.494983375078 z__(2) = - 4.949833750843 - 01
df_ref(1, 1) = - 8.416248614562 - 04 df_ref(1, 2) = - 8.291528123650 - 06
df_ref(1, 3) = 8.333255330228 - 08 df_ref(1, 4) = 8.377316618179 - 10
df_ref(1, 5) = - 8.498979298110 - 12
df_ref(2, 1) = 8.416248614518 - 04 df_ref(2, 2) = 8.291528119209 - 06
df_ref(2, 3) = - 8.333252932146 - 08 df_ref(2, 4) = - 8.377654125979 - 10
df_ref(2, 5) = 8.512301974406 - 12
Точные значения решения и производной такие:
y1 = 1.014949833745 , y2 = - 1.004949833748
y1' = 1.489950167079 , y2' = - 4.899501670798 - 01