Текст подпрограммы и версий
de41r_p.zip , de41e_p.zip 		Тексты тестовых примеров
tde41r_p.zip , tde41e_p.zip

Подпрограмма:  DE41R (модуль DE41R_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 ' (XN) = DYN   ,    DYN = ( y10', ..., yM0' )   ,   -

многошаговым предсказывающе - исправляющим методом пятого порядка точности. Суть метода состоит в том, что сначала по явной формуле Адамса - Штермера строятся предсказанные значения приближенного решения, которые затем исправляются по неявной формуле Адамса - Штермера.

Решение вычисляется в одной точке  XK, которая является концом интервала интегрирования. Интегрирование ведется с постоянным шагом, значение которого задается пользователем при обращении к подпрограмме.

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

Жоголев E.A. Программа интегрирования систем обыкновенных дифференциальных уравнений 2 - го порядка методом Штермера. Сб. "Вычислительные методы и программирование", вып. 1 , Изд - во МГУ, 1962.

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

procedure DE41R(F :Proc_F70; var M :Integer; XN :Real;
                var YN :Array of Real; var DYN :Array of Real;
                XK :Real; var H :Real; var Y :Array of Real;
                var DY :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 IERR :Integer);


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

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 - количество уравнений в системе (тип: целый);
         XN, YN -
         DYN   		начальные значения аргумента, решения и его производной, соответственно; в случае системы уравнений (т.е.  M ≠ 1)  YN и  DYN представляют одномерные массивы длиной  M (тип: вещественный);
XK - значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования).  XK может быть больше, меньше или pавно  XN (тип: вещественный);
H - вещественная переменная, содержащая значение шага интегрирования; может задаваться с учетом направления интегрирования, т.е. положительным, если  XK > XN, отрицательным, если  XK < XN, или без такого учета в виде абсолютной величины;
Y - искомое решение задачи Коши, вычисленное подпрограммой для значения аргумента  XK. Для системы уравнений (когда  M ≠ 1)  Y задается массивом длиной  M. B случае совпадения значений параметров  XN и  XK значение   Y полагается равным начальному значению  YN (тип: вещественный);
            DY -
       DELTY  
       RFN, RF  
         YP, DYP   		вещественные одномерные рабочие массивы длиной  M;
DF - двумерный вещественный рабочий массив размеpа  M * 5;
IERR - целая переменная, значение которой в pезультате работы подпрограммы полагается равным 65, если при заданном значении шага интегрирования не может быть достигнута сходимость численного решения к точному значению решения; в этом случае интегрирование системы можно повторить обращением к подпрограмме с новым (меньшим) значением параметра  H. 
   Версии
DE41E - вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений второго порядка в конце интервала интегрирования методом Штермера с постоянным шагом, при этом все вычисления над действительными числами ведутся с удвоенным числом значащих цифр. B этом случае параметры  XN, YN, DYN, XK, H, Y, DY, DELTY, DF, RFN, RF, YP, DYP и параметры  X, Y и  DY в подпрограмме  F должны иметь тип Extended. 
   Вызываемые подпрограммы
DE40R - выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений второго порядка методом Штермера без контроля точности; вызывается при работе подпрограммы DE41R.
DE40E - выполнение одного шага численного интегрирования системы обыкновенных дифференциальных уpавнений второго порядка методом Штермера без контроля точности, при этом все вычисления над действительными числами ведутся с удвоенным числом значащих цифр; вызывается при работе подпрограммы DE41E.
UTDE12 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE41R.
UTDE13 - подпрограмма выдачи диагностических сообщений при работе подпрограммы DE41E. 
   Замечания по использованию
 

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

Если после работы подпрограммы нет необходимости иметь начальные значения решения  YN и его производной  DYN, то параметры  YN и  Y, а в случае производной  -  параметры   DYN и  DY, при обращении к ней можно совместить.

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

     y1'' = - 6y1 - 7y2   ,
     y2'' = - 3y1 - 2y2 + 2x   ,
     y1 (0) = 0   ,   y2 (0) = 0   ,   y1' (0) = - 1 / 9   ,   y2' (0) = - 3   .

Приводятся подпрограмма вычисления значений правой части системы и фрагмент вызывающей программы, выполняющей дважды интегрирование системы, сначала на отрезке [0, 5] слева направо, затем на отрезке [-5, 0] справа налево, а также результаты счета. 

Unit TDE41R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FDE41R_p, DE41R_p;

function TDE41R: String;

implementation

function TDE41R: String;
var
M,K,IERR :Integer;
H,XN,XK,Y3,Y1,Y2 :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;
YP :Array [0..1] of Real;
DYP :Array [0..1] of Real;
label
_200;
begin
Result := '';  { результат функции }
H := 0.01;
M := 2;
XN := 0.0;
ХК := 5.0;
YN[0] := 0.0;
YN[1] := 0.0;
DYN[0] := -1.0/9.0;
DYN[1] := -3.0;
Y3 := (Exp(5.0)-Exp(-5.0))/3.0;
Y1 := Y3-7.0*Sin(15.0)/9.0+70.0/9.0;
Y2 := -Y3-Sin(15.0)/3.0-20.0/3.0;
for K:=1 to 2 do
 begin
  DE41R(FDE41R,M,XN,YN,DYN,XK,H,Y,DY,DELTY,DF,RFN,RF,YP,DYP,IERR);
  Result := Result + Format(' %20.16f %20.16f %20.16f %20.16f %20.16f %20.16f ',
 [XK,Y[0],Y[1],H,Y1,Y2]) + #$0D#$0A;
  ХК := -5.0;
  Y1 := -Y1;
  Y2 := -Y2;
  H := 0.01;
_200:
 end;
UtRes('TDE41R',Result);  { вывод результатов в файл TDE41R.res }
exit;
end;

end.

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

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

implementation

procedure fde41r(T :Real; var Y :Array of Real; var Z :Array of Real;
                M :Integer);
begin
Z[0] := -6.e0*Y[0]-7.e0*Y[1];
Z[1] := -3.e0*Y[0]-2.e0*Y[1]+2.e0*T;
end;

end.

Результаты:

после первого обращения к подпрограмме -

                   XK                           Y(1)                             Y(2)
       5.000000000 + 00      5.674080692 + 01     - 5.635223804 + 01

                   H                             Y1                                 Y2
       1.000000000 - 02      5.674080540 + 01     - 5.635223633 + 01

после второго обращения к подпрограмме -

                   XK                           Y(1)                             Y(2)
     - 5.000000000 + 00    - 5.674080382 + 01       5.635223450 + 01

                   H                             Y1                                 Y2
     - 1.000000000 - 02    - 5.674080540 + 01       5.635223633 + 01