|
Текст подпрограммы и версий (Паскаль) de46r_p.zip , de46e_p.zip |
Тексты тестовых примеров (Паскаль) tde46r_p.zip , tde46e_p.zip |
Выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений второго порядка с правой частью, зависящей от производной, методом Штермера с контролем точности.
Выполняется один шаг численного интегрирования системы М обыкновенных дифференциальных уравнений второго порядка
Y '' = F ( X, Y, Y ' ) ,
Y = ( y1, ..., yM ) ,
F = ( f1 ( X, y1, ..., yM, y1', ..., yM' ), ..., fM (X, y1, ..., yM, y '1, ..., y 'M ) )
многозначным предсказывающе - исправляющим методом.
При этом приближенное значение решения Y вычисляется методом Штермера пятого порядка точности, производной - методом Адамса пятого порядка точности.
Суть метода состоит в том, что сначала по явным формулам строятся предсказанные значения решения и производной, которые затем исправляются по неявным формулам.
По значениям решения YX и производной DYX в узле Xn вычисляются значения решения и производной в узле Xn + H.
Bместе с тем подпрограмма определяет также всю дополнительную информацию, необходимую для вычисления приближенных решения и производной в любом следующем узле, в частности, первую разностную производную решения назад, значение которой используется для вычисления приближенного решения YX.
Bсе компоненты решения вычисляются с контролем точности по мере погрешности, который заключается в следующем. Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.
Значение H может быть меньше или pавно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.
Бахвалов H.C., Численные методы, т.1, "Hаука", 1975.
procedure DE46R(F :Proc_F2_DE; M :Integer; var JSTART :Integer;
HMIN :Real; HMAX :Real; EPS :Real; P :Real;
var YX :Array of Real; var DYX :Array of Real;
var X :Real; var H :Real; var BUL :Boolean;
var Z :Array of Real; var DELTY :Array of Real;
var DF :Array of Real; var RFN :Array of Real;
var RF :Array of Real; var YP :Array of Real;
var DYP :Array of Real; var ZP :Array of Real;
var IERR :Integer);
Параметры
| F - | имя подпрограммы вычисления значений правой части системы. Первый оператор подпрограммы должен иметь вид: |
procedure F (X :Real; var Y :Array of Real; var DY :Array of Real;
var DY2 :Array of Real; M :Integer);
Здесь: X, Y, DY - значения независимой, зависимой переменных и производной решения, соответственно; вычисленное значение правой части должно быть помещено в D2Y. B случае системы уравнений, т.е. когда M ≠ 1 - параметры Y, DY, D2Y представляют одномерные массивы длины M (тип параметров X, Y, DY и D2Y: вещественный); |
| M - | количество уравнений в системе (тип: целый); |
| JSTART - | целый указатель режима использования подпрограммы, имеющий следующие значения: |
| 0 - | первое обращение к подпрограмме должно быть выполнено с нулевым значением JSTART; |
| +1 - | выполнить один шаг интегрирования системы дифференциальных уравнений для значений независимой и зависимой переменных, производной и шага интегрирования, заданных параметрами X, YX, DYX и H, соответственно; |
| -1 - | повторить последний шаг интегрирования с новыми значениями праметров H и / или HMIN; |
| на выходе из подпрограммы JSTART полагается равным + 1; |
| HMIN - | минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| HMAX - | максимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| EPS - | допустимая меpа погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
| P - | граница перехода, используемая при оценке меры погрешности приближенного решения (тип: вещественный); |
|
YX - DYX, X | задаваемые вещественные значения решения, производной и соответствующее им значение аргумента; в результате работы подпрограммы в X получается новое значение аргумента, а в YX и DYX соответствующие ему значения решения и производной; в случае системы уравнений, т.е. когда M ≠ 1, YX и DYX задаются одномерными массивами длины M; |
| H - | вещественная переменная, содержащая значение шага интегрирования; если для этого значения шага точность приближенного pешения достигается, то именно он и реализуется подрограммой, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность EPS; на выходе из подпрограммы H содержит pекомендуемое подпрограммой значение следующего шага интегрирования, определяемое ею с целью достижения более экономного способа интегрирования; |
| BUL - | логическая переменная, значение которой при обращении к подпрограмме полагается равным TRUE, если заданный в H шаг выводит в конец интервала интегрирования, и FALSE в противном случае; в результате работы подпрограммы BUL pавно FALSE, если вместо исходного шага интегрирования был реализован меньший шаг; в противном случае, т.е. когда был выполнен именно заданный при обращении в H шаг, значение паpаметpа BUL не меняется; |
|
Z - DELTY RFN, RF ZP | одномерные вещественные рабочие массивы длины M; |
| DF - | двумерный вещественный рабочий массив размеpа M * 5, в котоpом запоминаются значения правой части системы и ее разностей до четвертого порядка включительно, умноженные на коэффициент H / 12; |
| YP, DYP - | одномерные вещественные рабочие массивы длины M, представляющие значения решения и производной в предыдущем узле интегрирования; |
| IERR - | целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом: |
| IERR=65 - | когда интегрирование системы выполнено с заданным в HMIN минимальным шагом, но требуемая точность полученного при этом значения YX решения не достигается; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров H и HMIN и со значением JSTART = - 1; |
| IERR=66 - | когда решение не может быть вычислено с требуемой точностью EPS при первом обращении к подпрограмме (т.е. со значением JSTART = 0); в этом случае интегрирование системы можно начать сначала обращением к подпрограмме с новыми значениями параметpов H и HMIN и со значением JSTART = 0 . |
Версии
| DE46E - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений второго порядка с правой частью, зависящей от производной, методом Штермера с удвоенным числом значащих цифр. При этом параметры HMIN, HMAX, EPS, P, YX, DYX, X, H, Z, DELTY, DF, RFN, RF, YP, DYP, ZP и параметры X, Y, DY и D2Y в подпрограмме F должны иметь тип Extended. |
Вызываемые подпрограммы
|
DE38R - DE38E | построение начальных значений при интегрировании системы обыкновенных дифференциальных уравнений второго порядка с правой частью, зависящей от производной, методом Штермера с контролем точности. |
| UTDE16 - UTDE17 | подпрограммы выдачи диагностических сообщений. |
|
Подпрограммы DE38R и UTDE16 вызываются при работе подпрограммы DE46R, а подпрограммы DE38E и UTDE17 - при работе подпрограммы DE46E. Kpоме того, DE46R и DE46E используют рабочие подпрограммы DE28RS, DE48RP и DE28ES, DE48EP. |
Замечания по использованию
|
B общем случае заданая точность решения не гарантируется. Приближенное значение производной на точность не проверяется. При работе подпрограммы и ее версии значения параметров M, HMIN, HMAX, EPS, P сохраняются. При многократном использовании подпрограммы и ее версии для вычисления решения системы уравнений на отрезке значения параметров M, YX, DYX, X и значения рабочих массивов, задаваемых параметрами Z, DF, YP, DYP, ZP, не должны изменяться в вызывающей программе между последовательными обращениями к подпрограмме. Значение параметра HMAX не должно быть слишком большим, т.к. в противном случае величина шага интегрирования может достичь таких размеров, которые приведут к абсолютной неустойчивости численного метода. Tак как подпрограммы DE46R и DE46E используют глобальные записи (структуры данных) с именами _COM48R и _COM48D, соответственно, для хранения промежуточных значений, поэтому пользователь не должен портить элементы этих записей. |
y1'' = 0.5 ( y2 - y1 + y1' + y2' ) - 0.5 ,
y1 (0) = 1 , y1' (0) = 1.5 ;
y2'' = y1 - 0.5 x ,
y2 (0) = - 1 , y2' = - 0.5 .
Unit TDE46R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE46R_p, DE46R_p;
function TDE46R: String;
implementation
function TDE46R: String;
var
M,JSTART,IH,IERR,_i :Integer;
X,HMIN,HMAX,EPS,P,H,Y1,Y2,DY1,DY2 :Real;
BUL :Boolean;
YX :Array [0..1] of Real;
DYX :Array [0..1] of Real;
Z :Array [0..1] of Real;
DELTY :Array [0..1] of Real;
DF :Array [0..9] of Real;
RFN :Array [0..1] of Real;
RF :Array [0..1] of Real;
YP :Array [0..1] of Real;
DYP :Array [0..1] of Real;
ZP :Array [0..1] of Real;
label
_6,_102,_103,_104,_20;
begin
Result := ''; { результат функции }
M := 2;
YX[0] := 1.0;
YX[1] := -1.0;
DYX[0] := 1.5;
DYX[1] := -0.5;
X := 0.0;
JSTART := 0;
HMIN := 1.E-11;
НМАХ := 0.1;
EPS := 0.0001;
P := 100.0;
H := 0.01;
BUL :=False;
IH := 0;
_6:
DE46R(FDE46R,M,JSTART,HMIN,HMAX,EPS,P,YX,DYX,X,H,BUL,Z,DELTY,DF,
RFN,RF,YP,DYP,ZP,IERR);
Result := Result + Format('%s',[' IERR=']);
Result := Result + Format('%3d ',[IERR]) + #$0D#$0A;
if ( IERR <> 0 )
then goto _20;
IH := IH+1;
Y1 := Sin(X)+Cos(X)+0.5*X;
Y2 := -Sin(X)-Cos(X)+0.5*X;
DY1 := Cos(X)-Sin(X)+0.5;
DY2 := -Cos(X)+Sin(X)+0.5;
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 1)=0 )
then Result := Result + #$0D#$0A;
end;
Result := Result + #$0D#$0A;
Result := Result + Format('%20.16f %20.16f %20.16f',
[Y1,Y2,H]) + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 1 do
begin
Result := Result + Format('%20.16f ',[DYX[_i]]);
if ( ((_i+1) mod 3)=0 )
then Result := Result + #$0D#$0A;
end;
Result := Result + #$0D#$0A;
Result := Result + Format('%20.16f %20.16f ',
[DY1,DY2]) + #$0D#$0A;
case IH of
1: goto _6;
2: goto _6;
3: goto _102;
4: goto _103;
5: goto _104;
6: goto _6;
7: goto _20;
end;
_102:
JSTART := -1;
H := -0.005;
goto _6;
_103:
JSTART := -1;
H := 0.02;
goto _6;
_104:
JSTART := -1;
H := 0.01;
goto _6;
_20:
UtRes('TDE46R',Result); { вывод результатов в файл TDE46R.res }
exit;
end;
end.
Unit FDE46R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure FDE46R(X :Real; var Y :Array of Real;
var DY :Array of Real; var Z :Array of Real;
M :Integer);
implementation
procedure FDE46R(X :Real; var Y :Array of Real;
var DY :Array of Real; var Z :Array of Real;
M :Integer);
begin
Z[0] := 0.5*(Y[1]+DY[0]-Y[0]+DY[1])-0.5;
Z[1] := Y[0]-0.5*X;
end;
end.
Результаты:
после первого обращения к подпрограмме -
X YX(1) YX(2)
1.000000000001 - 02 1.014949833751 + 00 - 1.004949833752 + 00
Y1 Y2
1.014949833745 + 00 - 1.004949833748 + 00
H DYX(1) DYX(2)
1.000000000001 - 02 1.489950167079 + 00 - 4.899501670834 - 01
DY1 DY2
1.489950167079 + 00 - 4.899501670798-01
после второго обращения к подпрограмме -
X YX(1) YX(2)
2.000000000001 - 02 1.029798673357 + 00 - 1.009798673364 + 00
Y1 Y2
1.029798673355 + 00 - 1.009798673360 + 00
H DYX(1) DYX(2)
2.000000000001 - 02 1.479801339963 + 00 - 4.798013399754 - 01
DY1 DY2
1.479801339969 + 00 - 4.798013399704-01
после третьего обращения к подпрограмме -
X YX(1) YX(2)
4.000000000002 - 02 1.059189440843 + 00 - 1.019189440856 + 00
Y1 Y2
1.059189440844 + 00 - 1.019189440849 + 00
H DYX(1) DYX(2)
4.000000000002 - 02 1.459210772458 + 00 - 4.592107724779-01
DY1 DY2
1.459210772471 + 00 - 4.592107724729-01
после четвертого обращения к подпрограмме -
X YX(1) YX(2)
1.500000000000 - 02 1.022386939610 + 00 - 1.007386939622 + 00
Y1 Y2
1.022386939612 + 00 - 1.007386939615 + 00
H DYX(1) DYX(2)
- 1.000000000001 - 02 1.484888064588 + 00 - 4.848880646064 - 01
DY1 DY2
1.484888064599 + 00 - 4.848880646005 - 01
после пятого обращения к подпрограмме -
X YX(1) YX(2)
4.000000000002 - 02 1.059189440841 + 00 - 1.019189440856 + 00
Y1 Y2
1.059189440844 + 00 - 1.019189440849 + 00
H DYX(1) DYX(2)
4.000000000002 - 02 1.459210772457 + 00 - 4.592107724784 - 01
DY1 DY2
1.459210772471 + 00 - 4.592107724729 - 01
после шестого обращения к подпрограмме -
X YX(1) YX(2)
3.000000000000 - 02 1.044545533945 + 00 - 1.014545533959 + 00
Y1 Y2
1.044545533947 + 00 - 1.014545533950 + 00
H DYX(1) DYX(2)
2.000000000001 - 02 1.469554533529 + 00 - 4.695545335508 - 01
DY1 DY2
1.469554533543 + 00 - 4.695545335444 - 01
после седьмого обращения к подпрограмме -
X YX(1) YX(2)
5.000000000001 - 02 1.073729429656 + 00 - 1.023729429677 + 00
Y1 Y2
1.073729429661 + 00 - 1.023729429664 + 01
H DYX(1) DYX(2)
4.000000000002 - 02 1.448771091098 + 00 - 4.487710911303 - 01
DY1 DY2
1.448771091120 + 00 - 4.487710911217 - 01