|
Текст подпрограммы и версий de15r_c.zip , de15d_c.zip |
Тексты тестовых примеров tde15r_c.zip , tde15d_c.zip |
Выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом Рунге - Кутта - Фельберга.
Выполняется один шаг численного интегрирования системы М обыкновенных дифференциальных уравнений первого порядка
Y ' = F (X, Y) ,
Y = ( y1, ... , yM ) ,
F = ( f1 (X, y1, ... , yM), ... , fM (X, y1, ... , yM) )
методом Рунге - Кутта - Фельберга пятого порядка точности. По заданному значению решения YX в узле Xn вычисляется значение этого решения в узле Xn + H. Bсе компоненты решения вычисляются с контролем точности по меpе погрешности, который заключается в следующем.
Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы Р, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. B качестве абсолютной погрешности решения принимается оценка главного члена асимптотического разложения погрешности метода на шаге, получаемая при вычитании двух выражений, представляющих приближенные значения решения пятого и четвертого порядка точности. При этом на одном шаге интегрирования для определения решения и погрешности используются всего шесть вычислений правой части системы. Значение H может быть меньше или равно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.
Дж.Форсайт, М.Малькольм, К.Моулер, Машинные методы математических вычислений, "Мир", M., 1980.
int de15r_c (S_fp f, integer *m, integer *jstart, real *hmin,
real *eps, real *p, real *yx, real *x, real *h, logical *bul,
real *xp, real *yp, real *dy, real *r1, real *r2, real *r3, real *r4,
integer *ierr)
Параметры
| f - |
имя подпрограммы вычисления значений правой
части дифференциального уравнения. Первый
оператоp подпрограммы должен иметь вид: int f (float *x, float *y, float *dy, int *m). Здесь: x, y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в dy. B случае системы уравнений, т.е. когда m ≠ 1 , параметры y и dy представляют массивы длины m (тип параметров x, y и dy: вещественный); |
| m - | количество уравнений в системе (тип: целый); |
| jstart - | целый указатель режима использования подпрограммы, имеющий следующие значения: |
| 0,+1 - | выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных и шага интегрирования, заданных параметрами x, yx и h, соответственно; |
| -1 - | повторить последний шаг интегрирования с новыми значениями паpаметpов h и/или hmin; |
| на выходе из подпрограммы jstart полагается равным + 1; | |
| hmin - | минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| eps - | допустимая меpа погрешности, с которой тpебуется вычислить все компоненты решения (тип: вещественный); |
| p - | граница перехода, используемая при оценке меры погрешности решения (тип: вещественный); |
| yx, x - | заданные вещественные значения решения и соответствующее ему значение аргумента; в pезультате работы подпрограммы в x получается новое значение аргумента, а в yx - соответствующее значение решения; в случае системы уравнений, т.е. когда m ≠ 1, yx задается одномерным массивом длины m; |
| h - | вещественная переменная, содержащая значение шага интегрирования; если для этого значения шага точность приближенного решения достигается, то именно он и реализуется подпрограммой, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность eps; на выходе из подпрограммы h содержит рекомендуемое подпрограммой значение следущего шага интегрирования, определяемое ею с целью достижения более экономного способа интегрирования; |
| bul - | логическая переменная, значение которой при обращении к подпрограмме полагается равным TRUE_, если заданный в h шаг выводит в конец интервала интегрирования, и FALSE_ в противном случае; в результате работы подпрограммы bul pавно FALSE_, если вместо исходного шага интегрирования был реализован меньший шаг; в противном случае, т.е. когда был выполнен именно заданный при обращении в h шаг, значение параметра bul не меняется; |
| xp, yp - | вещественные рабочая переменная и одномерный рабочий массив длины m; значения параметpов xp, yp на выходе из подпрограммы pавны тем значениям, которые имели параметры x, yx, соответственно, при входе в нее (т.е. предыдущий узел и решение в нем); |
|
dy - r1, r2 r3, r4 | одномерные вещественные рабочие массивы длины m; |
| ierr - | целая переменная, значение которой в pезультате работы подпрограммы полагается равным 65, если какая - нибудь компонента pешения не может быть вычислена с требуемой точностью eps; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров h, hmin и значением jstart = - 1. |
Версии
| de15d_c - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений первого порядка методом Рунге - Кутта - Фельберга с удвоенным числом значащих цифр. При этом параметры hmin, eps, p, yx, x, h, xp, yp, dy, r1, r2, r3, r4 и параметры x, y и dy в подпрограмме f должны иметь тип double. |
Вызываемые подпрограммы
| utde16_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de15r_c. |
| utde17_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de15d_c. |
Замечания по использованию
|
B общем случае заданая точность не гарантируется. При работе подпрограммы значения параметров m, hmin, eps, p сохраняются. При работе подпрограммы счета правой части f значения параметров x, y и m не должны изменяться. При обращении к подпрограмме со значением jstart = - 1 в качестве исходных значений аргумента и решения принимаются значения параметров xp и yp, соответственно, т.е. те значения, которые эти параметры получили после самого последнего обращения к подпрограмме с неотрицательным значением jstart. |
y1' = -y2 - 0.1 x - 0.9
y2' = y1 - 0.1 x - 1.1
y1 = 1 , y2 (0) = -2 , 0 ≤ x ≤ π
Точное решение системы:
y1 = 0.1 x + sin x + 1. , y2 = -0.1 x - cos x - 1.
Приводятся подпрограмма вычисления значений правой части системы и фрагмент вызывающей программы, а также результаты счета.
int main(void)
{
/* Builtin functions */
double sin(double), cos(double);
/* Local variables */
extern int de15r_c(U_fp, int *, int *, float *, float *, float *,
float *, float *, float *, logical *, float *,
float *, float *, float *, float *, float *,
float *, int *);
static float hmin;
static int ierr;
static float h__;
static int m;
extern int ffelb_c();
static float p, x, r1[2], r2[2], r3[2], r4[2], y1, y2;
static int ih;
static float dy[2], xp, yp[2], yx[2];
static int jstart;
static logical bul;
static float eps;
m = 2;
jstart = 0;
hmin = 1e-12f;
eps = 1e-5f;
p = 100.f;
yx[0] = 1.f;
yx[1] = -2.f;
x = 0.f;
h__ = .01f;
bul = FALSE_;
ih = 0;
L6:
de15r_c((U_fp)ffelb_c, &m, &jstart, &hmin, &eps, &p, yx, &x, &h__, &bul,
&xp, yp, dy, r1, r2, r3, r4, &ierr);
printf("\n %5i \n", ierr);
if (ierr != 0) {
goto L20;
}
y1 = x * .1f + (float)sin(x) + 1.f;
y2 = x * -.1f - (float)cos(x) - 1.f;
printf("\n %16.7e \n", x);
printf("\n %16.7e %16.7e \n", yx[0], yx[1]);
printf("\n %16.7e \n", h__);
printf("\n %16.7e %16.7e \n", y1, y2);
++ih;
switch (ih) {
case 1: goto L6;
case 2: goto L6;
case 3: goto L103;
case 4: goto L104;
case 5: goto L105;
case 6: goto L6;
case 7: goto L20;
}
L103:
jstart = -1;
h__ = -.005f;
goto L6;
L104:
jstart = -1;
h__ = .02f;
goto L6;
L105:
jstart = -1;
h__ = .01f;
goto L6;
L20:
return 0;
} /* main */
int ffelb_c(float *x, float *y, float *z, int *m)
{
static float r__;
/* Parameter adjustments */
--z__;
--y;
/* Function Body */
r__ = *x * .1f;
z__[1] = -y[2] - r__ - .9f;
z__[2] = y[1] - r__ - 1.1f;
return 0;
} /* ffelb_c */
Результаты:
после первого обращения к подпрограмме -
x yx (1) yx (2)
9.999999999991-03 1.010999833334 + 00 -2.000950000420 + 00
h__ y1 y2
1.999999999998-02 1.010999833334 + 00 -2.000950000416 + 00
после второго обращения к подпрограмме -
x yx (1) yx (2)
2.999999999997-02 1.032995500202 + 00 -2.002550033754 + 00
h__ y1 y2
3.999999999996-02 1.032995500202 + 00 -2.002550033747 + 00
после третьего обращения к подпрограмме -
x yx (1) yx (2)
6.999999999994-02 1.076942847336 + 00 -2.004551000264 + 00
h__ y1 y2
7.999999999993-02 1.076942847336 + 00 -2.004551000253 + 00
после четвертого обращения к подпрограмме -
x yx (1) yx (2)
2.499999999998-02 1.027497395913 + 00 -2.002187516282 + 00
h__ y1 y2
-9.999999999991-03 1.027497395913 + 00 -2.002187516275 + 00
после пятого обращения к подпрограмме -
x yx (1) yx (2)
4.999999999995-02 1.054979169268 + 00 -2.003750260403 + 00
h__ y1 y2
3.999999999996-02 1.054979169270 + 00 -2.003750260392 + 00
после шестого обращения к подпрограмме -
x yx (1) yx (2)
3.999999999996-02 1.043989334185 + 00 -2.003200106668 + 00
h__ y1 y2
1.999999999998-02 1.043989334186 + 00 -2.003200106661 + 00
после седьмого обращения к подпрограмме -
x yx (1) yx (2)
5.999999999995-02 1.065964006477 + 00 -2.004200539945 + 00
h__ y1 y2
3.999999999996-02 1.065964006479 + 00 -2.004200539934 + 00