|
Текст подпрограммы и версий de08r_p.zip , de08e_p.zip |
Тексты тестовых примеров tde08r_p.zip , tde08e_p.zip |
Выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом Инглэнда.
Выполняется один шаг численного интегрирования системы M обыкновенных дифференциальных уравнений первого порядка
Y ' = F(X, Y) ,
Y = ( y1, ... , yM ) , F = ( f1(X, y1,..., yM), ... , fM(X, y1,..., yM) )
методом Инглэнда пятого порядка точности. По заданному значению решения YX в узле Xn вычисляется значение этого решения в узле Xn + H. Все компоненты решения вычисляются с контролем точности по мере погрешности следующим образом.
Если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной. В качестве абсолютной погрешности решения принимается оценка главного члена асимптотического разложения погрешности метода на шаге, получаемая при вычитании двух выражений, представляющих приближенные значения решения пятого и четвертого порядка точности. При этом на одном шаге интегрирования для определения решения и погрешности используются всего шесть вычислений правой части системы.
Значение H может быть меньше или равно значению шага интегрирования, задаваемому пользователем при обращении к подпрограмме.
| 1. | Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. - М.: Наука, 1987. |
| 2. | England R. Error estimates for Runge - Kutta type solutions to systems of ordinary differential equations. The Computer Journal. 1969 - v.12, No.2. |
procedure DE08R(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 R1 :Array of Real; var R2 :Array of Real;
var R3 :Array of Real; var IERR :Integer);
Параметры
| F - |
имя подпрограммы вычисления значений правой части
системы дифференциальных уравнений. Первый оператор
подпрограммы должен иметь вид procedure F (X :Real; var Y :Array of Real; var DY :Array of Real; M :Integer); Здесь X, Y - значения независимой и зависимой переменных, соответственно; вычисленное значение правой части должно быть помещено в DY. В случае ситемы уравнений, т.е. когда M ≠ 1, параметры Y и DY представляют одномерные массивы длины M (тип параметров X, Y и DY: вещественный); |
| 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, соответственно, при входе в нее (т.е. предыдущий узел и решение в нем); |
| R1 - | двумерный вещественный рабочий массив размера М*5; |
| R2, R3 - | одномерные вещественные рабочие массивы длины M; |
| IERR - | целая переменная, значение которой в результате работы подпрограммы полагается равным 65, если какая - нибудь компонента решения не может быть вычислена с требуемой точностью EPS; в этом случае последний шаг интегрирования системы можно повторить обращением к подпрограмме с новыми значениями параметров H, HMIN и значением JSTART = - 1 . |
Версии
| DE08E - | выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уравнений первого порядка методом Инглэнда с удвоенным числом значащих цифр. При этом параметры HMIN, EPS, P, YX, X, H, XP, YP, R1, R2, R3 и параметры X, Y и DY в подпрограмме F должны иметь тип Extended. |
Вызываемые подпрограммы
| UTDE20 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE08R. |
| UTDE21 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE08E. |
Замечания по использованию
|
В общем случае заданная точность не гарантируется. При работе подпрограммы значения параметров M, HMIN, EPS, P сохраняются. При работе подпрограммы счета правой части F значения параметров X, Y и M не должны изменяться. При обращении к подпрограмме со значением JSTART = - 1 в качестве исходных значений аргумента и решения принимаются значения параметров XP и YP, соответственно, т.е. те значения, которые эти параметры получили после самого последнего обращения к подпрограмме с неотрицательным значением JSTART. На выходе из подпрограммы в массиве R3 содержится оценка локальной погрешности метода на выполненном шаге. |
y1' = - y1 - 5y2 , y1(0) = 1
y2' = y1 + y2 , y2(0) = 1
Точное решение системы:
y1 = cos2x - 3 sin2x ,
y2 = cos2x + sin2x
Приводятся подпрограмма вычисления значений правой части системы и фрагмент вызывающей программы, а также результаты счета.
Unit TDE08R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE08R_p, DE08R_p;
function TDE08R: String;
implementation
function TDE08R: String;
var
M,JSTART,IH,IERR :Integer;
P,HMIN,H,EPS,X,Y1,Y2,ХР :Real;
BUL :Boolean;
YX :Array [0..1] of Real;
YP :Array [0..1] of Real;
R1 :Array [0..9] of Real;
R2 :Array [0..1] of Real;
R3 :Array [0..1] of Real;
label
_20,_103,_104,_105,_100;
begin
Result := ''; { результат функции }
M := 2;
P := 100.0;
HMIN := 1.E-12;
H := 1.E-2;
EPS := 1.E-5;
JSTART := 1;
X := 0.0;
YX[0] := 1.0;
YX[1] := 1.0;
IH := 0;
BUL := False;
_20:
DE08R (FDE08R,M,JSTART,HMIN,EPS,P,YX,X,H,BUL,XP,YP,
R1,R2,R3,IERR);
Y1 := Cos(2.0*X) - 3.0*Sin(2.0*X);
Y2 := Cos(2.0*X) + Sin(2.0*X);
IH := IH + 1;
Result := Result + Format(' %20.16f %20.16f %20.16f %20.16f %20.16f ',
[X,YX[0],YX[1],R3[0],R3[1]]) + #$0D#$0A;
Result := Result + Format(' %20.16f %20.16f %20.16f ',[H,Y1,Y2]);
case IH of
1: goto _20;
2: goto _20;
3: goto _103;
4: goto _104;
5: goto _105;
6: goto _20;
7: goto _100;
end;
_103:
JSTART := -1;
H := -0.005;
goto _20;
_104:
JSTART := -1;
H := 0.02;
goto _20;
_105:
EPS := 1.E-8;
H := 0.01;
goto _20;
_100:
UtRes('TDE08R',Result); { вывод результатов в файл TDE08R.res }
exit;
end;
end.
Unit fde08r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure fde08r(X :Real; var Y :Array of Real; var Z :Array of Real;
M :Integer);
implementation
procedure fde08r(X :Real; var Y :Array of Real; var Z :Array of Real;
M :Integer);
begin
Z[0] := -Y[0] - 5.0*Y[1];
Z[1] := Y[0] + Y[1];
end;
end.
Результаты:
после первого обращения к подпрограмме -
X YX(1) YX(2)
1.000000000001-02 9.398040065844-01 1.019798673355 + 00
H Y1 Y2
1.000000000001-02 9.398040065835-01 1.019798673356 + 00
контрольный член: -7.983658179000-11 2.673061771929-11
после второго обращения к подпрограмме -
X YX(1) YX(2)
2.000000000001-02 8.792321040974-01 1.039189440839 + 00
H Y1 Y2
1.000000000001-02 8.792321040992-01 1.039189440846 + 00
контрольный член: -8.037659426918-11 2.623323780426-11
после третьего обращения к подпрограмме -
X YX(1) YX(2)
3.00000000000-02 8.183085204928-01 1.058164546403 + 00
H Y1 Y2
1.00000000001-02 8.183085204937-01 1.058164546412 + 00
контрольный член: -8.088818503892-11 2.565069276093-11
после четвертого обращения к подпрограмме -
X YX(1) YX(2)
1.500000000000-02 9.095635331360-01 1.029545533938 + 00
H Y1 Y2
-1.000000000001-02 9.095635331387-01 1.029545533947 + 00
контрольный член: 2.458477865729-12 -8.100187187665-13
после пятого обращения к подпрограмме -
X YX(1) YX(2)
3.000000000000-02 8.183085204928-01 1.058164546403 + 00
H Y1 Y2
1.000000000001-02 8.183085204937-01 1.058164546412 + 00
контрольный член: -8.088818503892-11 2.565059276093-11
после шестого обращения к подпрограмме -
X YX(1) YX(2)
3.249999999997-02 8.030255271615-01 1.062842482492+00
H Y1 Y2
2.500000000001-03 8.030255271660-01 1.062842482504+00
контрольный член: -1.136868377216-13 1.953992523340-14
после седьмого обращения к подпрограмме -
X YX(1) YX(2)
3.499999999997-02 7.877224582326-01 1.067493847573+00
H Y1 Y2
2.500000000001-03 7.877224582389-01 1.067493847588+00
контрольный член: -9.237055564881-14 1.065814103640-14