БЧА НИВЦ МГУ. DE35R_P. Задачи Коши для системы уравнений второго порядка

Текст подпрограммы и версий
de35r_p.zip , de35e_p.zip

Тексты тестовых примеров
tde35r_p.zip , tde35e_p.zip

Подпрограмма: DE35R (модуль DE35R_p)

Назначение

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

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

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

          Y '' = F (X, Y)   ,
          Y = ( y₁, ..., y_M )   ,
          F = ( f₁ (X, y₁, ..., y_M), ..., f_M (X, y₁, ..., y_M) )   ,

с начальными условиями, заданными в точке XN:

          Y (XN) = YN ,         YN = ( y₁₀, ..., y_M0 )   ,
          Y ' (XN) = DYN ,     DYN = ( y₁₀', ..., y_M0' )   ,

отыскиваются необходимые при численном интегрировании этой системы многошаговым методом Штермера значение решения Y в первом (после начальной точки XN) узле интегрирования ХN + Н, его конечная разность первого порядка и конечные разности правой части системы до четвертого порядка включительно в том же узле XN + H. Четвертый порядок разностей соответствует шестому порядку точности используемого метода Штермера.

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

Используемый в подпрограмме метод разгона является итерационным способом, опирающимся на пару формул, одна из которых представляет значение решения в узле XN + H через начальные условия (т.е. через его значение и значение его производной в начальной точке XN), а также через конечные разности правой части системы в точке XN, а другая является экстраполяционной формулой, применяемой для предсказания приближенного значения решения в используемом для интегрирования методе Штермера.

Березин И.С., Жидков Н.П. Методы вычислений, т.2. Физматгиз, М., 1960.

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

procedure DE35R(F :Proc_F70; var M :Integer; IORDER :Integer;
                var XN :Real; var YN :Array of Real;
                var DYN :Array of Real; H :Real; var X :Real;
                var YX :Array of Real; var DY :Array of Real;
                var DF :Array of Real; var RFN :Array of Real;
                var RF :Array of Real; var R :Array of Real);

   Параметры

F -	имя подпрограммы вычисления значений правой части системы. Первый оператор подпрограммы должен иметь вид:

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 -	порядок точности без единицы того метода Штермера, для которого выполняется разгон и который будет использоваться при интегрировании данной системы уравнений (тип: целый);
XN, YN - DYN	начальные значения аргумента, решения и его производной; в случае системы уравнений (т.е. M ≠ 1) YN и DYN представляют одномерные массивы длиной M (тип: вещественный);
H -	вещественное значение шага интегрирования;
X -	вещественная переменная, значение которой на выходе из подпрограммы представляет первый (после начальной точки XN) узел интегрирования XN + H, в котоpом вычислены необходимые для интегрирования данной системы уравнений начальные значения;
YX, DY -	одномерные вещественные массивы длиной M, в которых запоминаются значения решения и его первой разности, вычисленные в узле X, при этом погрешность решения и разности имеет (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.

   Версии

DE35E -

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

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

Подпрограммы DE35R и DE35E используют рабочие подпрограммы DE28RS и DE28ES, соответственно.

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

Значение параметра H, задаваемое при обращении к подпрограмме, должно быть таким, чтобы узел XN + (IORDER - 1) * H не выходил за конец интервала интегрирования.

Значение параметра IORDER на входе в подпрограмму полагается равным 5. В дальнейшем предполагается расширить допустимое множество значений этого параметра.

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

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

      y₁'' = - y₁   ,
      y₂'' = - y₂   ,
      y₁ (3/4 * π) = √2 / 2   ,       y₂ (3/4 * π) = - √2 / 2   , 
      y₁' (3/4 * π) = - √2 / 2   ,    y₂' (3/4 * π) = - √2 / 2   .

Unit TDE35R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE35R_p, DE35R_p;

function TDE35R: String;

implementation

function TDE35R: String;
var
IORDER,M,_i,J,I :Integer;
H,XN,X :Real;
YN :Array [0..1] of Real;
DYN :Array [0..1] of Real;
Y :Array [0..1] of Real;
DY :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;
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];
DYN[0] := YN[1];
DYN[1] := -YN[0];
DЕ35R(FDE35R,M,IORDER,XN,YN,DYN,H,X,Y ,DY,DF,RFN,RF,DELTY);
Result := Result + Format(' %20.16f ',[X]) + #$0D#$0A;
RЕSULt := Result + #$0D#$0A;
for _i:=0 to 1 do
 begin
  Result := Result + Format('%20.16f ',[Y[_i]]);
  if ( ((_i+1) mod 2)=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 ',[DY[_i]]);
  if ( ((_i+1) mod 2)=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('TDE35R',Result);  { вывод результатов в файл TDE35R.res }
exit;
end;

end.

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

procedure fde35r(T :Real; var Y :Array of Real; var Z :Array of Real;
                M :Integer);

implementation

procedure fde35r(T :Real; var Y :Array of Real; var Z :Array of Real;
                M :Integer);
begin
Z[0] := -Y[0];
Z[1] := -Y[1];
end;

end.

Результаты:

       X = 2.366194490 + 00

       Y(1) = 7.000004762 - 01            Y(2) = - 7.141423761 - 01

       DY(1) = - 7.106305006 - 03          DY(2) = - 7.035594917 - 03

       DF(1, 1) = - 5.83333730151 - 06     DF(1, 2) =    5.92192083634 - 08
       DF(1, 3) =   5.89251453187 - 10      DF(1, 4) = - 5.86490578325 - 12
       DF(1, 5) = - 5.76760861293 - 14

       DF(2, 1) =   5.95118646749 - 06      DF(2, 2) =   5.86299576250 - 08
       DF(2, 3) = - 5.89250072347 - 10      DF(2, 4) = - 5.92374760355 - 12
       DF(2, 5) =   6.01324545713 - 14