|
Текст подпрограммы и версий de37r_c.zip , de37d_c.zip |
Тексты тестовых примеров tde37r_c.zip , tde37d_c.zip |
Вычисление решения задачи Коши для жесткой линейной системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования неявным методом Рунге - Кутта.
Решается задача Коши для жесткой линейной системы M обыкновенных дифференциальных уравнений первого порядка с переменными коэффициентами
Y ' (X) = A(X) * Y(X) + φ(X) ,
Y = ( y1, ..., yM ) ,
A(X) = ( ai j(X) ), i, j = 1, ..., M ,
φ(X) = ( φ1(X), ..., φM(X) ) .
с начальными условиями, заданными в точке XN:
Y (XN) = YN, YN = ( y10, ..., yM0 ) ,
трехстадийным А - устойчивым неявным методом Рунге - Кутта шестого порядка точности. Решение вычисляется в одной точке ХК, которая является концом интервала интегрирования.
Bсе компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Абсолютная погрешность приближенного решения оценивается по правилу Рунге.
Современные численные методы решения обыкновенных дифференциальных уравнений. Ред. Дж.Холл и Дж.Уатт. "Мир", M., 1979.
Butcher J.C. Implicit Runge - Kutta processes. Math. Comp., 18, 50 - 64, 1964.
int de37r_c (real *fa, real *fi, integer *m, real *xn,
real *yn, real *xk, real *hmin, real *eps, real *p, real *h,
real *y, real *rab, integer *ir, integer *ierr)
Параметры
| fa - |
подпрограмма вычисления матрицы системы a (x) в
любой точке x. Первый оператор подпрограммы
должен иметь вид: int fa (float *a, float *x, int *m) Здесь: a - двумерный массив размера m * m, в который помещается матрица системы, вычисленная при значении аргумента x (тип параметров a, x: вещественный); |
| fi - |
подпрограмма вычисления неоднородности правой
части системы φ (x)
в любой точке x. Первый
оператор подпрограммы должен иметь вид: int fi (float *g, float *x, int *m) Здесь g - одномерный массив длины m, в который помещается неоднородность правой части системы, вычисленная при значении аргумента x (тип параметров g, x: вещественный); |
| m - | количество уравнений в системе (тип: целый); |
| xn, yn - | начальные значения аргумента и решения; в случае системы уравнений (т.е. m ≠ 1) yn представляет одномерный массив длины m (тип: вещественный); |
| xk - | значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования); xk может быть больше, меньше или pавно xn (тип: вещественный); |
| hmin - | минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| eps - | допустимая меpа погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
| p - | граница перехода, используемая при оценке погрешности решения (тип: вещественный); |
| h - | вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если xn < xk, отрицательным, если XN > XK, или без всякого учета в виде абсолютной величины; |
| y - | искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента xk; для системы уравнений (когда m ≠ 1) задается одномерным массивом длины m. В случае совпадения значений параметров xn и xk значение y полагается равным начальному значению yn (тип: вещественный); |
| rab - | одномерный рабочий массив вещественного типа длины (10*m*m + 9*m + 1); |
| ir - | целый одномерный рабочий массив длины 3*m; |
| ierr - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая-нибудь компонента решения не может быть вычислена с требуемой точностью eps; в этом случае интегрирование системы прекращается; при желании интегрирование системы можно повторить обращением к подпрограмме с новыми значениями параметров hmin и h. |
Версии
| de37d_c - | вычисление решения задачи Коши для жесткой линейной системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования неявным методом Рунге - Кутта с удвоенным числом значащих цифр. При этом параметры xn, yn, xk, hmin, eps, p, h, y, rab и параметры a, g, x в подпрограммах fa и fi должны иметь тип double. |
Вызываемые подпрограммы
|
de36r_c - de36d_c | выполнение одного шага численного интегрирования жесткой линейной системы обыкновенных дифференциальных уравнений первого порядка неявным методом Рунге - Kутта. |
|
utde16_c - utde17_c | подпрограммы выдачи диагностических сообщений. |
| Подпрограммы de36r_c, utde16_c вызываются при работе подпрограммы de37r_c, а подпрограммы de36d_c, utde17_c - при pаботе de37d_c. |
Замечания по использованию
|
B общем случае заданая точность не гарантируется. При работе подпрограммы значения параметров m, xn, yn, xk, hmin, eps, p сохраняются. При работе подпрограмм fa и fi значения параметров x и m не должны изменяться. Если после работы подпрограммы нет необходимости иметь начальное значение решения yn, то параметры yn и y при обращении к ней можно совместить. При этом следует иметь в виду, что в случае аварийного выхода из подпрограммы, т.е. со значением ierr = 65, значение параметра yn будет испорчено. Tак как при интегрировании уравнений с помощью подпрограмм de37r_c и de37d_c используются (для версии фортран) общие блоки с именами com36r и com36d, соответственно, то пользователю не рекомендуется использовать для своих целей (для версии фортран) общие блоки с указанными именами. |
y1' = - 20y1 + y2 , y1 = 2 ,
y2' = 19y1 - 2y2 , y2 = 18 , 0 ≤ x ≤ 5 .
Точное решение системы:
y1 = e - x + e - 21x ,
y2 = 19e - x - e - 21x .
int main(void)
{
/* Builtin functions */
double exp(double);
/* Local variables */
extern int de37r_c(U_fp, U_fp, int *, float *, float *, float *,
float *, float *, float *, float *, float *,
float *, int *, int *);
static float hmin;
extern int fbut_c();
static int ierr;
static float h__;
static int k, m;
static float p, x, y[2], y1, y2;
static int ir[6];
static float xk, xn, yn[2];
extern int fbutfi_c();
static float rab[59], eps;
m = 2;
xn = 0.f;
yn[0] = 2.f;
yn[1] = 18.f;
xk = 5.f;
hmin = 1e-5f;
eps = 1e-5f;
p = 100.f;
for (k = 1; k <= 2; ++k) {
h__ = .01f;
de37r_c((U_fp)fbut_c, (U_fp)fbutfi_c, &m, &xn, yn, &xk, &hmin, &eps,
&p, &h__, y, rab, ir, &ierr);
printf("\n %5i \n", ierr);
x = xk;
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", xk);
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);
eps = 5e-6f;
/* l200: */
}
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 */
Результаты:
после первого обращения к подпрограмме -
xk y(1) y(2)
5.000000000000 + 00 6.737934370626 - 03 1.280207530417 - 01
h__ y1 y2
1.280000000001 + 00 6.737946999117 - 03 1.280209929830 - 01
после второго обращения к подпрограмме -
xk y(1) y(2)
5.000000000000 + 00 6.737946721380 - 03 1.280209877059 - 01
h__ y1 y2
6.400000000003 - 01 6.737946999117 - 03 1.280209929830 - 01