|
Текст подпрограммы и версий de15r_p.zip , de15e_p.zip |
Тексты тестовых примеров tde15r_p.zip , tde15e_p.zip |
Выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом Рунге - Кутта - Фельберга.
Выполняется один шаг численного интегрирования системы М обыкновенных дифференциальных уравнений первого порядка
Y ' = F (X, Y) ,
Y = ( y1, ... , yM ) ,
F = ( f1 (X, y1, ... , yM), ... , fM (X, y1, ... , yM) )
методом Рунге - Кутта - Фельберга пятого порядка точности. По заданному значению решения YX в узле Xn вычисляется значение этого решения в узле Xn + H. Bсе компоненты решения вычисляются с контролем точности по меpе погрешности, который заключается в следующем.
Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы Р, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. B качестве абсолютной погрешности решения принимается оценка главного члена асимптотического разложения погрешности метода на шаге, получаемая при вычитании двух выражений, представляющих приближенные значения решения пятого и четвертого порядка точности. При этом на одном шаге интегрирования для определения решения и погрешности используются всего шесть вычислений правой части системы. Значение H может быть меньше или равно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.
Дж.Форсайт, М.Малькольм, К.Моулер, Машинные методы математических вычислений, "Мир", M., 1980.
procedure DE15R(F :Proc_F_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 DY :Array of Real; var R1 :Array of Real;
var R2 :Array of Real; var R3 :Array of Real;
var R4 :Array of Real; var IERR :Integer);
Параметры
| F - |
имя подпрограммы вычисления значений правой
части дифференциального уравнения. Первый
оператоp подпрограммы должен иметь вид: procedure F(X :Real; var Y :Array of Real; var DY :Array of Real; M :Integer) Здесь: X, Y - значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY. B случае системы уравнений, т.е. когда M ≠ 1 , параметры Y и DY представляют массивы длины M (тип параметров X, Y и DY: вещественный); |
| M - | количество уравнений в системе (тип: целый); |
| JSTART - | целый указатель режима использования подпрограммы, имеющий следующие значения: |
| 0,+1 - | выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных и шага интегрирования, заданных параметрами X, YX и H, соответственно; |
| -1 - | повторить последний шаг интегрирования с новыми значениями паpаметpов H и/или HMIN; |
| на выходе из подпрограммы JSTART полагается равным + 1; | |
| HMIN - | минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| EPS - | допустимая меpа погрешности, с которой тpебуется вычислить все компоненты решения (тип: вещественный); |
| P - | граница перехода, используемая при оценке меры погрешности решения (тип: вещественный); |
| YX, X - | заданные вещественные значения решения и соответствующее ему значение аргумента; в pезультате работы подпрограммы в X получается новое значение аргумента, а в YX - соответствующее значение решения; в случае системы уравнений, т.е. когда M ≠ 1, YX задается одномерным массивом длины M; |
| H - | вещественная переменная, содержащая значение шага интегрирования; если для этого значения шага точность приближенного решения достигается, то именно он и реализуется подпрограммой, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность EPS; на выходе из подпрограммы H содержит рекомендуемое подпрограммой значение следущего шага интегрирования, определяемое ею с целью достижения более экономного способа интегрирования; |
| BUL - | логическая переменная, значение которой при обращении к подпрограмме полагается равным TRUE, если заданный в H шаг выводит в конец интервала интегрирования, и FALSE в противном случае; в результате работы подпрограммы BUL pавно FALSE, если вместо исходного шага интегрирования был реализован меньший шаг; в противном случае, т.е. когда был выполнен именно заданный при обращении в H шаг, значение параметра BUL не меняется; |
| XP, YP - | вещественные рабочая переменная и одномерный рабочий массив длины M; значения параметpов XP, YP на выходе из подпрограммы pавны тем значениям, которые имели параметры X, YX, соответственно, при входе в нее (т.е. предыдущий узел и решение в нем); |
|
DY - R1, R2 R3, R4 | одномерные вещественные рабочие массивы длины M; |
| IERR - | целая переменная, значение которой в pезультате работы подпрограммы полагается равным 65, если какая - нибудь компонента pешения не может быть вычислена с требуемой точностью EPS; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров H, HMIN и значением JSTART = - 1 . |
Версии
| DE15E - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений первого порядка методом Рунге - Кутта - Фельберга с удвоенным числом значащих цифр. При этом параметры HMIN, EPS, P, YX, X, H, XP, YP, DY, R1, R2, R3, R4 и параметры X, Y и DY в подпрограмме F должны иметь тип Extended. |
Вызываемые подпрограммы
| UTDE16 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE15R. |
| UTDE17 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE15E. |
Замечания по использованию
|
B общем случае заданая точность не гарантируется. При работе подпрограммы значения параметров M, HMIN, EPS, P сохраняются. При работе подпрограммы счета правой части F значения параметров X, Y и M не должны изменяться. При обращении к подпрограмме со значением JSTART = - 1 в качестве исходных значений аргумента и решения принимаются значения параметров XP и YP, соответственно, т.е. те значения, которые эти параметры получили после самого последнего обращения к подпрограмме с неотрицательным значением JSTART. |
y1' = -y2 - 0.1 x - 0.9
y2' = y1 - 0.1 x - 1.1
y1 = 1 , y2 (0) = -2 , 0 ≤ x ≤ π
Точное решение системы:
y1 = 0.1 x + sin x + , y2 = -0.1 x - cos x -
Приводятся подпрограмма вычисления значений правой части системы и фрагмент вызывающей программы, а также результаты счета.
Unit TDE15R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE15R_p, DE15R_p;
function TDE15R: String;
implementation
function TDE15R: String;
var
M,JSTART,IH,IERR,_i :Integer;
HMIN,EPS,P,X,H,Y1,Y2,ХР :Real;
BUL :Boolean;
YX :Array [0..1] of Real;
YP :Array [0..1] of Real;
DY :Array [0..1] 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;
label
_6,_103,_104,_105,_20;
begin
Result := ''; { результат функции }
M := 2;
JSTART := 0;
HMIN := 1.E-12;
EPS := 1.E-5;
P := 100.0;
YX[0] := 1.0;
YX[1] := -2.0;
X := 0.0;
H := 0.01;
BUL := False;
IH := 0;
_6:
DE15R(FDE15R,M,JSTART,HMIN,EPS,P,YX,X,H,BUL,XP,YP,DY,R1,
R2,R3,R4,IERR);
Result := Result + Format('%s',[' IERR=']);
Result := Result + Format('%3d ',[IERR]) + #$0D#$0A;
if ( IERR <> 0 )
then goto _20;
Y1 := 0.1*X+Sin(X)+1.0;
Y2 := -0.1*X-Cos(X)-1.0;
Result := Result + Format('%20.16f ',[X]) + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 1 do
begin
Result := Result + Format('%20.16f ',[YX[_i]]);
if ( ((_i+1) mod 2)=0 )
then Result := Result + #$0D#$0A;
end;
Result := Result + #$0D#$0A;
Result := Result + Format('%20.16f %20.16f %20.16f ',
[H,Y1,Y2]) + #$0D#$0A;
IH := IH+1;
case IH of
1: goto _6;
2: goto _6;
3: goto _103;
4: goto _104;
5: goto _105;
6: goto _6;
7: goto _20;
end;
_103:
JSTART := -1;
H := -0.005;
goto _6;
_104:
JSTART := -1;
H := 0.02;
goto _6;
_105:
JSTART := -1;
H := 0.01;
goto _6;
_20:
UtRes('TDE15R',Result); { вывод результатов в файл TDE15R.res }
exit;
end;
end.
Unit fde15r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure fde15r(X :Real; var Y :Array of Real; var Z :Array of Real;
M :Integer);
implementation
procedure fde15r(X :Real; var Y :Array of Real; var Z :Array of Real;
M :Integer);
var
R :Real;
begin
R := 0.1*X;
Z[0] := -Y[1]-R-0.9;
Z[1] := Y[0]-R-1.1;
end;
end.
Результаты:
после первого обращения к подпрограмме -
X YX (1) YX (2)
9.999999999991-03 1.010999833334 + 00 -2.000950000420 + 00
H Y1 Y2
1.999999999998-02 1.010999833334 + 00 -2.000950000416 + 00
после второго обращения к подпрограмме -
X YX (1) YX (2)
2.999999999997-02 1.032995500202 + 00 -2.002550033754 + 00
H Y1 Y2
3.999999999996-02 1.032995500202 + 00 -2.002550033747 + 00
после третьего обращения к подпрограмме -
X YX (1) YX (2)
6.999999999994-02 1.076942847336 + 00 -2.004551000264 + 00
H Y1 Y2
7.999999999993-02 1.076942847336 + 00 -2.004551000253 + 00
после четвертого обращения к подпрограмме -
X YX (1) YX (2)
2.499999999998-02 1.027497395913 + 00 -2.002187516282 + 00
H Y1 Y2
-9.999999999991-03 1.027497395913 + 00 -2.002187516275 + 00
после пятого обращения к подпрограмме -
X YX (1) YX (2)
4.999999999995-02 1.054979169268 + 00 -2.003750260403 + 00
H Y1 Y2
3.999999999996-02 1.054979169270 + 00 -2.003750260392 + 00
после шестого обращения к подпрограмме -
X YX (1) YX (2)
3.999999999996-02 1.043989334185 + 00 -2.003200106668 + 00
H Y1 Y2
1.999999999998-02 1.043989334186 + 00 -2.003200106661 + 00
после седьмого обращения к подпрограмме -
X YX (1) YX (2)
5.999999999995-02 1.065964006477 + 00 -2.004200539945 + 00
H Y1 Y2
3.999999999996-02 1.065964006479 + 00 -2.004200539934 + 00