|
Текст подпрограммы и версий de30r_c.zip , de30d_c.zip |
Тексты тестовых примеров tde30r_c.zip , tde30d_c.zip |
Вычисление решения задачи Коши для жестких линейных однородных систем обыкновенных дифференциальных уравнений первого порядка с постоянными коэффициентами экспоненциальным методом.
Решается задача Коши для жесткой линейной однородной системы М обыкновенных дифференциальных уравнений первого порядка с постоянными коэффициентами
Y ' = AY ,
здесь Y = ( y1, ..., yM ), А - квадратная вещественная матрица порядка M. Начальные значения заданы в точке XN :
Y(XN) = YN , YN = ( y10,..., yM0 ) ,
Предполагается, что матрица А имеет все собственные числа с отрицательными вещественными частями, среди которых встречаются большие по абсолютной величине. Решение вычисляется в одной точке XK (которая является концом интервала интегрирования) по следующей формуле:
(1) Y(XK) = [EXP( AH )]N*YN , H = ( XK-XN ) / N , N > 0.
Для аппроксимации матричной экспоненты EXP(AH) используется частичная сумма ряда
7
EXP(AH) = ∑ ( AH )k / k! ,
k=0
а значение N выбирается таким, чтобы выполнялось условие
(2) || A H || ≤ 1
(в качестве || A H || можно взять максимальную сумму модулей элементов матрицы AH по стpокам). Это значение N задается пользователем при обращении к подпрограмме.
С.Ф.Залеткин, O численном решении задачи Коши для обыкновенных линейных однородных дифференциальных уравнений на больших отрезках интегрирования, "Вычислительные методы и программирование", вып. 26, Изд - во МГУ, 1977.
int de30r_c (integer *m, real *xn, real *yn, real *xk,
integer *n, real *a, real *y, real *e, real *r, real *r1,
real *r2, integer *ierr)
Параметры
| m - | количество уравнений в системе (тип: целый); |
| xn, yn - | начальные значения аргумента и решения; yn представляет одномерный массив длины m (тип xn и yn: вещественный); |
| xk - | значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования). xk может быть больше, меньше или pавно xn (тип: вещественный); |
| n - | целое число равных частей, на которые делится отрезок интегрирования; при обращении к подпрограмме параметр n выбирается таким образом, чтобы || A *(xk - xh) / n || ≤ 1 (см. замечания по использованию); |
| a - | вещественный двумерный массив размера m*m, содержащий матрицу системы уравнений; |
| y - | искомое решение задачи Коши, вычисленное подпрограммой при значении аргумента xk; y представляет одномерный массив длины m. B случае совпадения значений параметров xn и xk значение y полагается равным начальному значению yn (тип: вещественный); |
|
e - r1, r2 | вещественные двумерные рабочие массивы размеpа m*m; |
| r - | вещественный одномерный рабочий массив длины m; |
| ierr - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если значение параметра n меньше 1. B этом случае интегрирование системы прекращается; при этом значения параметров m, xn, yn, xk, n, a не изменяются. |
Версии
| de30d_c - | вычисление решения задачи Коши для линейной однородной системы обыкновенных дифференциальных уравнений первого порядка с постоянными коэффициентами экспоненциальным методом с повышенной точностью. При этом параметры xn, yn, xk, a, y, e, r1, r2, r должны иметь тип double. |
Вызываемые подпрограммы
| utde10_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de30r_c. |
| utde11_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de30d_c. |
Замечания по использованию
|
Хотя подпрограмма предназначена для интегрирования линейных систем дифференциальных уравнений с матрицей, имеющей все собственные числа с отрицательными вещественными частями, опыт ее эксплуатации показал, что она является весьма эффективной и для тех систем уравнений, которые имеют собственные числа с нулевой вещественной частью. Экспоненциальный метод, реализованный в данной программе, является пока, по - существу, единственным практически пригодным методом численного интегрирования сильно жестких систем дифференциальных уравнений указанного вида. Значение параметра n, задаваемое при обращении к подпрограмме таким образом, что || a h || ≤ 1, есть, по - существу, число шагов, которые необходимо было бы выполнить, чтобы получить численную устойчивость приближенного решения данной задачи, если ее решать каким - нибудь методом типа Рунге - Кутта. При работе подпрограммы значения параметров m, xn, yn, xk, n сохраняются. На месте матрицы A помещается матрица AH. Если после работы подпрограммы нет необходимости иметь начальное значение решения yn, то параметры yn и y при обращении к ней можно совместить. |
Использование подпрограммы иллюстрируется на примере:
y1' = -20 y1 + y2
y2' = -y1 - 20 y2
y3' = -21 y1 - 19 y2 , 0 < x < 10 ,
y1(0) = 10 , y2(0) = 0 , y3(0) = 0.
Приводится фрагмент вызывающей программы и результаты счета
int main(void)
{
/* Initialized data */
static float a[9] /* was [3][3] */ = { -20.f,-1.f,-21.f,1.f,-20.f,-19.f,
0.f,0.f,0.f };
/* Local variables */
extern int de30r_c(int *, float *, float *, float *, int *, float *,
float *, float *, float *, float *, float *, int *);
static int ierr;
static float e[9] /* was [3][3] */;
static int m, n;
static float r__[3], y[3], r1[9] /* was [3][3] */,
r2[9] /* was [3][3] */, xk, xn, yn[3];
m = 3;
yn[0] = 10.f;
yn[1] = 0.f;
yn[2] = 0.f;
xn = 0.f;
xk = 10.f;
n = 200;
de30r_c(&m, &xn, yn, &xk, &n, a, y, e, r__, r1, r2, &ierr);
printf("\n %5i \n", ierr);
printf("\n %15.7e \n", xk);
printf("\n %15.7e %15.7e %15.7e \n", y[0], y[1], y[2]);
return 0;
} /* main */
Результаты:
y(1) = 0 , y(2) = 0 , y(3) = -10
ierr = 0