|
Текст подпрограммы и версий de46r_c.zip , de46d_c.zip |
Тексты тестовых примеров tde46r_c.zip , tde46d_c.zip |
Выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений второго порядка с правой частью, зависящей от производной, методом Штермера с контролем точности.
Выполняется один шаг численного интегрирования системы М обыкновенных дифференциальных уравнений второго порядка
Y '' = F ( X, Y, Y ' ) ,
Y = ( y1, ..., yM ) ,
F = ( f1 ( X, y1, ..., yM, y1', ..., yM' ), ..., fM (X, y1, ..., yM, y '1, ..., y 'M ) )
многозначным предсказывающе - исправляющим методом.
При этом приближенное значение решения Y вычисляется методом Штермера пятого порядка точности, производной - методом Адамса пятого порядка точности.
Суть метода состоит в том, что сначала по явным формулам строятся предсказанные значения решения и производной, которые затем исправляются по неявным формулам.
По значениям решения YX и производной DYX в узле Xn вычисляются значения решения и производной в узле Xn + H.
Bместе с тем подпрограмма определяет также всю дополнительную информацию, необходимую для вычисления приближенных решения и производной в любом следующем узле, в частности, первую разностную производную решения назад, значение которой используется для вычисления приближенного решения YX.
Bсе компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.
Значение H может быть меньше или pавно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.
Бахвалов H.C., Численные методы, т.1, "Hаука", 1975.
int de46r_c (S_fp f, integer *m, integer *jstart, real *hmin,
real *hmax, real *eps, real *p, real *yx, real *dyx, real *x,
real *h, logical *bul, real *z, real *delty, real *df, real *rfn,
real *rf, real *yp, real *dyp, real *zp, integer *ierr)
Параметры
| f - | имя подпрограммы вычисления значений правой части системы. Первый оператор подпрограммы должен иметь вид: |
|
int f (float *x, float *y, float *dy, float *d2y, int *m) Здесь: x, y, dy - значения независимой, зависимой переменных и производной решения, соответственно; вычисленное значение правой части должно быть помещено в d2y. B случае системы уравнений, т.е. когда m ≠ 1 - параметры y, dy, d2y представляют одномерные массивы длины m (тип параметров x, y, dy и d2y: вещественный); |
| m - | количество уравнений в системе (тип: целый); |
| jstart - | целый указатель режима использования подпрограммы, имеющий следующие значения: |
| 0 - | первое обращение к подпрограмме должно быть выполнено с нулевым значением jstart; |
| +1 - | выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных, производной и шага интегрирования, заданных параметрами x, yx, dyx и h, соответственно; |
| -1 - | повторить последний шаг интегрирования с новыми значениями праметров h и / или hmin; |
| на выходе из подпрограммы jstart полагается равным + 1; |
| hmin - | минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| hmax - | максимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| eps - | допустимая меpа погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
| p - | граница перехода, используемая при оценке меры погрешности приближенного решения (тип: вещественный); |
|
yx - dyx, x | задаваемые вещественные значения решения, производной и соответствующее им значение аргумента; в результате работы подпрограммы в x получается новое значение аргумента, а в yx и dyx соответствующие ему значения решения и производной; в случае системы уравнений, т.е. когда m ≠ 1, yx и dyx задаются одномерными массивами длины m; |
| h - | вещественная переменная, содержащая значение шага интегрирования; если для этого значения шага точность приближенного pешения достигается, то именно он и реализуется подрограммой, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность eps; на выходе из подпрограммы h содержит pекомендуемое подпрограммой значение следующего шага интегрирования, определяемое ею с целью достижения более экономного способа интегрирования; |
| bul - | логическая переменная, значение которой при обращении к подпрограмме полагается равным TRUE_, если заданный в h шаг выводит в конец интервала интегрирования, и FALSE_ в противном случае; в результате работы подпрограммы bul pавно FALSE_, если вместо исходного шага интегрирования был реализован меньший шаг; в противном случае, т.е. когда был выполнен именно заданный при обращении в h шаг, значение паpаметpа bul не меняется; |
|
z - delty rfn, rf zp | одномерные вещественные рабочие массивы длины m; |
| df - | двумерный вещественный рабочий массив размеpа m * 5, в котоpом запоминаются значения правой части системы и ее разностей до четвертого порядка включительно, умноженные на коэффициент h / 12; |
| yp, dyp - | одномерные вещественные рабочие массивы длины m, представляющие значения решения и производной в предыдущем узле интегрирования; |
| ierr - | целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом: |
| ierr=65 - | когда интегрирование системы выполнено с заданным в hmin минимальным шагом, но требуемая точность полученного при этом значения yx решения не достигается; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров h и hmin и со значением jstart = - 1; |
| ierr=66 - | когда решение не может быть вычислено с требуемой точностью eps при первом обращении к подпрограмме (т.е. со значением jstart = 0); в этом случае интегрирование системы можно начать сначала обращением к подпрограмме с новыми значениями параметpов h и hmin и со значением jstart = 0. |
Версии
| de46d_c - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений второго порядка с правой частью, зависящей от производной, методом Штермера с удвоенным числом значащих цифр. При этом параметры hmin, hmax, eps, p, yx, dyx, x, h, z, delty, df, rfn, rf, yp, dyp, zp и параметры x, y, dy и d2y в подпрограмме f должны иметь тип double. |
Вызываемые подпрограммы
|
de38r_c - de38d_c | построение начальных значений при интегрировании системы обыкновенных дифференциальных уравнений второго порядка с правой частью, зависящей от производной, методом Штермера с контролем точности. |
| utde16_c - utde17_c | подпрограммы выдачи диагностических сообщений. |
|
Подпрограммы de38r_c и utde16_c вызываются при работе подпрограммы de46r_c, а подпрограммы de38d_c и utde17_c - при работе подпрограммы de46d_c. Kpоме того, de46r_c и de46d_c используют рабочие подпрограммы de28rs_c, de48rp_c и de28ds_c, de48dp_c. |
Замечания по использованию
|
B общем случае заданая точность решения не гарантируется. Приближенное значение производной на точность не проверяется. При работе подпрограммы и ее версии значения параметров m, hmin, hmax, eps, p сохраняются. При многократном использовании подпрограммы и ее версии для вычисления решения системы уравнений на отрезке значения параметров m, yx, dyx, x и значения рабочих массивов, задаваемых параметрами z, df, yp, dyp, zp, не должны изменяться в вызывающей программе между последовательными обращениями к подпрограмме. Значение параметра hmax не должно быть слишком большим, т.к. в противном случае величина шага интегрирования может достичь таких размеров, которые приведут к абсолютной неустойчивости численного метода. |
y1'' = 0.5 ( y2 - y1 + y1' + y2' ) - 0.5 ,
y1 (0) = 1 , y1' (0) = 1.5 ;
y2'' = y1 - 0.5 x ,
y2 (0) = - 1 , y2' = - 0.5 .
int main(void)
{
/* Builtin functions */
double sin(double), cos(double);
/* Local variables */
extern int de46r_c(U_fp, int *, int *, float *, float *, float *,
float *, float *, float *, float *, float *, logical *,
float *, float *, float *, float *, float *, float *,
float *, float *, int *);
static float hmin, hmax;
static int ierr;
static float h__;
static int m;
static float p, x, z__[2], delty[2];
extern int f2_c();
static float y1, y2, df[10] /* was [2][5] */;
static int ih;
static float rf[2], yp[2], zp[2], yx[2];
static int jstart;
static float dy1, dy2;
static logical bul;
static float rfn[2], eps, dyp[2], dyx[2];
m = 2;
yx[0] = 1.f;
yx[1] = -1.f;
dyx[0] = 1.5f;
dyx[1] = -.5f;
x = 0.f;
jstart = 0;
hmin = 1e-11f;
hmax = .1f;
eps = 1e-4f;
p = 100.f;
h__ = .01f;
bul = FALSE_;
ih = 0;
l6:
de46r_c((U_fp)f2_c, &m, &jstart, &hmin, &hmax, &eps, &p, yx, dyx, &x,
&h__, &bul, z__, delty, df, rfn, rf, yp, dyp, zp, &ierr);
printf("\n %5i \n", ierr);
if (ierr != 0) {
goto l20;
}
++ih;
y1 = (float)sin(x) + (float)cos(x) + x * .5f;
y2 = -(float)sin(x) - (float)cos(x) + x * .5f;
dy1 = (float)cos(x) - (float)sin(x) + .5f;
dy2 = -(float)cos(x) + (float)sin(x) + .5f;
printf("\n %16.7e \n", x);
printf("\n %16.7e %16.7e \n", yx[0], yx[1]);
printf("\n %16.7e %16.7e \n", y1, y2);
printf("\n %16.7e \n", h__);
printf("\n %16.7e %16.7e \n", dyx[0], dyx[1]);
printf("\n %16.7e %16.7e \n", dy1, dy2);
switch (ih) {
case 1: goto l6;
case 2: goto l6;
case 3: goto l102;
case 4: goto l103;
case 5: goto l104;
case 6: goto l6;
case 7: goto l20;
}
l102:
jstart = -1;
h__ = -.005f;
goto l6;
l103:
jstart = -1;
h__ = .02f;
goto l6;
l104:
jstart = -1;
h__ = .01f;
goto l6;
l20:
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 */
Результаты:
после первого обращения к подпрограмме -
x yx(1) yx(2)
1.000000000001 - 02 1.014949833751 + 00 - 1.004949833752 + 00
y1 y2
1.014949833745 + 00 - 1.004949833748 + 00
h__ dyx(1) dyx(2)
1.000000000001 - 02 1.489950167079 + 00 - 4.899501670834 - 01
dy1 dy2
1.489950167079 + 00 - 4.899501670798-01
после второго обращения к подпрограмме -
x yx(1) yx(2)
2.000000000001 - 02 1.029798673357 + 00 - 1.009798673364 + 00
y1 y2
1.029798673355 + 00 - 1.009798673360 + 00
h__ dyx(1) dyx(2)
2.000000000001 - 02 1.479801339963 + 00 - 4.798013399754 - 01
dy1 dy2
1.479801339969 + 00 - 4.798013399704-01
после третьего обращения к подпрограмме -
x yx(1) yx(2)
4.000000000002 - 02 1.059189440843 + 00 - 1.019189440856 + 00
y1 y2
1.059189440844 + 00 - 1.019189440849 + 00
h__ dyx(1) dyx(2)
4.000000000002 - 02 1.459210772458 + 00 - 4.592107724779-01
dy1 dy2
1.459210772471 + 00 - 4.592107724729-01
после четвертого обращения к подпрограмме -
x yx(1) yx(2)
1.500000000000 - 02 1.022386939610 + 00 - 1.007386939622 + 00
y1 y2
1.022386939612 + 00 - 1.007386939615 + 00
h__ dyx(1) dyx(2)
- 1.000000000001 - 02 1.484888064588 + 00 - 4.848880646064 - 01
dy1 dy2
1.484888064599 + 00 - 4.848880646005 - 01
после пятого обращения к подпрограмме -
x yx(1) yx(2)
4.000000000002 - 02 1.059189440841 + 00 - 1.019189440856 + 00
y1 y2
1.059189440844 + 00 - 1.019189440849 + 00
h__ dyx(1) dyx(2)
4.000000000002 - 02 1.459210772457 + 00 - 4.592107724784 - 01
dy1 dy2
1.459210772471 + 00 - 4.592107724729 - 01
после шестого обращения к подпрограмме -
x yx(1) yx(2)
3.000000000000 - 02 1.044545533945 + 00 - 1.014545533959 + 00
y1 y2
1.044545533947 + 00 - 1.014545533950 + 00
h__ dyx(1) dyx(2)
2.000000000001 - 02 1.469554533529 + 00 - 4.695545335508 - 01
dy1 dy2
1.469554533543 + 00 - 4.695545335444 - 01
после седьмого обращения к подпрограмме -
x yx(1) yx(2)
5.000000000001 - 02 1.073729429656 + 00 - 1.023729429677 + 00
y1 y2
1.073729429661 + 00 - 1.023729429664 + 01
h__ dyx(1) dyx(2)
4.000000000002 - 02 1.448771091098 + 00 - 4.487710911303 - 01
dy1 dy2
1.448771091120 + 00 - 4.487710911217 - 01