|
Текст подпрограммы и версий de30r_p.zip , de30e_p.zip |
Тексты тестовых примеров tde30r_p.zip , tde30e_p.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.
procedure DE30R(M :Integer; XN :Real; var YN :Array of Real;
var XK :Real; N :Integer; var A :Array of Real;
var Y :Array of Real; var E :Array of Real;
var R :Array of Real; var R1 :Array of Real;
var R2 :Array of Real; var IERR :Integer);
Параметры
| M - | количество уравнений в системе (тип: целый); |
| XN, YN - | начальные значения аргумента и решения; YN представляет одномерный массив длины M (тип XN и YN: вещественный); |
| XK - | значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования). XK может быть больше, меньше или pавно XN (тип: вещественный); |
| N - | целое число равных частей, на которые делится отрезок интегрирования; при обращении к подпрограмме параметр N выбирается таким образом, чтобы || A (XK - XN) / 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, А не изменяются. |
Версии
| DE30E - | вычисление решения задачи Коши для линейной однородной системы обыкновенных дифференциальных уравнений первого порядка с постоянными коэффициентами экспоненциальным методом с расширенной (Extended) точностью. При этом параметры XN, YN, XK, A, Y, E, R1, R2, R должны иметь тип Extended. |
Вызываемые подпрограммы
| UTDE10 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE30R. |
| UTDE11 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE30E. |
Замечания по использованию
|
Хотя подпрограмма предназначена для интегрирования линейных систем дифференциальных уравнений с матрицей, имеющей все собственные числа с отрицательными вещественными частями, опыт ее эксплуатации показал, что она является весьма эффективной и для тех систем уравнений, которые имеют собственные числа с нулевой вещественной частью. Экспоненциальный метод, реализованный в данной программе, является пока, по - существу, единственным практически пригодным методом численного интегрирования сильно жестких систем дифференциальных уравнений указанного вида. Значение параметра N, задаваемое при обращении к подпрограмме таким образом, что || A H || ≤ 1, есть, по - существу, число шагов, которые необходимо было бы выполнить, чтобы получить численную устойчивость приближенного решения данной задачи, если ее решать каким - нибудь методом типа Рунге - Кутта. При работе подпрограммы значения параметров M, XN, YN, XK, N сохраняются. На месте матрицы А помещается матрица 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.
Приводится фрагмент вызывающей программы и результаты счета
Unit TDE30R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, DE30R_p;
function TDE30R: String;
implementation
function TDE30R: String;
var
M,N,_i,IERR :Integer;
XN,ХК :Real;
YN :Array [0..2] of Real;
Y :Array [0..2] of Real;
E :Array [0..8] of Real;
R :Array [0..2] of Real;
R1 :Array [0..8] of Real;
R2 :Array [0..8] of Real;
const
A :Array [0..8] of Real = ( -20.0,-1.0,-21.0,1.0,-20.0,-19.0,0.0,0.0,0.0 );
begin
Result := ''; { результат функции }
M := 3;
YN[0] := 10.0;
YN[1] := 0.0;
YN[2] := 0.0;
XN := 0.0;
ХК := 10.0;
N := 200;
DE30R(M,XN,YN,XK,N,A,Y,E,R,R1,R2,IERR);
Result := Result + Format('%s',[' XK=']);
Result := Result + Format('%20.16f ',[XK]);
Result := Result + Format('%s',[#$0D#$0A + ' Y=']);
Result := Result + #$0D#$0A;
for _i:=0 to 2 do
begin
Result := Result + Format('%20.16f ',[Y[_i]]);
if ( ((_i+1) mod 3)=0 )
then Result := Result + #$0D#$0A;
end;
Result := Result + #$0D#$0A;
UtRes('TDE30R',Result); { вывод результатов в файл TDE30R.res }
exit;
end;
end.
Результаты:
Y(1) = 0 , Y(2) = 0 , Y(3) = -10
IERR = 0