|
Текст подпрограммы и версий de38r_p.zip , de38e_p.zip |
Тексты тестовых примеров tde38r_p.zip , tde38e_p.zip |
Построение начальных значений при интегрировании системы обыкновенных дифференциальных уравнений второго порядка с правой частью, зависящей от производной, методом Штермера с контролем точности.
Для системы M обыкновенных дифференциальных уравнений второго порядка
Y '' = F( X, Y, Y ' ) ,
Y = ( y1, ..., yM ) ,
F = ( f1( X, y1, ..., yM, y1', ..., yM' ), ..., fM( X, y1, ..., yM, y1', ..., yM' ) ) ,
с начальными условиями, заданными в точке XN:
Y(XN) = YN, YN = ( y10, ..., yM0 ) ,
Y ' (XN) = DYN, DYN = ( y10', ..., yM0' ) ,
отыскиваются необходимые при численном интегрировании этой системы методом Штермера значения решения Y и производной Y ' в первом (после начальной точки XN) узле интегрирования ХN + Н, первая разностная производная назад и конечные разности правой части системы до четвертого порядка включительно в том же узле XN + H.
Четвертый порядок разностей соответствует пятому порядку точности используемого метода Штермера.
Указанная процедура вычисления необходимых для счета начальных значений называется разгоном.
Используемый в подпрограмме метод разгона является итерационным способом, опирающимся на три формулы.
Первая формула представляет значение решения в узле ХN + Н через начальные условия (т.е. через его значение и значение его производной в начальной точке XN), а также через конечные расности правой части системы в точке XN.
Вторая формула является предсказывающей формулой Адамса, по которой вычисляется производная решения.
Третья формула является экстраполяционной формулой, применяемой для предсказания приближенного значения решения в используемом для интегрирования методе Штермера.
При этом приближенные значения вычисляются через значения его разностной производной назад.
B качестве критерия точности при разгоне принята оценка первого отброшенного члена экстраполяционной формулы Штермера, при этом все компоненты решения проверяются на точность по мере погрешности, который заключается в следующем: если некоторая компонента приближенного решения по абсолютной величине не меньше некоторой наперед заданной константы P, то контроль точности для этой компоненты ведется по относительной погрешности, иначе - по абсолютной.
Березин И.С., Жидков Н.П. Методы вычислений, т.2. Физматгиз, М., 1960.
Бахвалов H.C., Численные методы, т.1, "Hаука", 1975.
procedure DE38R(F :Proc_F2_DE; M :Integer; IORDER :Integer; var XN :Real;
var YN :Array of Real; var DYN :Array of Real;
HMIN :Real; EPS :Real; P :Real; var H :Real;
var X :Real; var YX :Array of Real;
var DYX :Array of Real; var Z :Array of Real;
var DF :Array of Real; var RFN :Array of Real;
var RF :Array of Real; var R :Array of Real;
var IERR :Integer);
Параметры
| F - |
имя подпрограммы вычисления значений правой
части системы. Первый оператор подпрограммы
должен иметь вид:
procedure F (X :Real; var Y :Array of Real; var DY :Array of Real; var D2Y :Array of Real; M :Integer); Здесь: X, Y, DY - значения независимой, зависимой переменных и производной решения, соответственно; вычисленное значение правой части должно быть помещено в D2Y. B случае системы уравнений, т.е. когда M ≠ 1, параметры Y, DY и D2Y представляют одномерные массивы длины M (тип параметров X, Y, DY и D2Y: вещественный); |
| M - | количество уравнений в системе (тип: целый); |
| IORDER - | порядок точности без единицы того метода Штермера, для которого выполняется разгон и который будет использоваться при интегрировании данной системы уравнений (тип: целый); |
|
XN, YN - DYN | начальные значения аргумента, решения и его производной; в случае системы уравнений (т.е. M ≠ 1) YN и DYN представляют одномерные массивы длиной M (тип: вещественный); |
| HMIN - | минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный); |
| EPS - | допустимая меpа погрешности, с которой тpебуется вычислить все компоненты решения (тип: вещественный); |
| P - | граница перехода, используемая при оценке меры погрешности решения (тип: вещественный); |
| H - | вещественная переменная, содержащая значение шага интегрирования. Если для этого значения шага точность при разгоне достигается, то именно он и реализуется на разгоне, иначе этот шаг уменьшается подпрограммой до тех пор, пока не будет достигнута заданная точность; |
| X - | вещественная переменная, значение которой на выходе из подпрограммы представляет первый (после начальной точки XN) узел интегрирования XN + H, в котоpом вычислены необходимые для интегрирования данной системы уравнений начальные значения; |
|
YX - DYX, Z | одномерные вещественные массивы длины M, в которых запоминаются значения решения, его производной и первой разностной производной назад, вычисленные в узле X; при этом погрешность решения имеет (IОRDЕR + 2) - й порядок по H, а погрешность DYX и Z - (IОRDЕR + 1) - й порядок по H; |
| DF - | двумерный вещественный массив размера M * IORDER, в котоpом запоминаются значения правой части системы и ее разностей до (IОRDЕR - 1) - го порядка включительно, вычисленные в узле X и умноженные на коэффициент H / 12; при этом элемент DF (I, 1) этого массива содержит значение правой части I - ого уравнения системы, а DF (I, J + 1) - ее J - ю разность, погрешность которой имеет (IОRDЕR + 1) - й порядок по H; |
|
RFN - RF, R | одномерные вещественные рабочие массивы длины M. |
| IERR - | целая переменная, значение которой в pезультате работы подпрограммы полагается равным 65, если на разгоне не может быть достигнута требуемая точность EPS; в этом случае разгон можно начать сначала обращением к подпрограмме с новыми значениями паpаметpов H и HMIN. |
Версии
| DE38E - | построение начальных значений при интегрировании системы обыкновенных дифференциальных уpавнений второго порядка, с правой частью, зависящей от производной, методом Штермера с контpолем точности, при этом все промежуточные вычисления выполняются с удвоенным числом значащих цифр. B этом случае параметры XN, YN, DYN, HMIN, EPS, P, H, X, YX, DYX, Z, DF, RFN, RF, R и параметры X, Y, DY и D2Y в подпрограмме F должны иметь тип Extended. |
Вызываемые подпрограммы
| UTDE16 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE38R. |
| UTDE17 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE38E. |
Замечания по использованию
|
B общем случае требуемая точность при разгоне не гарантируется. Значение параметра H, задаваемое при обращении к подпрограмме, должно быть таким, чтобы узел XN + (IORDER - 1) * H не выходил за конец интервала интегрирования. Значение параметра IORDER на входе в подпрограмму полагается равным 5. B дальнейшем предполагается расширить допустимое множество значений этого параметра. При работе подпрограммы и ее версии значения параметров M, IORDER, XN, YN, DYN, HMIN, EPS, P сохраняются. |
y1'' = 0.5 ( y2 - y1 + y1' + y2' ) - 0.5 , y1(0) = 1 , y1' (0) = 1.5 ,
y2'' = y1 - 0.5x , y2(0) = - 1 , y2' (0) = - 0.5 .
Unit TDE38R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE38R_p, DE38R_p;
function TDE38R: String;
implementation
function TDE38R: String;
var
_i,J,I,IERR :Integer;
XN,P,H,HMIN,EPS,X :Real;
YN :Array [0..1] of Real;
DYN :Array [0..1] of Real;
YX :Array [0..1] of Real;
DYX :Array [0..1] of Real;
Z :Array [0..1] of ReАL;
DF :Array [0..9] of Real;
RFN :Array [0..1] of Real;
RF :Array [0..1] of Real;
R :Array [0..1] of Real;
begin
Result := ''; { результат функции }
YN[0] := 1.0;
YN[1] := -1.0;
DYN[0] := 1.5;
DYN[1] := -0.5;
XN := 0.0;
P := 100.0;
H := 0.01;
HMIN := 1.E-11;
EPS := 0.0001;
DE38R(FDE38R,2,5,XN,YN,DYN,HMIN,EPS,P,H,X,YX,DYX,Z,
DF,RFN,RF,R,IERR);
Result := Result + Format('%s',[' IERR=']);
Result := Result + Format('%6d ',[IERR]) + #$0D#$0A;
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 + #$0D#$0A;
for _i:=0 to 1 do
begin
Result := Result + Format('%20.16f ',[DYX[_i]]);
if ( ((_i+1) mod 1)=0 )
then Result := Result + #$0D#$0A;
end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 1 do
begin
Result := Result + Format('%20.16f ',[Z[_i]]);
if ( ((_i+1) mod 4)=0 )
then Result := Result + #$0D#$0A;
end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for J:=1 to 5 do
begin
for I:=1 to 2 do
begin
Result := Result + Format('%20.16f ',[DF[(I-1)+(J-1)*2]]) + #$0D#$0A;
end;
end;
Result := Result + #$0D#$0A;
UtRes('TDE38R',Result); { вывод результатов в файл TDE38R.res }
exit;
end;
end.
Unit FDE38R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;
procedure FDE38R(X :Real; var Y :Array of Real;
var DY :Array of Real; var Z :Array of Real;
M :Integer);
implementation
procedure FDE38R(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.
Результаты:
IERR = 0
X = 1.000000000001 - 02
YX(1) = 1.014949833751 YX(2) = - 1.004949833752
DYX(1) = 1.489950167079 DYX(2) = - 4.899501670834 - 01
Z(1) = 1.494983375078 Z(2) = - 4.949833750843 - 01
DF(1, 1) = - 8.416248614562 - 04 DF(1, 2) = - 8.291528123650 - 06
DF(1, 3) = 8.333255330228 - 08 DF(1, 4) = 8.377316618179 - 10
DF(1, 5) = - 8.498979298110 - 12
DF(2, 1) = 8.416248614518 - 04 DF(2, 2) = 8.291528119209 - 06
DF(2, 3) = - 8.333252932146 - 08 DF(2, 4) = - 8.377654125979 - 10
DF(2, 5) = 8.512301974406 - 12
Точные значения решения и производной такие:
y1 = 1.014949833745 , y2 = - 1.004949833748
y1' = 1.489950167079 , y2' = - 4.899501670798 - 01