|
Текст подпрограммы и версий de31r_c.zip , de31d_c.zip |
Тексты тестовых примеров tde31r_c.zip , tde31d_c.zip |
Вычисление решения задачи Коши для жестких линейных однородных систем обыкновенных дифференциальных уравнений первого порядка с переменными коэффициентами экспоненциальным методом.
Решается задача Коши для жесткой линейной однородной системы М обыкновенных дифференциальных уравнений первого порядка с переменными коэффициентами
(1) Y ' = A(x) Y
здесь Y = ( y1, ..., yM ), А(x) - квадратная вещественная матрица порядка M, элементами которой являются трижды непрерывно дифференцируемые функции. Начальные значения заданы в точке XN :
Y(XN) = YN , YN = ( y10,..., yM0 ) ,
Предполагается, что матрица A (x) имеет все собственные числа с отрицательными вещественными частями, среди которых встречаются большие по абсолютной величине. Решение вычисляется в одной точке XK (которая является концом интервала интегрирования).
Интервал интегрирования разбивается на NS равных частей длины T, на каждой из которых исходная система уравнений (1) аппроксимируется системой уравнений с постоянными коэффициентами
~ ~ ~ ~ (2) Y ' = A Y , A = CONST.
Решение системы (2) вычисляется с помощью матричной экспоненты
~
(3) EXP( A T ) =
~ N
= [ EXP( A H ) ] , H = T/N , N > 0
Значения N и NS задаются пользователем при обращении к подпрограмме, причем значение N выбирается таким образом, чтобы выполнялось условие
(4) max || A(x) H || ≤ 1
XN ≤ x ≤ XK
(в качестве || A(x) H || можно взять максимальную сумму модулей элементов матрицы A (x) H по стpокам).
С.Ф.Залеткин, O численном решении задачи Коши для обыкновенных линейных однородных дифференциальных уравнений на больших отрезках интегрирования, "Вычислительные методы и программирование", вып.26, Изд-во МГУ, 1977.
int de31r_c (S_fp fa, integer *m, real *xn, real *yn, real *xk,
integer *n, integer *ns, real *y, real *a, real *e, real *r,
real *r1, real *r2, real *r3, integer *ierr)
Параметры
| fa - |
подпрограмма вычисления матрицы A(x) в
любой точке x интервала интегрирования.
Первый оператор подпрограммы должен иметь вид: int fa (float *a, float *x, int *m). Здесь a - двумерный массив размера m*m, в который помещается матрица системы, вычисленная при значении аргумента x, т.е. A (x) (тип параметров a и x: вещественный); |
| m - | количество уравнений в системе (тип: целый); |
| xn, yn - | начальные значения аргумента и решения; yn представляет одномерный массив длины m (тип: вещественный); |
| xk - | значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования). xk может быть больше, меньше или pавно xn (тип: вещественный); |
| n - |
целое число равных частей, на которые
делится отрезок интегрирования системы
линейных дифференциальных уравнений с
постоянными коэффициентами; при обращении к
подпрограмме параметр n выбирается таким образом, чтобы || a(x)*( ((xk - xn)/ns)/n )|| ≤ 1 для xn ≤ x ≤ xk (см. замечания по использованию); |
| ns - | число равных частей ,на которые разбивается интервал интегрирования. ns должно быть таким, чтобы система уравнений (2) достаточно хорошо аппроксимировала исходную систему уравнений с переменными коэффициентами; ns > 0 (тип: целый); |
| y - | искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента xk; y представляет одномерный массив длины m. В случае совпадения значений параметров xn и xk значение y полагается равным начальному значению yn (тип: вещественный); |
|
a, e - r1, r2 | вещественные двумерные рабочие массивы размеpа m*m; |
| r, r3 - | вещественные одномерные рабочие массивы длины m; |
| ierr - | целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом: |
| ierr=65 - | когда значение параметра n меньше 1; |
| ierr=66 - | когда значение параметра ns меньше 1. |
| В каждом из этих случаев интегрирование системы прекращается. |
Версии
| de31d_c - | вычисление решения задачи Коши для линейной однородной системы обыкновенных дифференциальных уравнений первого порядка с переменными коэффициентами экспоненциальным методом с повышенной точностью. При этом параметры xn, yn, xk, y, a, e, r, r1, r2, r3, а также параметры a, x в подпрограмме fa должны иметь тип double. |
Вызываемые подпрограммы
| utde10_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de31r_c. |
| utde11_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de31d_c. |
Замечания по использованию
|
Значение параметра n, задаваемое при обращении к
подпрограмме таким образом, чтобы |
Использование подпрограммы иллюстрируется на примере:
y1' = - (2 + x) y1 / (1 + x) + 20 x y2 ,
y2' = -20 x y1 - (2 + x) y2 / (1 + x) ,
y1(0) = 0 , y2(0) = 1 , 0 ≤ x ≤ 6 .
Приводятся подпрограмма вычисления матрицы системы и фрагмент вызывающей программы, а также результаты счета.
int main(void)
{
/* Local variables */
extern int de31r_c(U_fp, int *, float *, float *, float *, int *, int *,
float *, float *, float *, float *, float *, float *,
float *, int *);
static int ierr;
static float a[4] /* was [2][2] */, e[4] /* was [2][2] */;
static int m, n;
static float r__[2], y[2], r1[4] /* was [2][2] */,
r2[4] /* was [2][2] */, r3[2];
extern int fa_c();
static float xk, xn, yn[2];
static int is1;
m = 2;
xn = 0.f;
yn[0] = 0.f;
yn[1] = 1.f;
xk = 6.f;
n = 128;
is1 = 128;
de31r_c((U_fp)fa_c, &m, &xn, yn, &xk, &n, &is1, y, a, e, r__, r1, r2, r3,
&ierr);
printf("\n %15.7e \n", xk);
printf("\n %15.7e %15.7e \n", y[0], y[1]);
return 0;
} /* main */
int fa_c(float *a, float *x, int *m)
{
/* System generated locals */
int a_dim1, a_offset;
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
/* Parameter adjustments */
a_dim1 = *m;
a_offset = 1 + a_dim1 * 1;
a -= a_offset;
/* Function Body */
a_ref(1, 1) = -(*x + 2.f) / (*x + 1.f);
a_ref(*m, *m) = a_ref(1, 1);
a_ref(*m, 1) = *x * -20.f;
a_ref(1, *m) = -a_ref(*m, 1);
return 0;
} /* fa_c */
Результаты:
y(1) = 0.3395896401 * 10-3
y(2) = -0.1004661347 * 10-3
ierr = 0