|
Текст подпрограммы и версий de07r_c.zip , de07d_c.zip |
Тексты тестовых примеров tde07r_c.zip , tde07d_c.zip |
Интегрирование системы линейных неоднородных уравнений с постоянными коэффициентами и специальным свободным членом экспоненциальным методом.
Решается система Коши для линейной неоднородной системы М обыкновенных дифференциальных уравнений первого порядка с постоянными коэффициентами
Y ' (x) = A * Y(x) + φ(x) ,
Здесь Y = ( y1(x), ... , yM(x) ) ,
A - квадратная вещественная матрица порядка M ,
φ(x) = ( α1 e b1, ... , αM e bM ) ,
где bi = βi x , i = 1, ... , M .
Начальные условия заданы в точке XN:
Y(XN) = YN , YN = ( y1 0, ... , yM 0 ) .
Решение вычисляется в заданной точке XK = XN + T по следующей формуле
Y(XK) = Y(XN+T) = e AT Y(XN) +
M
+ ∑ αm e bm g m ( T ) ,
m=1
где bm = βm XK ,
g m(T) - m-ый столбец матрицы
T
Φm (T) = ∫ e λm dt , λm = ( A - βm T )( T - t )
0
Матричная экспонента e AT и интегралы от матричной экспоненты g m (T) вычисляются с помощью рекуррентных соотношений, которым удовлетворяют матрица e AH и векторы g m (H) для значений скалярного аргумента H = H j , где
H j = ( T / K )K , K = 2 j , 1 ≤ j ≤ N .
Значение N выбирается таким, чтобы выполнялось неpавенство
|| A HN || ≤ 1
(в качестве || A HN || можно взять максимальную сумму модулей элементов матрицы A HN по стpокам). Это значение N задается пользователем при обращении к подпрограмме.
Залеткин С.Ф. Численное решение системы линейных обыкновенных дифференциальных уравнений с постоянными коэффициентами на больших промежутках интегрирования. - В сб.: Вопросы конструирования библиотек программ. M.: Изд - во МГУ, 1985.
int de07r_c (integer *m, real *xn, real *yn, real *xk,
integer *n, integer *q, real *a, real *alp, real *bet,
integer *jstart, real *y, real *e, real *g, real *a1,
real *r, real *r1, integer *ierr)
Параметры
| m - | количество уравнений в системе (тип: целый); |
| xn, yn - | начальные значения аргумента и решения; yn представляет собой одномерный массив длины m (тип xn и yn: вещественный); |
| xk - | значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования); xk может быть больше, меньше или pавно xn (тип: вещественный); |
| n - | показатель степени 2n, равный числу частей, на которые делится отрезок интегрирования [ xn, xk ]; при обращении к подпрограмме параметp n задается таким образом, чтобы || A (xk - xn) / 2n || ≤ 1 (тип: целый); |
| q - |
число членов ряда для аппроксимации матричной экспоненты
q
eAH = ∑ ( AH ) i / i ! , H = ( xk - xn ) / 2n ;
i=0
задается при обращении к подпрограмме (тип: целый);
|
| a - | вещественный двумерный массив порядка m на m, содержащий матрицу коэффициентов системы уpавнений; |
| alp - | вещественный одномерный массив длины m, содержащий коэффициенты α i : alp (i) = α i; |
| bet - | вещественный одномерный массив длины m, содержащий показатели степени β i ; bet (i) = β i; |
| jstart - | целый показатель режима использования подпрограммы, принимающий перед обращением к подпрограмме значение 0; в том случае, когда требуется вычислить решение задачи Коши на равномерной сетке узлов xj = xn + j t, j = 1, 2, ..., то первое обращение к подпрограмме для вычисления решения в первом узле сетки x1 = xn + t производится со значением jstart = 0, а для получения решения в каждом следующем узле x j , j ≥ 2, рекомендуется повтоpно обращаться к подпрограмме со значением jstart = 1, при этом необходимо положить xn = xj - 1, yn = y ( xj - 1 ) (см. также замечания по использованию); |
| y - | вещественный одномерный массив длины m, содержащий вычисленное решение задачи в точке xk; |
|
e, g - a1 | вещественные двумерные рабочие массивы размеpа m на m; |
| r, r1 - | вещественные одномерные рабочие массивы длины m; |
| ierr - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если n < 1, и 66, если q < 1. В этом случае решение задачи прекращается; при этом значения параметров m, xn, yn, xk, n, q, a, alp, bet не изменяются. |
Версии
| de07d_c - | то же, что и de07r_c, но решение задачи вычисляется с удвоенной значностью. При этом параметры xn, yn, xk, a, alp, bet, y, e, g, a1, r, r1 должны иметь тип double. |
Вызываемые подпрограммы
| utde18_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de07r_c. |
| utde19_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы de07d_c. |
Замечания по использованию
|
Значения параметров m, xn, yn, xk, n, q, a, alp, bet при работе подпрограмм сохраняются. При многократном использовании подпрограммы для вычисления решения задачи Коши на равномерной сетке узлов со значением jstart = 1 значения рабочих массивов, задаваемых параметрами a1 и g, не должны изменяться в вызывающей программе между последовательными обращениями к подпрограмме. Если матрица A ( xk - xn ) имеет собственные значения с положительными вещественными частями, то работа программы при вычислении матричной экспоненты eA (xk - xn) при больших значениях | xk - xn | может привести к переполнению арифметического устройства. Использование подпрограммы особенно эффективно для интегрирования жестких систем уравнений. |
y1' = - y1 + 2 y2 - 2 e - k , k = 2x , y1(0) = 5.
y2' = - 3 y1 - 2 y2 + 15 e - x , y2(0) = 1.
Решение вычисляется в двух точках x1 = 5,
x2 = 10
(точное решение задачи
y1 = 5 e- x,
y2 = e - k,
k = 2 x).
int main(void)
{
/* Builtin functions */
double exp(double);
/* Local variables */
extern int de07r_c(int *, float *, float *, float *, int *, int *,
float *, float *, float *, int *, float *, float *,
float *, float *, float *, float *, int *);
static int ierr;
static float a[4] /* was [2][2] */,
e[4] /* was [2][2] */,
g[4] /* was [2][2] */;
static int m, n, q;
static float r__[2], y[2],
a1[4] /* was [2][2] */, r1[2], xk, xn, yn[2];
static int jstart;
static float yt1, yt2, bet[2], alp[2];
#define a_ref(a_1,a_2) a[(a_2)*2 + a_1 - 3]
m = 2;
xn = 0.f;
yn[0] = 5.f;
yn[1] = 1.f;
xk = 5.f;
n = 7;
q = 10;
a_ref(1, 1) = -1.f;
a_ref(1, 2) = 2.f;
a_ref(2, 1) = -3.f;
a_ref(2, 2) = -2.f;
alp[0] = -2.f;
alp[1] = 15.f;
bet[0] = -2.f;
bet[1] = -1.f;
jstart = 0;
de07r_c(&m, &xn, yn, &xk, &n, &q, a, alp, bet, &jstart, y, e, g, a1, r__,
r1, &ierr);
printf("\n %5i \n", ierr);
printf("\n %16.7e \n", xk);
printf("\n %16.7e %16.7e \n", y[0], y[1]);
yt1 = (float)exp((float)(-xk)) * 5.f;
yt2 = (float)exp((float)(xk * -2.f));
printf("\n %16.7e \n", xk);
printf("\n %16.7e %16.7e \n", yt1, yt2);
yn[0] = y[0];
yn[1] = y[1];
xn = xk;
xk += xk;
jstart = 1;
de07r_c(&m, &xn, yn, &xk, &n, &q, a, alp, bet, &jstart, y, e, g, a1, r__,
r1, &ierr);
printf("\n %5i \n", ierr);
printf("\n %16.7e \n", xk);
printf("\n %16.7e %16.7e \n", y[0], y[1]);
yt1 = (float)exp((float)(-xk)) * 5.f;
yt2 = (float)exp((float)(xk * -2.f));
printf("\n %16.7e \n", xk);
printf("\n %16.7e %16.7e \n", yt1, yt2);
return 0;
} /* main */
Результаты:
ierr = 0
xk = 5.00000000000 + 00 y(1) = 3.36897349947-02
xk = 5.00000000000 + 00 yt1 = 3.36897349956-02
y(2) = 4.53999302508-05
yt2 = 4.53999297625-05
ierr = 0
xk = 1.00000000000 + 01 y(1) = 2.26999648806-04
xk = 1.00000000000 + 01 yt1 = 2.26999648813-04
y(2) = 2.06115688405-09
yt2 = 2.06115362244-09