Текст подпрограммы и версий
de33r_p.zip , de33e_p.zip
Тексты тестовых примеров
tde33r_p.zip , tde33e_p.zip

Подпрограмма:  DE33R (модуль DE33R_p)

Назначение

Построение начальных значений при интегрировании системы обыкновенных дифференциальных уравнений первого порядка методом Адамса без контроля точности.

Математическое описание

Для системы M обыкновенных дифференциальных уравнений первого порядка

          Y ' = F (X, Y) ,

          Y = ( y1, ..., yM ) ,
          F = ( f1 (X, y1,..., yM),..., fM (X, y1,..., yM) )
 с начальными условиями, заданными в точке XN:
          Y(XN) = YN ,     YN = ( y10,..., yM0 ) , 

с помощью модифицированного приема А.Н.Крылова отыскиваются необходимые при численном интегрировании этой системы многошаговым методом Адамса значения решения Y в нескольких первых узлах и конечные разности правой части F (X,Y) системы до некоторого порядка в точке XN. Количество таких узлов и максимальный порядок получаемых разностей определяются в подпрограмме автоматически и соответствуют порядку точности используемого метода Адамса, который задается пользователем при обращении к подпрограмме.

Указанная процедура вычисления необходимых для счета начальных значений, или "фронта Адамса", называется разгоном. Используемый в подпрограмме метод разгона является итерационным способом, опирающимся на экстраполяционную формулу Адамса, записанную в разностной форме.

Беленький В.З. Стандартная программа для интегрирования системы дифференциальных уравнений методом Адамса. Сб. "Вычислительные методы и программирование", вып. 3, Изд-во МГУ, 1965.

Использование

procedure DE33R(F :Proc_F70; var M :Integer; IORDER :Integer; XN :Real;
                var YN :Array of Real; H :Real; var DF :Array of Real;
                var X :Array of Real; var Y :Array of Real;
                var RS :Array of Real; var RF :Array of Real;
                var R :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 - количество уравнений в системе (тип: целый);
IORDER - порядок точности того метода Адамса, для которого выполняется разгон и который будет использоваться при интегрировании данной системы уравнений; IORDER должен быть не больше 10 (тип: целый);
XN, YN - начальные значения аргумента и решения; в случае системы уравнений (т.е. M ≠ 1) YN представляет одномерный массив длины M (тип: вещественный);
H - вещественное значение шага интегрирования;
DF - двумерный вещественный массив размера M*IORDER, в котоpом запоминаются значения правой части системы и ее разностей до порядка (IORDER - 1) включительно, вычисленные в точке XN, при этом элемент DF (I, 1) этого массива содержит значение правой части I - ого уравнения, а DF (I,J + 1) - ее J - ю разность, погрешность вычисления которой имеет (IОRDЕR + 1) - й порядок по H;
X - одномерный вещественный массив длины IORDER, содержащий на выходе из подпрограммы IORDER узлов, включая XN, в которых вычисляются при разгоне приближенные значения решения, при этом элемент X (J) этого массива содержит узел, равный XN + (J - 1) * H;
Y - двумерный вещественный массив размера M*IORDER, в котоpом запоминаются приближенные значения решения, вычисленные для значений аргумента, хранящихся в массиве X, а именно, значению аргумента в X (J) соответствует приближенное решение Y (I,J), при этом погрешность этого решения имеет порядок IORDER по H;
            RS -
         RF, R  
одномерные вещественные рабочие массивы длины M;
IERR - целая переменная, значение которой в результате работы подпрограммы полагается равным 1, если неправильно задан параметр IORDER, т.е. IORDER > 10; в этом случае разгон выполняется для IORDER = 10.

Версии

DE33E - построение начальных значений при интегрировании системы обыкновенных дифференциальных уpавнений первого порядка методом Адамса без контроля точности, при этом все вычисления с действительными числами выполняются с удвоенным числом значащих цифр; в этом случае параметры XN, YN, H, DF, X, Y, RS, RF, R и параметры X, Y и DY в подпрограмме F должны иметь тип Extended.

Вызываемые подпрограммы

UTDE12 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE33R.
UTDE13 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE33E.
  Kpоме того, подпограммы DE33R и DE33E используют рабочие подпрограммы DE28RS и DE28ES, соответственно.

Замечания по использованию

 

При работе подпрограммы и ее версии значения параметров M, IORDER, XN, YN, H сохраняются.

Подпрограмма и ее версия могут использоваться как для разгона, так и непосредственно для численного решения системы уравнений на всем интервале интегрирования.

Пример использования

           y1'  =   y2
           y2'  =  -y1

          y1 (3/4 π)  =   √2 /2
          y2 (3/4 π)  =  -√2 /2

Unit TDE33R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE33R_p, DE33R_p;

function TDE33R: String;

implementation

function TDE33R: String;
var
IORDER,M,J,I,IERR :Integer;
H,XN,R1,R2 :Real;
YN :АRray [0..1] of Real;
DF :Array [0..9] of Real;
Y :Array [0..9] of Real;
X :Array [0..4] of Real;
RS :Array [0..1] of Real;
RF :Array [0..1] of Real;
R :Array [0..1] of Real;
begin
Result := '';  { результат функции }
IORDER := 5;
M := 2;
H := 0.01;
XN := 0.75*3.14159265359;
YN[0] := Sqrt(2.0)/2.0;
YN[1] := -YN[0];
DE33R(FDE33R,M,IORDER,XN,YN,H,DF,X,Y,RS,RF,R,IERR);
XN := X[4];
R1 := Sin(XN);
R2 := Cos(XN);
Result := Result + Format(' %20.16f %20.16f %20.16f  %20.16f  %20.16f ',
 [X[4],Y[8],Y[9],R1,R2]) + #$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('TDE33R',Result);  { вывод результатов в файл TDE33R.res }
exit;
end;

end.

Unit fde33r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;

procedure fde33r(X :Real; var Y :Array of Real; var DY :Array of Real;
                M :Integer);

implementation

procedure fde33r(X :Real; var Y :Array of Real; var DY :Array of Real;
                M :Integer);
begin
DY[0] := Y[1];
DY[1] := -Y[0];
end;

end.

Результаты:

          IERR  =  0
    
          X(5)  =  2.396194490 + 00
          Y(1, 5)  =   6.782644418-01
          Y(2, 5)  =  -7.348179005-01
   
          DF(1, 1)  =  -7.07106781186-03
          DF(1, 2)  =  -7.10630480683-05
          DF(1, 3)  =   6.99988802921-07
          DF(1, 4)  =   7.18290493751-09
          DF(1, 5)  =  -7.27808924239-11

          DF(2, 1)  =  -7.07106781186-03
          DF(2, 2)  =   7.03559472512-05
          DF(2, 3)  =   7.14142103675-07
          DF(2, 4)  =  -6.97064450605-09
          DF(2, 5)  =  -6.86739554112-11