|
Текст подпрограммы и версий de36r_c.zip , de36d_c.zip |
Тексты тестовых примеров tde36r_c.zip , tde36d_c.zip |
Выполнение одного шага интегрирования жесткой линейной системы обыкновенных дифференциальных уравнений первого порядка неявным методом Рунге - Кутта.
Выполняется один шаг численного интегрирования жесткой линейной системы M обыкновенных дифференциальных уравнений первого порядка с переменными коэффициентами
Y ' (X) = A(X) * Y(X) + φ(X) ,
Y = ( y1, ..., yM ) ,
A(X) = ( a i j(X) ) , i, j = 1, ..., M ,
φ(X) = ( φ1(X), ..., φM(X) ) .
По заданному значению решения YX в узле Xn вычисляется значение этого решения в узле Xn + H. Вычисление производится по неявной трехчленной формуле Рунге - Кутта шестого порядка точности. Значение H может быть меньше или pавно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.
Bсе компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Абсолютная погрешность приближенного решения оценивается по правилу Рунге.
Современные численные методы решения обыкновенных дифференциальных уравнений. Ред. Дж.Холл и Дж.Уатт. "Мир", M., 1979.
Butcher J.C. Implicit Runge - Kutta processes. Math. Comp., 18, 50 - 64, 1964.
int de36r_c (S_fp fa, S_fp fi, integer *m, integer *jstart,
real *hmin, real *eps, real *p, real *yx, real *x, real *h,
logical *bul, real *xp, real *yp, integer *ir, real *r1, real *r2,
real *r3, real *r4, real *r5, integer *ierr)
Параметры
| fa - | подпрограмма вычисления матрицы системы a (x) в любой точке x. Первый оператор подпрограммы должен иметь вид: |
|
int fa (float *a, float *x, int *m) Здесь: a - двумерный массив размера m * m, в котоpом помещается матрица системы, вычисленная при значении аргумента x (тип параметров a, x: вещественный); |
| fi - | подпрограмма вычисления неоднородности правой части системы φ (x) в любой точке x. Первый оператор подпрограммы должен иметь вид: |
|
int fi (float *g, float *x, int *m) Здесь g - одномерный массив длины m, в который помещается неоднородность правой части системы, вычисленная при значении аргумента x (тип параметров g, x: вещественный); |
| m - | количество уравнений в системе (тип: целый); |
| jstart - | целый указатель режима использования подпрограммы, имеющий следующие значения: |
| 0 - | первое обращение к подпрограмме должно быть выполнено с нулевым значением jstart; |
| +1 - | выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных и шага интегрирования, заданных параметрами x, yx и h, соответственно; |
| -1 - | повторить последний шаг интегрирования с новыми значениями параметров h и/или hmin; |
| на выходе из подпрограммы jstart полагается равным + 1; |
| hmin - | минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| eps - | допустимая меpа погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
| p - | граница перехода, используемая при оценке меры погрешности решения (тип: вещественный); |
| yx, x - | заданные вещественные значения решения и соответствующее ему значение аргумента; в результате работы подпрограммы в x получается новое значение аргумента, а в yx - соответствующее значение решения; в случае системы уравнений, т.е. когда m ≠ 1, yx задается одномерным массивом длины m; |
| h - | вещественная переменная, содержащая значение шага интегрирования; если для этого значения шага точность приближенного решения достигается, то именно он и реализуется подпрограммой, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность eps; на выходе из подпрограммы h содержит рекомендуемое подпрограммой значение следущего шага интегрирования, определяемое ею с целью достижения более экономного способа интегрирования; |
| bul - | логическая переменная, значение которой при обращении к подпрограмме полагается равным TRUE_, если заданный в h шаг выводит в конец интервала интегрирования, и FALSE_ в противном случае; в результате работы подпрограммы bul pавно FALSE_, если вместо исходного шага интегрирования был реализован меньший шаг; в противном случае, т.е. когда был выполнен именно заданный при обращении в h шаг, значение параметра bul не меняется; |
| xp, yp - | вещественная рабочая переменная и одномерный рабочий массив длины m, соответственно; значения паметров xp, yp на выходе из подпрограммы равны тем значениям, которые имели параметры x, yx при входе в нее (т.е. предыдущий узел и решение в нем); |
| ir - | целый одномерный рабочий массив длины 3 * m; |
| r1, r2 - | вещественные одномерные рабочие массивы длины ( 3m * (3m + 1) ) и 3m, соответственно; |
| r3 - | вещественный двумерный рабочий массив размера m * m; |
| r4, r5 - | вещественные одномерные рабочие массивы длины m; |
| ierr - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью eps; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров h, hmin и значением jstart = -1. |
Версии
| de36d_c - | выполнение одного шага численного интегрирования жесткой линейной системы обыкновенных дифференциальных уравнений первого порядка неявным методом Рунге - Кутта с удвоенным числом значащих цифр. При этом параметры hmin, eps, p, yx, x, h, xp, yp, r1, r2, r3, r4, r5 и параметры a, g, x в подпрограммах fa и fi должны иметь тип double. |
Вызываемые подпрограммы
| utde16_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de36r_c. |
| utde17_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de36d_c. |
Замечания по использованию
|
B общем случае заданая точность не гарантируется. При работе подпрограммы значения параметров m, hmin, eps, p сохраняются. При работе подпрограмм fa и fi значения параметров x и m не должны изменяться. При обращении к подпрограмме со значением jstart = -1 в качестве исходных значений аргумента и решения принимаются значения параметров xp и yp, соответственно, т.е. те значения, которые эти параметры получили после самого последнего обращения к подпрограмме с неотрицательным значением jstart. |
y1' = - 20y1 + y2 ,
y2' = 19y1 - 2y2 ,
Точное решение системы:
y1 = e - x + e - 21x ,
y2 = 19e - x - e - 21x .
int main(void)
{
/* Builtin functions */
double exp(double);
/* Local variables */
extern int de36r_c(U_fp, U_fp, int *, int *, float *, float *, float *,
float *, float *, float *, logical *, float *,
float *, int *, float *, float *, float *, float *,
float *, int *);
static float hmin;
extern int fbut_c();
static int ierr;
static float h__;
static int m;
static float p, x, r1[42], r2[6], r3[4], r4[2], r5[2], y1, y2;
static int ih, ir[6];
static float xp, yp[2];
extern int fbutfi_c();
static float yx[2];
static int jstart;
static logical bul;
static float eps;
ih = 0;
m = 2;
x = 0.f;
yx[0] = 2.f;
yx[1] = 18.f;
hmin = 1e-6f;
eps = 1e-4f;
p = 100.f;
jstart = 0;
h__ = .01f;
bul = FALSE_;
L6:
de36r_c((U_fp)fbut_c, (U_fp)fbutfi_c, &m, &jstart, &hmin, &eps, &p, yx,
&x, &h__, &bul, &xp, yp, ir, r1, r2, r3, r4, r5, &ierr);
printf("\n %5i \n", ierr);
++ih;
y1 = (float)exp((float)(-x)) + (float)exp((float)(x * -21.f));
y2 = (float)exp((float)(-x)) * 19.f - (float)exp((float)(x * -21.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);
switch (ih) {
case 1: goto L6;
case 2: goto L6;
case 3: goto L104;
case 4: goto L105;
case 5: goto L106;
case 6: goto L6;
case 7: goto L20;
}
L104:
jstart = -1;
h__ = .005f;
goto L6;
L105:
jstart = -1;
h__ = .02f;
goto L6;
L106:
jstart = -1;
h__ = .03f;
goto L6;
L20:
return 0;
} /* main */
int fbut_c(float *a, float *x, int *m)
{
#define a_ref(a_1,a_2) a[(a_2)*2 + a_1]
/* Parameter adjustments */
a -= 3;
/* Function Body */
a_ref(1, 1) = -20.f;
a_ref(1, 2) = 1.f;
a_ref(2, 1) = 19.f;
a_ref(2, 2) = -2.f;
return 0;
} /* fbut_c */
int fbutfi_c(float *fi, float *x, int *m)
{
/* Parameter adjustments */
--fi;
/* Function Body */
fi[1] = 0.f;
fi[2] = 0.f;
return 0;
} /* fbutfi_c */
Результаты:
после первого обращения к подпрограмме -
x yx(1) yx(2)
9.999999999991 - 03 1.800634079711 + 00 1.800036259522 + 01
h__ y1 y2
1.999999999998 - 02 1.800634079724 + 00 1.800036259528 + 01
после второго обращения к подпрограмме -
x yx(1) yx(2)
2.999999999997 - 02 1.503037334354 + 00 1.790587333651 + 01
h__ y1 y2
3.999999999996 - 02 1.503037334556 + 00 1.790587333636 + 01
после третьего обращения к подпрограмме -
x yx(1) yx(2)
6.999999999994 - 02 1.162319294410 + 00 1.748555710353 + 01
h__ y1 y2
7.999999999993 - 02 1.162319305093 + 00 1.748555709302 + 01
после четвертого обращения к подпрограмме -
x yx(1) yx(2)
3.499999999997 - 02 1.445110875044 + 00 1.786699744989 + 01
h__ y1 y2
9.999999999991 - 03 1.445110875235 + 00 1.786699744992 + 01
после пятого обращения к подпрограмме -
x yx(1) yx(2)
4.999999999995 - 02 1.301167173347 + 00 1.772342131648 + 01
h__ y1 y2
3.999999999996 - 02 1.301167173613 + 00 1.772342131642 + 01
после шестого обращения к подпрограмме -
x yx(1) yx(2)
5.999999999995 - 02 1.225418558231 + 00 1.760987211324 + 01
h__ y1 y2
5.999999999995 - 02 1.225418560085 + 00 1.760987211161 + 01
после седьмого обращения к подпрограмме -
x yx(1) yx(2)
1.199999999999 - 01 9.673799790717 - 01 1.677102875494 + 01
h__ y1 y2
1.199999999999 - 01 9.673800434666 - 01 1.677102869086 + 01