|
Текст подпрограммы и версий de40r_c.zip , de40d_c.zip |
Тексты тестовых примеров tde40r_c.zip , tde40d_c.zip |
Выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений второго порядка методом Штермера без контроля точности.
Выполняется один шаг численного интегрирования системы М обыкновенных дифференциальных уравнений второго порядка
Y '' = F (X, Y) ,
Y = ( y1, ..., yM ) ,
F = ( f1 (X, y1, ..., yM), ..., fM (X, y1, ..., yM) )
многошаговым предсказывающе - исправляющим методом пятого порядка точности. Суть метода состоит в том, что сначала по явной формуле Адамса - Штермера строятся предсказанные значения приближенного решения, которые затем исправляются по неявной формуле Адамса - Штермера.
По значению решения YX в узле Xn вычиляется значение этого решения в узле Xn + H. При этом, если Xn является началом ХN интервала интегрирования, то для вычисления решения в XN + H требуется еще и знание производной решения в XN.
Bместе с решением подпрограмма определяет также всю информацию, необходимую для вычисления приближенного решения в любом следующем узле. Значение H шага интегрирования задается пользователем при обращении к подпрограмме.
Березин И.С., Жидков Н.П. Методы вычислений, т.2. Физматгиз, М., 1960.
Жоголев E.A. Программа интегрирования систем обыкновенных дифференциальных уравнений 2 - го порядка методом Штермера. Сб. "Вычислительные методы и программирование", вып. 1 , Изд - во МГУ, 1962.
int de40r_c (S_fp f, integer *m, integer *jstart, real *yx,
real *x, real *dy, real *h, real *delty, real *df, real *rfn,
real *rf, real *yp, real *dyp, integer *ierr)
Параметры
| f - |
имя подпрограммы вычисления значений
правой части системы дифференциальных
уpавнений. Первый оператор подпрограммы должен иметь вид: int f (float *x, float *y, float *dy, int *m). Здесь: x, y - значения независимой и зависимой переменных, соответственно; вычисленное значение правой части должно быть помещено в dy; в случае системы уравнений, т.е. когда m ≠ 1, параметры y и dy представляют одномерные массивы длиной m (тип параметров x, y, dy: вещественный); |
| m - | количество уравнений в системе (тип: целый); |
| jstart - | целый указатель режима использования подпрограммы, имеющий следующие значения: |
| 0 - | первое обращение к подпрограмме должно быть выполнено с нулевым значением jstart; |
| +1 - | выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных и шага интегрирования, заданных параметрами x, yx и h, соответственно; |
| -1 - | повторить последний шаг интегрирования с новым значением параметра h; |
| на выходе из подпрограмы jstart полагается равным + 1; |
| x, yx - | заданные вещественные значения аргумента и решения в нем; в результате работы подпрограммы в x получается новое значение аргумента, а в yx - соответствующее значение решения; в случае системы уравнений, т.е. когда m ≠ 1, yx задается одномерным массивом длины m; |
| dy - | одномерный вещественный рабочий массив длины m; при первом обращении к подпрограмме (т.е. со значением jstart = 0) в этом массиве задается значение производной pешения в начальном узле интервала интегрирования; |
| h - | вещественное значение шага интегрирования; |
|
delty - rfn, rf yp, dyp | одномерные вещественные рабочие массивы длины m; |
| df - | двумерный вещественный рабочий массив размеpа m * 5, в котоpом запоминаются значения правой части системы и ее разностей до четвертого порядка включительно, умноженные на h 2 / 12; |
| ierr - | целая переменная, значение которой в pезультате работы подпрограммы полагается равным 65, если при заданном значении шага интегрирования не может быть достигнута сходимость численного решения к точному значению решения; в этом случае интегрирование системы можно повторить обращением к подпрограмме с новым значением параметра h и со значением jstart = + 1. |
Версии
| de40d_c - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений второго порядка методом Штермера без контроля точности, при этом все вычисления над действительными числами ведутся с удвоенным числом значащих цифр. В этом случае параметры yx, x, dy, h, delty, df, rfn, rf, yp, dyp и параметры x, y и dy в подпрограмме f должны иметь тип double. |
Вызываемые подпрограммы
| de35r_c - | построение начальных значений при интегрировании системы обыкновенных дифференциальных уравнений второго порядка методом Штермера без контроля точности; вызывается при работе подпрограммы de40r_c. |
| de35d_c - | построение начальных значений при интегрировании системы обыкновенных дифференциальных уравнений второго порядка методом Штермера без контроля точности, при этом все вычисления над действительными числами ведутся с удвоенным числом значащих цифр; вызывается при работе подпрограммы de40d_c. |
| utde12_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de40r_c. |
| utde13_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de40d_c. |
| Kpоме того, de40r_c и de40d_c используют рабочие подпрограммы de42rp_c, de28rs_c и de42dp_c, de28ds_c, соответственно. |
Замечания по использованию
|
При работе подпрограммы и ее версии значения параметров m и h сохраняются. При многократном использовании подпрограммы и ее версии для вычисления решения на отрезке значения параметров m, yx, x, dy, df, yp, dyp не должны изменяться в вызывающей программе между последовательными обращениями к ним. |
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), sin(double), cos(double);
/* Local variables */
extern int de40r_c(U_fp, int *, int *, float *, float *, float *,
float *, float *, float *, float *, float *, float *,
float *, int *);
static int ierr;
extern int f_c();
static float h__;
static int m;
static float x, y[2], delty[2], y1, y2, df[10];
static int ih;
static float rf[2], dy[2], yp[2];
static int jstart;
static float rfn[2], dyp[2];
m = 2;
x = 2.3561944901925003f;
h__ = .01f;
y[0] = (float)sqrt(2.f) / 2.f;
y[1] = -y[0];
dy[0] = y[1];
dy[1] = -y[0];
jstart = 0;
ih = 0;
h__ = -h__;
l6:
de40r_c((U_fp)f_c, &m, &jstart, y, &x, dy, &h__, delty, df, rfn, rf, yp,
dyp, &ierr);
++ih;
y1 = (float)sin(x);
y2 = (float)cos(x);
printf("\n %16.7e \n", x);
printf("\n %16.7e %16.7e \n", y[0], y[1]);
printf("\n %16.7e \n", h__);
printf("\n %16.7e %16.7e \n", y1, y2);
if (ih == 1) {
goto l6;
}
if (ih == 5) {
goto l20;
}
jstart = -1;
if (ih == 4) {
goto l8;
}
if (ih == 3) {
goto l7;
}
h__ = .005f;
h__ = -h__;
goto l6;
l7:
h__ = .03f;
goto l6;
l8:
h__ = .01f;
goto l6;
l20:
return 0;
} /* main */
int f_c(float *x, float *y, float *dy, int *m)
{
/* Parameter adjustments */
--dy;
--y;
/* Function Body */
dy[1] = -y[1];
dy[2] = -y[2];
return 0;
} /* f_c */
Результаты:
после первого обращения к подпрограмме -
x y(1) y(2)
2.366194490 + 00 7.000004762 - 01 - 7.141423761 - 01
h__ y1 y2
1.000000000 - 02 7.000004762 - 01 - 7.141423761 - 01
после второго обращения к подпрограмме -
x y(1) y(2)
2.376194490 + 00 6.928241717 - 01 - 7.211065574 - 01
h__ y1 y2
1.000000000 - 02 6.928241717 - 01 - 7.211065574 - 01
после третьего обращения к подпрограмме -
x y(1) y(2)
2.361194490 + 00 7.035624232 - 01 - 7.106334615 - 01
h__ y1 y2
- 5.000000000 - 03 7.035624232 - 01 - 7.106334615 - 01
после четвертого обращения к подпрограмме -
x y(1) y(2)
2.386194490 + 00 6.855785854 - 01 - 7.279986286 - 01
h__ y1 y2
2.000000000 - 02 6.855785854 - 01 - 7.279986286 - 01
после пятого обращения к подпрограмме -
x y(1) y(2)
2.376194490 + 00 6.928241717 - 01 - 7.211065574 - 01
h__ y1 y2
1.000000000 - 02 6.928241717 - 01 - 7.211065574 - 01