|
Текст подпрограммы и версий de04r_p.zip , de04e_p.zip |
Тексты тестовых примеров tde04r_p.zip , tde04e_p.zip |
Выполнение одного шага численного интегрирования устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона.
Выполняется один шаг численного интегрирования линейной устойчивой системы M обыкновенных дифференциальных уравнений первого порядка с переменными коэффициентами
(1) Y ' (X) = A(X) * Y(X) + φ(X) ,
Y = ( y1,..., yM ) ,
A(X) = ( ai j(X) ) , i, j = 1, ..., M ,
φ(X) = ( φ1(X), ... , φM(X) )
Предполагается, что среди характеристических корней матрицы A(X) имеются большие по модулю корни, а функция φ(x) является достаточно малой. По заданному значению решения YX в узле xn вычисляется значение этого решения в узле xn + H . Вычисление производится по методу Лоусона.
Метод Лоусона заключается в следующем. Исходная система уравнений с помощью замены искомой функции Y (x) на [ xn, xn + H ] по формуле
Y(x) = exp [ ( x - xn ) A0 ] Z(x),
где A0 - некоторая постоянная матрица, преобразуется в систему уравнений относительно новой неизвестной функции Z (X):
(2) Z ' (x) = A1(x) Z(x) + φ1(x) = exp [ - ( x - xn ) A0 ] { A(x) - A0 }
exp [ ( x - xn ) A0 ] Z(x) + exp [ ( - ( x - xn ) A0 ] φ(x)
Данное преобразование выполняется самой подпрограммой. В качестве матрицы A0 подпрограмма выбирает матрицу A0 = A (xn + H /2). Если шаг H достаточно мал, то преобразование позволяет уменьшить характеристические корни матрицы A1 (x) по сравнению с характеристическими корнями исходной матрицы A (x). Это приводит к уменьшению константы Липшица системы (2) по сравнению с константой Липшица системы (1). Для решения системы (2) применяются формулы классического метода Рунге - Кутта четвертого порядка точности, причем одновременно с решением (2) производится обратное преобразование от функции Z (x) к функции Y (x).
Значение H может быть меньше или равно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме. Все компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. Абсолютная погрешность приближенного решения оценивается по правилу Рунге.
J.Douglas Louson. Generalized Runge - Kutta processes for stable systems with large Lipshitz constants // SIAM Journal on Numerical Analisys. - 1967. - Vol 4, No 3.
procedure DE04R(FA :Proc_F3_DE; FI :Proc_F3_DE; M :Integer; var JSTART :Integer;
HMIN :Real; EPS :Real; P :Real; var YX :Array of Real;
var X :Real; var H :Real; var BUL :Boolean;
var XP :Real; var YP :Array of Real;
var A :Array of Real; var E1 :Array of Real;
var E2 :Array of Real; var E3 :Array of Real;
var W1 :Array of Real; var R1 :Array of Real;
var R2 :Array of Real; var R3 :Array of Real;
var R4 :Array of Real; var R5 :Array of Real;
var R6 :Array of Real; var YD :Array of Real;
var IERR :Integer);
Параметры
| FA - |
подпрограмма вычисления матрицы системы A(X) в
точке X. Первый оператор подпрограммы должен
иметь вид: procedure A (var A :Array of Real; X :Real; M :Integer); Здесь A - двумерный массив размера M*M, в котором помещается матрица системы, вычисленная при значении аргумента X (тип параметров A, X: вещественный); |
| FI - |
подпрограмма вычисления неоднородности правой части системы
φ (X) в любой точке X. Первый
оператор подпрограммы должен иметь вид: procedure FI (var G :Array of Extended; X :Extended; M :Integer); Здесь G - одномерный массив длины M, в который помещается неоднородность правой части системы, вычисленная при значении аргумента X (тип параметров G, X: вещественный); |
| M - | количество уравнений в системе (тип: целый); |
| JSTART - | целый указатель режима использования подпрограммы, имеющий следующие значения: |
| 0,+1 - | выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных и шага интегрирования, заданных параметрами X, YX, и H, соответственно; |
| -1 - | повторить последний шаг интегрирования с новыми значениями параметров H и/или HMIN; на выходе из подпрограммы JSTART полагается равным + 1; |
| HMIN - | минимальное значение абсолютной величины шага, которое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| EPS - | допустимая мера погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
| P - | граница перехода, используемая при оценке меры погрешности решения (тип: вещественный); |
| YX, X - | заданные вещественные значение решения и соответствующее ему значение аргумента; в результате работы подпрограммы в X получается новое значение аргумента, в YX - соответствующее значение решения; в случае системы уравнений, т.е. когда M ≠ 1, YX задается одномерным массивом длины M; |
| H - | вещественная переменная, содержащая значение шага интегрирования; если для этого значения шага точность приближенного решения достигается, то именно он и реализуется подпрограммой, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность EPS; на выходе из подпрограммы H содержит рекомендуемое подпрограммой значение следующего шага интегрирования, определяемое ею с целью достижения более экономного способа интегрирования; |
| BUL - | логическая переменная, значение которой при обращении к подпрограмме полагается равным TRUE, если заданный в H шаг выводит в конец интервала интегрирования, и FALSE в противном случае; в результате работы подпрограммы BUL равно FALSE, если вместо исходного шага интегрирования был реализован меньший шаг; в противном случае, т.е. когда был выполнен именно заданный при обращении в H шаг, значение параметры BUL не меняется; |
| XP, YP - | вещественная рабочая переменная и одномерный рабочий массив длины M, соответственно; значения параметров XP, YP на выходе из подпрограммы равны тем значениям, которые имели параметры X, YX при входе в нее (т.е. предыдущий узел и решение в нем); |
|
A, E1 - E2, E3, W1 | вещественные двумерные рабочие массивы размера M*M; |
|
R1, R2 - R3, R4, R5, R6, YD | вещественные одномерные рабочие массивы длины M; |
| IERR - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров H, HMIN и значением JSTART = - 1 . |
Версии
| DE04E - | выполнение одного шага численного интегрирования линейной устойчивой системы обыкновенных дифференциальных уравнений первого порядка с большой константой Липшица методом Лоусона с удвоенным числом значащих цифр. При этом параметры HMIN, EPS, P, YX, X, H, XP, YP, A, E1, E2, E3, W1, R1, R2, R3, R4, R5, R6, YD и параметры A, G, X в подпрограммах FA и FI должны иметь тип Extended. |
Вызываемые подпрограммы
| AME2R - | подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы DE04R. |
| AME2E - | подпрограмма вычисления матричной экспоненты; вызывается при работе подпрограммы DE04E. |
| UTDE20 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE04R. |
| UTDE21 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE04D. |
|
Кроме того при работе подпрограмм DE04R и DE04E вызываются подпрограммы DE04RD, DE04RM, DE04RP, DE04RQ и DE04ED, DE04EM, DE04EP, DE04EQ соответственно. |
Замечания по использованию
|
Данная подпрограмма предназначена для интегрирования линейных систем, имеющих малую неоднородность φ (x). В общем случае заданная точность не гарантируется. При работе подпрограммы значения параметров M, HMIN, EPS, P сохраняются. При работе подпрограмм FA и FI значения параметров X и M не должны изменяться. При обращении к подпрограмме со значением JSTART = - 1 в качестве исходных значений аргумента и решения принимаются значения параметров XP, YP, соответственно, т.е. те значения, которые эти параметры получили после самого последнего обращения к подпрограмме с неотрицательным значением JSTART. После работы подпрограммы в массиве R1 содержится значение оценки погрешности на шаге, вычисленной по правилу Рунге. Данная подпрограмма и ее версия предназначены также для интегрирования жестких дифференциальных уравнений (1). |
y1' = - ( 2 + x ) y1 /( 1 + x ) + 20 x y2 ,
y1(0) = 0 ,
y2' = -20 x y1 + ( 2 + x ) y2 /( 1 + x ) ,
y2(0) = 1 .
Приводятся подпрограммы вычисления матрицы системы и неоднородной части, фрагмент вызывающей программы, выполняющей несколько шагов из одной точки, а также результаты счета.
Unit TDE04R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FADE04R_p, FIDE04R_p, DE04R_p;
function TDE04R: String;
implementation
function TDE04R: String;
var
M,JSTART,IH,IERR :Integer;
X,HMIN,EPS,P,H,ХР :Real;
BUL :Boolean;
A :Array [0..3] of Real;
E1 :Array [0..3] of Real;
E2 :Array [0..3] of Real;
E3 :Array [0..3] of Real;
W1 :Array [0..3] of Real;
R1 :Array [0..1] of Real;
R2 :Array [0..1] of Real;
R3 :Array [0..1] of Real;
R4 :Array [0..1] of Real;
R5 :Array [0..1] of Real;
R6 :Array [0..1] of Real;
YD :Array [0..1] of Real;
YP :Array [0..1] of Real;
YX :Array [0..1] of Real;
label
_1000,_101,_102,_103,_104,_105,_106,_107;
begin
Result := ''; { результат функции }
M := 2;
X := 0.0;
YX[0] := 0.0;
YX[1] := 1.0;
HMIN := 1.E-10;
EPS := 1.E-7;
P := 100.0;
JSTART := 1;
H := 0.01;
IH := 0;
_1000:
IH := IH + 1;
DE04R (FADE04R,FIDE04R,M,JSTART,HMIN,EPS,P,YX,X,H,BUL,XP,YP,A,E1,E2,
E3,W1,R1,R2,R3,R4,R5,R6,YD,IERR);
Result := Result + Format(' %20.16f %20.16f %20.16f %20.16f ',
[YX[0],YX[1],X,H]) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f ',[R1[0],R1[1]]);
case IH of
1: goto _101;
2: goto _102;
3: goto _103;
4: goto _104;
5: goto _105;
6: goto _106;
7: goto _107;
end;
_101:
EPS := 1.E-7;
goto _1000;
_102:
EPS := 1.0;
goto _1000;
_103:
JSTART := -1;
H := 0.01;
EPS := 1.E-5;
goto _1000;
_104:
JSTART := -1;
H := -0.0001;
EPS := 1.E-2;
goto _1000;
_105:
H := 1.0;
EPS := 1.E-10;
goto _1000;
_106:
EPS := 1.E-4;
goto _1000;
_107:
UtRes('TDE04R',Result); { вывод результатов в файл TDE04R.res }
exit;
end;
end.
Unit fade04r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure fade04r(var A :Array of Real; X :Real; M :Integer);
implementation
procedure fade04r(var A :Array of Real; X :Real; M :Integer);
begin
A[0] := -(2.0 + X)/(1.0 + X);
A[2] := 20.0*X;
A[1] := -A[2];
A[3] := A[0];
end;
end.
Unit fide04e_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure fide04e(var R1 :Array of Extended; X :Extended; M :Integer);
implementation
procedure fide04e(var R1 :Array of Extended; X :Extended; M :Integer);
begin
R1[0] := 0.e0;
R1[1] := 0.e0;
end;
end.
Результаты:
После первого обращения к подпрограмме
X YX(1) YX(2)
1.000000000001-02 9.802471967628-04 9.802468700200-01
H R1(1) R1(2)
2.000000000001-02 0.000000000000+00 6.063298011814-14
после второго обращения к подпрограмме
X YX(1) YX(2)
3.000000000000-02 8.479506692396-03 9.421419716109-01
H R1(1) R1(2)
4.000000000002-02 1.231607408649-14 1.333925562599-12
после третьего обращения к подпрограмме
X YX(1) YX(2)
6.999999999994-02 4.268132414620-02 8.703501915807-01
H R1(1) R1(2)
8.000000000004-02 1.970571853839-12 3.565219230945-11
после четвертого обращения к подпрограмме
X YX(1) YX(2)
3.999999999996-02 1.478074532271-02 9.237177506848-01
H R1(1) R1(2)
2.000000000001-02 0.000000000000+00 6.063298011814-14
после пятого обращения к подпрограмме
X YX(1) YX(2)
2.990000000000-02 8.424732657545-03 9.423281849631-01
H R1(1) R1(2)
-2.000000000004-04 0.000000000000+00 6.063298011814-14
после шестого обращения к подпрограмме
X YX(1) YX(2)
6.115000000000-02 3.314040588361-02 8.858545422663-01
H R1(1) R1(2)
3.125000000000-02 4.395891058563-13 1.097456940137-11
после седьмого обращения к подпрограмме
X YX(1) YX(2)
9.240000000000-02 7.117143016819-02 8.315812963056-01
H R1(1) R1(2)
6.250000000000-02 7.80649190203-13 8.852415097246-12